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

**Hugo Lebeau**

Université Grenoble Alpes, CNRS, Inria, Grenoble INP, LIG
Grenoble, 38000, France

HUGO.LEBEAU@UNIV-GRENOBLE-ALPES.FR

**Florent Chatelain**

Université Grenoble Alpes, CNRS, Grenoble INP, GIPSA-lab
Grenoble, 38000, France

FLORENT.CHATELAIN@GRENOBLE-INP.FR

**Romain Couillet**

Université Grenoble Alpes, CNRS, Inria, Grenoble INP, LIG
Grenoble, 38000, France

ROMAIN.COUILLET@UNIV-GRENOBLE-ALPES.FR

## 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., 2022), recommender systems (Karatzoglou et al., 2010; Rendle and Schmidt-Thieme, 2010; Frolov and Oseledets, 2017) and topic modeling (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 × pixel × wavelength × sample* in hyperspectral imaging (Zhang et al., 2013; Kanatsoulis et al., 2018), *time × spatial scale × electrode* in the EEG analysis by Acar et al.

©2025 Hugo Lebeau, Florent Chatelain and Romain Couillet.

License: CC-BY 4.0, see https://creativecommons.org/licenses/by/4.0/.
(2007) or neuron × time × 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 × time × patient
array, whereas Williams et al. (2018) use a low-rank tensor approximation on a neuron ×
time × 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 approxima-
tions 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 infor-
mation + 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 applica-
tions, Bi et al. (2021); 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 $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 $M$ in a compact way as the sum of $R$ rank-1 terms, $M = \sum_{i=1}^{R} \sigma_i u_i v_i^\top = U \Sigma V^\top$ where $U, V$ are respectively $n_1 \times R$ and $n_2 \times R$ matrices with orthonormal columns and $\Sigma$ is the $R \times R$ diagonal matrix of singular values. The rank of $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 tensor<fnref n="1" /> $\mathcal{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, $\mathcal{T} = \sum_{i=1}^{R} \sigma_i a_i \otimes b_i \otimes c_i$. What we have just described is the canonical polyadic decomposition (CPD) of $\mathcal{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 $a_i$ (resp. $b_i$, $c_i$) are, in general, no longer orthonormal. Conversely, retaining the orthonormality property inevitably results in the loss of the diagonality property, $\mathcal{T} = \sum_{i=1}^{r_1} \sum_{j=1}^{r_2} \sum_{k=1}^{r_3} G_{i,j,k} u_i \otimes v_j \otimes w_k$. This latter decomposition is called a Tucker decomposition and dates back to Tucker (1966). In fact, the best way to represent $\mathcal{T}$ with a Tucker decomposition is to choose the $u_i$ (resp. $v_i$, $w_i$) as the left singular vectors of the unfolding of $\mathcal{T}$ along mode 1 (resp. 2, 3)<fnref n="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) \le$

<footnote n="1">1. It is chosen of order 3 for simplicity of exposure.</footnote>

<footnote n="2">2. This is properly defined in Section 2.2.</footnote>
$R \le \min(r_1r_2, r_2r_3, r_1r_3)$. See, e.g., Sidiroopoulos 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 $\mathcal{T} \in \mathbb{R}^{n_1 \times \dots \times n_d}$ of possibly very high rank, we are interested in finding a low-rank approximation, i.e., an $n_1 \times \dots \times n_d$ tensor $\mathcal{X}$ which minimizes the distance $\|\mathcal{T} - \mathcal{X}\|_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_{\text{rank}(\mathcal{X}) \le (r_1, \dots, r_d)} \|\mathcal{T} - \mathcal{X}\|_F^2. \qquad (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{\mathcal{T}}$ remains a very good “first guess” as it verifies $\|\mathcal{T} - \mathcal{T}_*\|_F \le \|\mathcal{T} - \hat{\mathcal{T}}\|_F \le \sqrt{d}\|\mathcal{T} - \mathcal{T}_*\|_F$ where $\mathcal{T}_*$ 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 $\mathcal{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, $\mathcal{T} = \beta x^{\otimes d} + \mathcal{N}$ with $\|x\| = 1$, $\mathcal{N}$ Gaussian noise and $\beta$ a parameter controlling the signal-to-noise ratio (SNR). They show that estimation of $x$ from $\mathcal{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 is comprehensively studied by Benaych-Georges and Nadakuditi
(2012). Contrary to the matrix case however, Montanari and Richard (2014) make the dis-
turbing observation that none of the polynomial-time estimation algorithms among tensor
unfolding, power iteration and approximate message passing (AMP) succeed unless $\beta$ di-
verges as the dimensions of the tensor grow large. The results of Hopkins et al. (2015, 2017)
suggest that no polynomial-time algorithm can succeed unless $\beta \gtrsim N^{\frac{d-2}{4}}$, where $N$ scales
as the dimensions of the data tensor. While Perry et al. (2020) show that the information-
theoretic threshold is of order 1, this indicates 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 is studied by Ben Arous et al. (2019), who show 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. (2020) show<fnref n="3" /> 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$. Relying on the Kac-Rice method, Ros et al. (2019) thoroughly study such high-dimensional landscapes and classify the different behaviors and phase transitions which can occur.

So far, we have only referred to works dealing with the rank-1 symmetric case, but
there are also some studies on 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. The statistical inference of
finite rank tensors is studied by Chen et al. (2022) who identify the limit free energy
of the model in terms of a variational formula. Zhang and Xia (2018) consider a gen-
eral *low-multilinear-rank signal* $\mathcal{P}$ + *Gaussian noise* $\mathcal{N}$ model and bring to light the same
statistical-to-computational gap: if $\|\mathcal{P}\|_F$ is above a statistical threshold of order 1 then
Problem (1) has a solution which is aligned with the signal but is computationally in-
tractable unless $\|\mathcal{P}\|_F$ is above a computational threshold of order $N^{\frac{d-2}{4}}$. In this strong
SNR regime, the higher-order orthogonal iteration (HOOI) algorithm (De Lathauwer et al.,
2000a) is minimax-optimal. In fact, it is also proved by Ben Arous et al. (2020) that, with
Langevin dynamics and gradient descent, the algorithmic threshold behaves like $N^{\alpha}$ with
$\alpha > 0$. Other algorithmic thresholds have been shown as well for semi-definite and spec-
tral relaxations of the maximum likelihood problem (Hopkins et al., 2015, 2016; Kim et al.,
2017). AMP and tensor power iteration algorithms achieve $N^{\frac{d-1}{2}}$ (Lesieur et al., 2017;
Huang et al., 2022) while tensor unfolding methods (truncated MLSVD and HOOI algo-
rithm) achieve $N^{\frac{d-2}{4}}$ as already conjectured by Montanari and Richard (2014) and later
proven by Hopkins et al. (2015); Ben Arous et al. (2023) in the rank-1 case. The conver-
gence of the HOOI algorithm towards a local maximum for a *sufficiently close* initialization

<footnote n="3">3. The setting considered by Jagannath et al. (2020) is more general than the one of Ben Arous et al. (2019) because the noise in their model is not necessarily symmetric but the perturbation is still a rank-1 symmetric tensor.</footnote>
is proven by Xu (2018) and Feldman and Donoho (2023) show that it achieves exact recovery of a rank-1 perturbation in the large N regime when it is initialized with the dominant singular vectors of the unfoldings.

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) study the rank-1 symmetric case and are able to recover explicitly the same $\beta_s$ threshold as Jagannath et al. (2020) as well as to precisely quantify the alignment between a solution to Problem (1) and the signal. A similar analysis is 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). In particular, the results of Ben Arous et al. (2023) and Feldman (2023) on long random matrices, similar to those considered in this work, provide instructive insight into the recovery performance of tensor unfolding methods.

## 1.3 Summary of Contributions

In low-rank tensor approximation, tensor unfolding methods achieve the best known performance among polynomial-time algorithms. Motivated by several works suggesting that such method could actually reach the computational threshold (Hopkins et al., 2015, 2017; Zhang and Xia, 2018; Wein et al., 2019), we propose a thorough random matrix analysis of the low-multilinear-rank tensor approximation problem.

Consider the general spiked tensor model,

$$ \mathcal{T} = \mathcal{P} + \frac{1}{\sqrt{N}} \mathcal{N} \quad \in \mathbb{R}^{n_1 \times \dots \times n_d}, \quad \mathcal{N}_{i_1, \dots, i_d} \stackrel{\text{i.i.d.}}{\sim} \mathcal{N}(0, 1), \quad (2) $$

where $d \ge 3$ is the order of the tensor, $\mathcal{N}$ is an additive Gaussian noise, $N$ is a parameter controlling the size of the tensor such that the ratio $n_\ell/N$ is constant$^4$ (at least for $N$ above a certain threshold value $N_0$) for all $\ell \in \{1, \dots, d\}$ (for instance, $N = n_1$ or $N = \sum_{\ell=1}^d n_\ell$) and $\mathcal{P}$ is a low-multilinear-rank deterministic tensor, i.e., which can be decomposed as

$$ \mathcal{P} = \sum_{q_1=1}^{r_1} \dots \sum_{q_d=1}^{r_d} \mathcal{H}_{q_1, \dots, q_d} [\mathbf{x}_{q_1}^{(1)} \otimes \dots \otimes \mathbf{x}_{q_d}^{(d)}] \stackrel{\text{def}}{=} \left[ \mathcal{H}; \mathbf{X}^{(1)}, \dots, \mathbf{X}^{(d)} \right], \quad (3) $$

with $\mathcal{H} \in \mathbb{R}^{r_1 \times \dots \times r_d}$ and $\mathbf{X}^{(\ell)}$ an $n_\ell \times r_\ell$ matrix with orthonormal columns $\mathbf{x}_{q_\ell}^{(\ell)}$. The range of $\mathbf{X}^{(\ell)}$ is the $\ell$-th singular subspace of $\mathcal{P}$. This decomposition is illustrated for the case $d = 3$ in Figure 1. Model (2) with $\mathcal{P}$ as in Equation (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\}} \le R \le \min_{\ell\{\prod_{\ell' \ne \ell} r_{\ell'}\}}$ (Sidiropoulos et al., 2017).

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 $\mathcal{P}$ from $\mathcal{T}$ with a truncated MLSVD, which serves as initialization of the HOOI algorithm. In particular,

<footnote n="4">4. This ensures that the spectral norm of $\frac{1}{\sqrt{N}}\mathcal{N}$ is of order 1 (Tomioka and Suzuki, 2014).</footnote>
<figure><img src="image_1.png" /><figcaption>Figure 1: Illustration of the Tucker decomposition (3) of an $n_1 \times n_2 \times n_3$ tensor $\mathcal{P}$ with multilinear rank $(r_1, r_2, r_3)$. $\mathcal{H}$ is the $r_1 \times r_2 \times r_3$ core tensor and $\mathbf{X}^{(1)}, \mathbf{X}^{(2)}, \mathbf{X}^{(3)}$ are matrices with orthonormal columns spanning the singular subspaces of $\mathcal{P}$.</figcaption></figure>

we reveal that the interesting *non-trivial* regime is characterized by the $\Theta(N^{\frac{d-2}{2}})$ quantity<fnref n="5" />
$\sigma_N = \frac{1}{N} \prod_{\ell=1}^{d} \sqrt{n_\ell}$.

* If $\|\mathcal{P}\|_F^2 / \sigma_N \xrightarrow[N \to +\infty]{} 0$, then the noise completely masks the signal, and truncated MLSVD fails to recover $\mathcal{P}$.

* If $\|\mathcal{P}\|_F^2 / \sigma_N \xrightarrow[N \to +\infty]{} +\infty$, then the signal clearly stands out from the noise, and reconstruction of $\mathcal{P}$ with a truncated MLSVD is easy.

* If $\|\mathcal{P}\|_F^2/\sigma_N = \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 $\mathcal{P}$.

It is the analysis of this last regime which is of practical interest. Given the low-multilinear-rank approximation $\hat{\mathcal{T}} = [[\hat{\mathcal{G}}; \hat{\mathcal{U}}^{(1)}, \dots, \hat{\mathcal{U}}^{(d)}]]$ obtained with a truncated MLSVD of $\mathcal{T}$, we quantify how well $\hat{\mathcal{T}}$ reconstructs $\mathcal{P}$ in this non-trivial regime. To do so, we study the spectral properties of the unfoldings (i.e., matricizations) of the tensor $\mathcal{T}$, i.e., the $n_\ell \times \prod_{\ell' \neq \ell} n_{\ell'}$ matrices $\mathbf{T}^{(\ell)}$ whose columns are mode-$\ell$ fibers of $\mathcal{T}$ and the columns of $\hat{\mathcal{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. (2023) 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.

<footnote n="5">5. Note that this corresponds to the regime of the algorithmic threshold but, here, our measure of the signal power is $\|\mathcal{P}\|_F^2$ (and not $\|\mathcal{P}\|_F$) hence the $\Theta(N^{\frac{d-2}{2}})$ instead of $\Theta(N^{\frac{d-1}{2}})$.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 2: Alignments between singular subspaces (see Section 3.2) of the observation $\mathcal{T} = \sqrt{\omega} \mathcal{P}_o + \frac{1}{\sqrt{N}} \mathcal{N}$ and of the signal $\mathcal{P}_o$, with $\|\mathcal{P}_o\|_F^2 = \frac{\sqrt{n_1 n_2 n_3}}{N}$, as a function of the signal-to-noise ratio $\omega$. Theoretical alignments (Theorem 9) 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)$.</figcaption></figure>

Although the spectrum of $T^{(\ell)}T^{(\ell)\top}, \ell \in \{1, \dots, d\}$, diverges as $N \to +\infty$, we show that its eigenvalues (i.e., the squared singular values of $T^{(\ell)}$) gather in an interval $[\mu_N^{(\ell)} \pm 2\sigma_N]$ with $\mu_N^{(\ell)} = \frac{1}{N} \prod_{\ell' \neq \ell} n_{\ell'} = \Theta(N^{d-2})$. More precisely, the empirical spectral distribution of the centered-and-scaled matrix $\frac{1}{\sigma_N}[T^{(\ell)}T^{(\ell)\top} - \mu_N^{(\ell)}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 occurs: for each singular value of $P^{(\ell)}$ (the unfolding of $\mathcal{P}$ along mode $\ell$) which is above the threshold $\sqrt{\sigma_N}$, an eigenvalue of $T^{(\ell)}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 $T^{(\ell)}$) is aligned with the corresponding singular subspace of $\mathcal{P}$. The position of the isolated eigenvalue and this alignment are efficiently predicted by Theorem 9 (see also Figure 3).

As a result, Figure 2 plots, for an order-3 tensor, as a function of the signal-to-noise ratio (SNR) $\omega = \|\mathcal{P}\|_F^2/\sigma_N$, the alignments between the singular subspace of the signal $\mathcal{P}$ spanned by $\mathbf{X}^{(\ell)}$ and the dominant singular subspace of the observation $\mathcal{T}$ spanned by $\hat{\mathbf{U}}^{(\ell)}$. Solid curves are the alignments predicted by Theorem 9 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 $\mathcal{P}$ — the signal is masked by the noise. When it exceeds a critical value (see Theorem 9 and Section
3.2 for details), a phase transition phenomenon occurs$^6$: the alignment starts to grow —
i.e., truncated MLSVD now partially recovers $\mathcal{P}$ — and converges to 1 as $\omega \to +\infty$.

Besides, Figure 2 also plots the empirical alignments between the singular subspaces
of $\mathcal{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 15 that
the HOOI algorithm converges to a solution to Problem (2) as soon as its initialization suf-
ficiently 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 \|\mathcal{P}\|_F \ll N^{d/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 can be summarized as follows.

* We characterize, in the large $N$ limit, the behavior of the singular values of the unfoldings of the tensor $\mathcal{T}$ — denoted $T^{(\ell)}$, $\ell \in \{1, \dots, 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 $T^{(\ell)}T^{(\ell)\top}$ using standard tools from the theory of large random matrices.

* We give a precise condition, depending on a $\Theta(N^{d/4})$ threshold on the signal, for the detectability of a principal direction of the $\ell$-th singular subspace of $\mathcal{P}$ from the unfolding $T^{(\ell)}$. This corresponds to the presence of an isolated eigenvalue in the spectrum of $T^{(\ell)}T^{(\ell)\top}$ with associated eigenvector aligned with the sought singular subspace. We find similar formulae as Feldman (2023)$^7$ for the asymptotic position of this isolated eigenvalue, as well as the quality of the alignment (Theorem 9).

* Relying on our random matrix analysis, we characterize the performance of truncated MLSVD in the reconstruction of the signal $\mathcal{P}$ from the observation $\mathcal{T}$ (Section 3.2).

* We show that exact reconstruction of $\mathcal{P}$ from $\mathcal{T}$ is possible in the large $N$ regime with the HOOI algorithm (De Lathauwer et al., 2000a) as long as $\|\mathcal{P}\|_F \gg 1$ and it is initialized in the right basin of attraction (Theorem 15). Without prior information, this depends on the detectability of $\mathcal{P}$ in the truncated MLSVD of $\mathcal{T}$, which is only possible above the $\Theta(N^{d/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

<footnote n="6">6. In fact, we will see in Section 3.2 that there is one phase transition for each principal direction of the singular subspaces of $\mathcal{P}$, resulting in $\sum_{\ell=1}^d r_\ell$ phase transitions. Their positions correspond to sudden changes of slope in the solid curves of Figure 2.</footnote>

<footnote n="7">7. Feldman (2023) studies the spiked model associated with long random matrices within a similar framework as Benaych-Georges and Nadakuditi (2011). Our approach is different in that it mostly relies on a *deterministic equivalent* which is introduced in Theorem 5.</footnote>
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, **a**, **A** and **A** respectively denote a scalar, a vector, a matrix and a tensor. $a_i$, $A_{i,j}$ and $A_{i_1, \dots, 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, \dots, 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 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 \ge 0}$, its convergence in distribution to $\mathcal{L}$ is denoted $X_n \xrightarrow[n \to +\infty]{} \mathcal{L}$ and its almost sure convergence to $L$ is denoted $X_n \xrightarrow[n \to +\infty]{} \mathcal{L}$. The normal distribution with mean $\mu$ and variance $\sigma^2$ is denoted $N(\mu, \sigma^2)$. The span of an $n_1 \times n_2$ matrix **A** is Span **A** = $\{\mathbf{A}x \mid x \in \mathbb{R}^{n_2}\} \subset \mathbb{R}^{n_1}$. The singular values of **A** in non-increasing order are denoted $s_1(\mathbf{A}) \ge s_2(\mathbf{A}) \ge \dots \ge 0$. Given an $n \times n$ matrix **B**, its trace is $\mathrm{Tr} \, \mathbf{B} = \sum_{i=1}^n B_{i,i}$ and its spectrum, $\mathrm{Sp} \, \mathbf{B}$, is the set of all its eigenvalues. $\|\cdot\|$ denotes the standard Euclidean norm for vectors and the corresponding operator norm (spectral norm) for matrices.

Given two real sequences $(u_n)_{n \ge 0}$ and $(v_n)_{n \ge 0}$, we write $u_n = \mathcal{O}_{n \to +\infty}(v_n)$ if there exists a constant $C > 0$ and an integer $n_0$ such that $|u_n| \le C|v_n|$ as soon as $n \ge n_0$. If $u_n = \mathcal{O}_{n \to +\infty}(v_n)$ and $v_n = \mathcal{O}_{n \to +\infty}(u_n)^8$, then we write $u_n = \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 $|u_n| \le \varepsilon|v_n|$ as soon as $n \ge n_0$. When it is clear from context that $n \to +\infty$, we simply write $u_n = \Theta(v_n)$, $u_n = \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: $\mathcal{T} \in \mathbb{R}^{n_1 \times \cdots \times n_d}$ is a collection of elements $\mathcal{T}_{i_1, \ldots, i_d} \in \mathbb{R}$ with $i_\ell \in [n_\ell]$, $\ell \in [d]$. An $\ell$-fiber of $\mathcal{T}$ is the vector of $\mathbb{R}^{n_\ell}$ obtained by fixing all the indices of $\mathcal{T}$ but the $\ell$-th. This is the generalization of columns and rows of matrices, which are respectively 1- and 2-fibers. Unfolding (or matricization) is the process

<footnote n="8">8. I.e., if there exist two constants $c, C > 0$ and an integer $n_0$ such that $c|v_n| \le |u_n| \le C|v_n|$ as soon as $n \ge n_0$.</footnote>
by which a matrix is built from a tensor — $\mathcal{T}^{(\ell)}$ is the unfolding of $\mathcal{T}$ along mode $\ell$, i.e., the $n_\ell \times \prod_{\ell' \neq \ell} n_{\ell'}$ matrix whose columns are $\ell$-fibers of $\mathcal{T}^9$.

Given two $n_1 \times \dots \times n_d$ tensors $\mathcal{A}, \mathcal{B}$, their scalar product is

$$ \langle \mathcal{A}, \mathcal{B} \rangle_{\mathrm{F}} \stackrel{\mathrm{def}}{=} \sum_{i_1=1}^{n_1} \dots \sum_{i_d=1}^{n_d} \mathcal{A}_{i_1, \dots, i_d} \mathcal{B}_{i_1, \dots, i_d}. $$

The Frobenius norm of $\mathcal{A}$ is $\|\mathcal{A}\|_{\mathrm{F}} \stackrel{\mathrm{def}}{=} \sqrt{\langle \mathcal{A}, \mathcal{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, $x y^\top \stackrel{\mathrm{def}}{=} x \otimes 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., $\mathcal{T}_{i_1, \dots, i_d} = x_{i_1}^{(1)} \dots x_{i_d}^{(d)} \iff \mathcal{T} \stackrel{\mathrm{def}}{=} \otimes_{\ell=1}^d x^{(\ell)}$.

Given two matrices $\mathcal{A}$ and $\mathcal{B}$ of respective sizes $n_1 \times n_2$ and $p_1 \times p_2$, their Kronecker product, denoted $\mathcal{A} \boxtimes \mathcal{B}$, is the $n_1 p_1 \times n_2 p_2$ matrix such that $[\mathcal{A} \boxtimes \mathcal{B}]_{p_1(i_1-1)+j_1, p_2(i_2-1)+j_2} = \mathcal{A}_{i_1, i_2} \mathcal{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,

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

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

Because of their multidimensional 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).** $\mathcal{T}$ is decomposed as a sum of $R$ rank-1 terms, $\mathcal{T} = \sum_{i=1}^R \sigma_i \otimes_{\ell=1}^d a_i^{(\ell)}$, with $a_i^{(\ell)}$ 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).** $\mathcal{T}$ is written as the multilinear tensor-matrix product of an $r_1 \times \dots \times r_d$ core tensor $\mathcal{G}$ with $d$ factors $U^{(\ell)}, \ell \in [d]$, of respective sizes $n_\ell \times r_\ell$ and orthonormal columns $u_{q_\ell}^{(\ell)}$,

$$ \mathcal{T} = \sum_{q_1=1}^{r_1} \dots \sum_{q_d=1}^{r_d} \mathcal{G}_{q_1, \dots, q_d} \otimes_{\ell=1}^{d} u_{q_\ell}^{(\ell)} \stackrel{\mathrm{def}}{=} \left[ \mathcal{G}; U^{(1)}, \dots, U^{(d)} \right]. $$

The integer $r_\ell$ is the dimension of the subspace spanned by the $\ell$-fibers of $\mathcal{T}$, whose an orthonormal basis is formed by the columns of $U^{(\ell)}$. This decomposition, introduced by Tucker (1966), is also called higher-order SVD (HOSVD, De Lathauwer et al., 2000b). From

<footnote n="9">9. The order of the columns does not matter, as long as other operations, such as the Kronecker product, are defined in a consistent manner.</footnote>
this decomposition, the unfoldings of $\mathcal{T}$ can be expressed with matrix and Kronecker products between the factor matrices and the unfoldings of the core tensor $\mathcal{G}$. For instance, if $\mathcal{T} = [[\mathcal{G}; U, V, W]]$,

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

Note how these expressions generalize (4).

The contraction of $\mathcal{T}$ on $A^{(1)}, \dots, A^{(d)}$ with $A^{(\ell)} \in \mathbb{R}^{n_\ell \times p_\ell}$, $\ell \in [d]$, is the $p_1 \times \dots \times p_d$ tensor $\mathcal{T}(A^{(1)}, \dots, A^{(d)})$ whose $(j_1, \dots, j_d)$-entry is $\sum_{i_1=1}^{n_1} \dots \sum_{i_d=1}^{n_d} \mathcal{T}_{i_1, \dots, i_d} \prod_{\ell=1}^d A^{(\ell)}_{i_\ell, j_\ell}$. If $\mathcal{T} = [[\mathcal{G}; U^{(1)}, \dots, U^{(d)}]]$ then $\mathcal{G} = [[\mathcal{T}; U^{(1)\top}, \dots, U^{(d)\top}]] = \mathcal{T}(U^{(1)}, \dots, U^{(d)})$. The tensor contraction generalizes the matrix operation $A^{(1)\top}MA^{(2)} = M(A^{(1)}, A^{(2)})$.

**Remark 1 (Uniqueness of the MLSVD up to isometries)** For all $r_\ell \times r_\ell$ orthogonal matrices $O^{(\ell)}$, $\ell \in [d]$, we have the equivalent decomposition

$$[[\mathcal{G}; U^{(1)}, \dots, U^{(d)}]] = [[\check{\mathcal{G}}; U^{(1)}O^{(1)}, \dots, U^{(d)}O^{(d)}]]$$

where $\check{\mathcal{G}} = \mathcal{G}(O^{(1)}, \dots, O^{(d)})$ is the contraction of $\mathcal{G}$ on $O^{(1)}, \dots, O^{(d)}$. Nevertheless, up to isometries, the multilinear singular value decomposition is unique (De Lathauwer et al., 2000b, Property 4).

The actual rank of a tensor is usually understood as its CPD-rank $R$ whereas its MLSVD-rank $(r_1, \dots, r_d)$ is called *multilinear rank*. These two ranks are constrained by the following inequalities,

$$\max_{1 \le l \le d} r_l \le R \le \min_{1 \le l \le d} \prod_{l' \ne l} r_{l'}, \quad (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 $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 $M$, i.e., $\frac{1}{n} \sum_{\lambda \in \text{Sp } 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 $M$ is symmetric, its eigenvalues and eigenvectors are real and a key tool to answer the preceding questions is its resolvent, $Q_M(z) = (M - zI_n)^{-1}$, defined for all $z \in \mathbb{C} \setminus \text{Sp } M$. Indeed, $\frac{1}{n} \text{Tr } Q_M(z) = \frac{1}{n} \sum_{\lambda \in \text{Sp } 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 \text{Supp } \mu$ as $m_\mu(z) = \int_{\mathbb{R}} \frac{1}{1-z} d\mu(t)$.
The knowledge of $m_μ$ is equivalent to the knowledge of $μ$ thanks to the inversion formula $μ([a, b]) = \frac{1}{\pi} \lim_{η \downarrow 0} \int_a^b \mathfrak{F}[m_μ(x+iη)] dx$ for any continuity points $a, b$ of $μ$ (Couillet and Liao, 2022, Theorem 2.1). Hence, following the behavior of $\frac{1}{n} \text{Tr} \mathcal{Q}_M(z)$ as $n \to +\infty$ gives insight into the limiting spectral distribution of $M$ via its Stieltjes transform. To do so, we seek first a deterministic equivalent of $\mathcal{Q}_M$.

**Definition 3 (Deterministic equivalent)** Let $X$ be a random $n \times n$ matrix and $\bar{X}$ be a deterministic $n \times n$ matrix. We write $X \leftrightarrow \bar{X}$ if, for all deterministic matrices $A \in \mathbb{R}^{n \times n}$ and vectors $a, b \in \mathbb{R}^n$ of bounded norms (spectral and Euclidean norms respectively)$^{10}$,

$$ \frac{1}{n} \text{Tr} A(X - \bar{X}) \xrightarrow{n \to +\infty} 0 \quad \text{and} \quad a^\top(X - \bar{X})b \xrightarrow{n \to +\infty} 0. $$

The matrix $Y$ is called a deterministic equivalent of $X$. For more details on this notion, see Couillet and Liao (2022, §2.1.4).

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{C}$ be a continuously differentiable function. When the following expectations exist, $\mathbb{E}[Zf(Z)] = \mathbb{E}[f'(Z)]$.

Once a deterministic equivalent $\bar{Q}_M$ of $\mathcal{Q}_M$ is found, the limiting spectral distribution of $M$ is accessible through the limit of $\frac{1}{n} \text{Tr} \bar{Q}_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 $\mathcal{Q}_M$. Hence, the asymptotic analysis of the resolvent provides equations yielding their almost sure asymptotic position. Moreover, thanks to Cauchy's integral formula, for all $a \in \mathbb{R}^n$, $|a^\top u|^2 = \frac{1}{2i\pi} \oint_a^b a^\top \mathcal{Q}_M(z)a dz$ where $\gamma$ is a positively-oriented simple closed complex contour circling around an isolated eigenvalue (assuming it has multiplicity 1) with all other eigenvalues left outside and $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, see Appendix A for an extension. 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 Equation (2) using tools presented in Section 2.3. We give precise results on the

<footnote n="10">10. Strictly 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.</footnote>
spectral behavior of the unfoldings of the observed tensor $\mathcal{T}$, and specify the achievable performance in the estimation of the underlying signal $\mathcal{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 \dots \times n_d$ tensor $\mathcal{T} = \mathcal{P} + \frac{1}{\sqrt{N}}\mathcal{N}$ of order $d \ge 3$, modeling a low-rank deterministic signal $\mathcal{P}$ corrupted by an additive Gaussian noise tensor $\mathcal{N}$ whose entries are independent $\mathcal{N}(0, 1)$ random variables$^{11}$. We denote by $(r_1, \dots, r_d)$ the multilinear rank of $\mathcal{P}$ and study this model in the asymptotic regime where $N \to +\infty$ with $n_\ell = \Theta(N)$ and $r_\ell = \Theta(1)$, $\ell \in [d]$.

The estimation of $\mathcal{P}$ with a truncated MLSVD on $\mathcal{T}$ is simply the computation of the dominant singular subspaces of $\mathcal{T}$. Specifically, $\hat{U}^{(\ell)} \in \mathbb{R}^{n_\ell \times r_\ell}$ gathers the $r_\ell$ dominant left singular vectors of $T^{(\ell)}$ — and thus, $\hat{U}^{(\ell)\top}\hat{U}^{(\ell)} = I_{r_\ell}$. Then, a low-multilinear-rank approximation of $\mathcal{T}$ is $\hat{\mathcal{T}} = [[\hat{\mathcal{G}}; \hat{U}^{(1)}, \dots, \hat{U}^{(d)}]]$ with an $r_1 \times \dots \times r_d$ core tensor $\hat{\mathcal{G}} = \mathcal{T}(\hat{U}^{(1)}, \dots, \hat{U}^{(d)})$. An equivalent expression is $\hat{\mathcal{T}} = [[\mathcal{T}; \hat{U}^{(1)}\hat{U}^{(1)\top}, \dots, \hat{U}^{(d)}\hat{U}^{(d)\top}]]$, which explicitly shows that $\hat{\mathcal{T}}$ is the projection of $\mathcal{T}$ on its dominant singular subspaces$^{12}$. Thus, the quality of this estimation hinges upon the alignments between the singular subspaces of $\mathcal{P}$ and the dominant singular subspaces of $\mathcal{T}$. Namely, denoting $\hat{u}_{q_\ell}^{(\ell)}$, for $q_\ell \in [r_\ell]$, the columns of $\hat{U}^{(\ell)} = [\hat{u}_1^{(\ell)} \dots \hat{u}_{r_\ell}^{(\ell)}]$ and given that $\mathcal{P} = [[\mathcal{H}; X^{(1)}, \dots, X^{(d)}]]$, the quantities of interest are $\|X^{(\ell)\top}\hat{u}_{q_\ell}^{(\ell)}\|^2$ since they represent how much of $\hat{u}_{q_\ell}^{(\ell)}$ is in the $\ell$-th singular subspace of the signal $\mathcal{P}$.

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

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

<footnote n="11">11. We highlight the fact that the Gaussian noise assumption is not restrictive. Firstly, the universality result of Gurau (2014) shows that, as $N \to +\infty$, the distribution of a random tensor with i.i.d. entries has the same limit than that of a tensor with i.i.d. Gaussian entries. Moreover, our results can be extended to non-Gaussian noise up to a control on the moments of the distribution with the “interpolation trick” of Lytova and Pastur (2009, Corollary 3.1). These are technical aspects which go beyond the scope of our work. Examples of results relying on this method are Pastur and Shcherbina (2011, Theorems 18.4.2 and 19.2.1); Merlevède et al. (2015); Banna et al. (2015).</footnote>

<footnote n="12">12. In the more familiar matrix case ($d = 2$), the expression $[[T; \hat{U}^{(1)}\hat{U}^{(1)\top}, \hat{U}^{(2)}\hat{U}^{(2)\top}]]$ is equivalent to $\hat{U}^{(1)}\hat{U}^{(1)\top}T\hat{U}^{(2)}\hat{U}^{(2)\top}$.</footnote>
Hence, we need to consider instead a centered-and-scaled version of our random matrix $T^{(\ell)}T^{(\ell)\top}$ to properly study the behavior of its spectrum. The quantities $\mu_N^{(\ell)}$ and $\sigma_N$ introduced in Theorem 5 below are such that the eigenvalues of $\frac{1}{\sigma_N}[T^{(\ell)}T^{(\ell)\top} - \mu_N^{(\ell)}I_{n_\ell}]$ neither diverge nor vanish but stay at a $\Theta(1)$ scale as $N \to +\infty$.

**Theorem 5 (Deterministic equivalent)** For $\ell \in [d]$, define the following quantities,

$$ \mu_N^{(\ell)} = \frac{1}{N} \prod_{\ell \neq \ell'} n_{\ell'}, \qquad \sigma_N = \frac{1}{N} \sqrt{\prod_{\ell \in [d]} n_{\ell}}. $$

As $N \to +\infty$, if the ratio $\|\mathcal{P}\|_F^2/\sigma_N$ is bounded, then the resolvent of the centered-and-scaled matrix $\frac{1}{\sigma_N}[T^{(\ell)}T^{(\ell)\top} - \mu_N^{(\ell)}I_{n_\ell}]$ has the following deterministic equivalent (Definition 3), defined for all $\tilde{z} \in \mathbb{C} \setminus \mathbb{R}$,

$$ \tilde{Q}^{(\ell)}(\tilde{z}) \stackrel{\text{def}}{=} \left( \frac{1}{\sigma_N} [T^{(\ell)}T^{(\ell)\top} - \mu_N^{(\ell)}I_{n_\ell}] - \tilde{z}I_{n_\ell} \right)^{-1} \\ \leftarrow \tilde{Q}^{(\ell)}(\tilde{z}) \stackrel{\text{def}}{=} \left( \frac{1}{\sigma_N} P^{(\ell)}P^{(\ell)\top} + \frac{1}{\tilde{m}(\tilde{z})} I_{n_\ell} \right)^{-1} $$

where $\tilde{m}(\tilde{z}) \stackrel{\text{def}}{=} \lim_{N \to +\infty} \frac{1}{n_\ell} \text{Tr} \tilde{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. \qquad (6) $$

**Proof** See Appendix B. ■

This theorem is fundamental. It gives a deterministic equivalent $\tilde{Q}^{(\ell)}(\tilde{z})$ of the resolvent of the centered-and-scaled matrix $\frac{1}{\sigma_N}[T^{(\ell)}T^{(\ell)\top} - \mu_N^{(\ell)}I_{n_\ell}]$, 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_N^{(\ell)}$ and $\sigma_N$ are respectively $\Theta(N^{d-2})$ and $\Theta(N^{\frac{d-2}{2}})$, meaning that the eigenvalues of $T^{(\ell)}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 Equation (6) characterizes $\tilde{m}$, the Stieltjes transform of the limiting spectral distribution (LSD) of $\frac{1}{\sigma_N}[T^{(\ell)}T^{(\ell)\top} - \mu_N^{(\ell)}I_{n_\ell}]$. We see that this LSD is the same regardless of the low-rank perturbation $P^{(\ell)}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 $P^{(\ell)} = 0_{n_1 \times \dots \times n_d}$, we simply have $\tilde{Q}^{(\ell)}(\tilde{z}) = \tilde{m}(\tilde{z})I_{n_\ell}$. As a corollary, we recover the limiting spectral distribution of long random matrices, which was first characterized by Bai and Yin (1988).

**Corollary 6 (Limiting spectral distribution)** As $N \to +\infty$, the empirical spectral distribution of the centered-and-scaled matrix $\frac{1}{\sigma_N}[T^{(\ell)}T^{(\ell)\top} - \mu_N^{(\ell)}I_{n_\ell}]$ converges weakly almost surely to $\mu_{\text{SC}}$, the semicircle distribution on $[-2, +2]$,

$$ d\mu_{\text{SC}}(x) \stackrel{\text{def}}{=} \frac{1}{2\pi} \sqrt{[4 - x^2]^+} dx. $$
<figure><img src="image_3.png" /><figcaption>Figure 3: **Top:** empirical spectral distribution (ESD) of $T^{(\ell)}T^{(\ell)\top}$. The orange curve is the density of the stretched semicircle on $[\mu_N^{(\ell)} \pm 2\sigma_N]$ (Corollary 6). Green dashed lines represent asymptotic positions of spikes $\mu_N^{(\ell)} + \sigma_N \tilde{\xi}_{q_\ell}^{(\ell)}$ (Theorem 9). **Bottom:** Observed alignments between the dominant eigenvectors of $T^{(\ell)}T^{(\ell)\top}$ and $P^{(\ell)}P^{(\ell)\top}$ (purple bars) with their predicted asymptotic values $[\zeta_{q_\ell}^{(\ell)}]^+$ (red curve, Theorem 9). **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 $\|P\|_F^2/\sigma_N = 15$.</figcaption></figure>

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

This result states that the limiting spectral distribution of $\frac{1}{\sigma_N}[T^{(\ell)}T^{(\ell)\top} - \mu_N^{(\ell)}I_{n_\ell}]$ 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 $T^{(\ell)}T^{(\ell)\top}$ is a stretched semicircle on $[\mu_N^{(\ell)} \pm 2\sigma_N]$. This phenomenon is illustrated in the first row of Figure 3, where the empirical spectral distribution (ESD) of $T^{(\ell)}T^{(\ell)\top}$ is represented with the corresponding stretched semicircle for every mode $\ell$ of an order-3 tensor of size 300 × 500 × 700 following the general spiked tensor model (2).

**Remark 7 (From Marčenko-Pastur to Wigner)** Given a random matrix $X \in \mathbb{R}^{p_1 \times p_2}$ with i.i.d. $N(0, \frac{1}{p_2})$ entries, it is well known that the empirical spectral distribution of $XX^\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_{SC}$ is known to be the limiting spectral distribution of symmetric $p \times p$ random matrices with i.i.d. (up to symmetry) $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, Corollary 6 shows that if $p_2$ grows at a faster polynomial rate than $p_1$, the matrix $X X^\top$ behaves asymptotically (up to a deterministic rescaling and shift) 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 = \Theta(N)$), e.g., $(n_1, n_2, n_3) = (1000, 40, 40)$, then the empirical spectral distribution of $T^{(1)}T^{(1)\top}$ is better modeled by a Marčenko-Pastur distribution than by a Wigner semicircle.

**Remark 8 (Confinement of the spectrum)** The weak convergence of the empirical spectral distribution to $\mu_{SC}$ stated in Corollary 6 could allow for a negligible amount of eigenvalues to stay outside the support of the limiting spectral distribution. In fact, in Appendix B, we show an even more precise statement on the global behavior of the eigenvalues of $\frac{1}{N}N^{(\ell)}N^{(\ell)\top}$ (the model without signal): for all $\varepsilon > 0$, there exists an integer $N_0$ such that $\text{Dist}(\frac{1}{\sigma_N}[\lambda - \mu_N^{(\ell)}], [-2, 2]) \le \varepsilon$ almost surely for all $\lambda \in \text{Sp} \frac{1}{N}N^{(\ell)}N^{(\ell)\top}$ as soon as $N \ge N_0$. This means that no eigenvalue of $\frac{1}{\sigma_N}[\frac{1}{N}N^{(\ell)}N^{(\ell)\top} - \mu_N^{(\ell)}I_{n_\ell}]$ stays outside the support of the semicircle distribution $[-2, 2]$ almost surely as $N \to +\infty$.

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 $\mathcal{P}$ which, in this setting, has multilinear rank (3, 4, 5). The estimate of $\mathcal{P}$ given by a truncated MLSVD on $\mathcal{T}$ has its singular subspaces spanned by the dominant eigenvectors of $T^{(\ell)}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 9 (Spike behavior)** For $\ell \in [d]$ and $q_\ell \in [r_\ell]$, define the following quantities$^{13}$,

$$ \rho_{q_\ell}^{(\ell)} = \frac{s_{q_\ell}^2(\mathbf{P}^{(\ell)})}{\sigma_N}, \quad \tilde{\xi}_{q_\ell}^{(\ell)} = \rho_{q_\ell}^{(\ell)} + \frac{1}{\rho_{q_\ell}^{(\ell)}} \quad \text{and} \quad \zeta_{q_\ell}^{(\ell)} = 1 - \frac{1}{[\rho_{q_\ell}^{(\ell)}]^2}. $$

As $N \to +\infty$, if the ratio $\|\mathcal{P}\|_F^2/\sigma_N$ is bounded and $\rho_{q_\ell}^{(\ell)} > 1$, then

$$ \frac{1}{\sigma_N} [s_{q_\ell}^2(\mathbf{T}^{(\ell)}) - \mu_N^{(\ell)}] \xrightarrow{\text{a.s.}} \tilde{\xi}_{q_\ell}^{(\ell)} \quad \text{and} \quad \|\mathbf{X}^{(\ell)\top}\hat{\mathbf{u}}_{q_\ell}^{(\ell)}\|^2 \xrightarrow{\text{a.s.}} \zeta_{q_\ell}^{(\ell)} $$

where $\hat{\mathbf{u}}_{q_\ell}^{(\ell)}$ is the $q_\ell$-th dominant left singular vector of $\mathbf{T}^{(\ell)}$.

**Proof** See Appendix C.

<footnote n="13">13. We recall that the notation $s_i(\mathbf{A})$ denotes the $i$-th singular value of $\mathbf{A}$ in non-increasing order.</footnote>
The first quantity defined in this theorem, $\rho_{q\ell}^{(\ell)}$, should be understood as a signal-to-noise ratio (SNR). Indeed, the squared $q\ell$-th singular value of $\mathbf{P}^{(\ell)}$ (i.e., the $q\ell$-th eigenvalue of $\mathbf{P}^{(\ell)}\mathbf{P}^{(\ell)\top}$), $s_{q\ell}^2(\mathbf{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}_{q\ell}^{(\ell)}$ and $\zeta_{q\ell}^{(\ell)}$ depend only on the value of this SNR and indicate respectively the position of an isolated eigenvalue in the spectrum of $\mathbf{T}^{(\ell)}\mathbf{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_{q\ell}^{(\ell)} > 1$, an eigenvalue of $\mathbf{T}^{(\ell)}\mathbf{T}^{(\ell)\top}$ isolates itself from the semicircle$^{14}$ around $\mu_N^{(\ell)} + \sigma_N \tilde{\xi}_{q\ell}^{(\ell)}$. Moreover, recalling that $\mathcal{P} = [\mathcal{H}; \mathbf{X}^{(1)}, \dots, \mathbf{X}^{(d)}]$, the eigenvector associated with this isolated eigenvalue is aligned with the subspace spanned by $\mathbf{X}^{(\ell)}$, which is the $\ell$-th singular subspace of $\mathcal{P}$. The quality of this alignment is given by $0 < \zeta_{q\ell}^{(\ell)} \le 1$.

Most importantly, this result reveals the non-trivial regime for the estimation of $\mathcal{P}$ with a truncated MLSVD. Since $\sigma_N = \Theta(N^{\frac{d-2}{2}})$, it shows that $\|\mathcal{P}\|_F^2 = \sum_{q\ell=1}^{r_\ell} s_{q\ell}^2(\mathbf{P}^{(\ell)})$ must also be of the same order. Indeed, if $\|\mathcal{P}\|_F^2 \ll N^{\frac{d-2}{2}}$, then $\rho_{q\ell}^{(\ell)} \to 0$, the SNR is too small and no signal can be recovered, whereas if $\|\mathcal{P}\|_F^2 \gg N^{\frac{d-2}{2}}$, then $\rho_{q\ell}^{(\ell)} \to +\infty$, the SNR is very high and recovery of $\mathcal{P}$ is easy. It is precisely between these two regimes, i.e., $\|\mathcal{P}\|_F = \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. (2023) and Zhang and Xia (2018). In this non-trivial regime, the quantities $\zeta_{q\ell}^{(\ell)}$ given in Theorem 9 precisely quantify how well the dominant eigenvectors of $\mathbf{T}^{(\ell)}\mathbf{T}^{(\ell)\top}$ are aligned with the sought signal, i.e., the singular subspaces of $\mathcal{P}$. In section 3.2 below, this result is used to study the reconstruction performance of truncated MLSVD.

**Remark 10 (Spiked Wigner model)** The reader familiar with the spiked Wigner model may have recognized the expressions of $\tilde{\xi}_{q\ell}^{(\ell)}$ and $\zeta_{q\ell}^{(\ell)}$ given in Theorem 9. Indeed, given a symmetric $p \times p$ random matrix $\mathbf{W}$ with i.i.d. (up to symmetry) $N(0, \frac{1}{p})$ entries, the spectrum of $[\rho\mathbf{x}\mathbf{x}^\top + \mathbf{W}]$ with $\|\mathbf{x}\| = 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 $\mathbf{u}$ is such that $|\mathbf{x}^\top\mathbf{u}|^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 a deterministic rescaling and shift, $\mathbf{T}^{(\ell)}\mathbf{T}^{(\ell)\top}$ asymptotically behaves like a spiked Wigner matrix.

Theorem 9 is illustrated in Figure 3. In the first row, asymptotic positions of isolated eigenvalues $\mu_N^{(\ell)} + \sigma_N \tilde{\xi}_{q\ell}^{(\ell)}$ are represented by the green dashed lines. In our experiment, $\mathcal{P}$ has multilinear rank (3, 4, 5). Hence 3, 4 and 5 isolated eigenvalues are expected in the spectrum of $\mathbf{T}^{(\ell)}\mathbf{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(\mathbf{P}^{(3)})$ is not “energetic enough” to extricate itself from the bulk of eigenvalues, i.e, the SNR $\rho_5^{(3)}$ is below the phase transition threshold. Hence the 5-th dominant left singular vector of $\mathbf{T}^{(\ell)}$

<footnote n="14">14. Indeed, note that $\rho_{q\ell}^{(\ell)} > 1 \Rightarrow \tilde{\xi}_{q\ell}^{(\ell)} > 2$.</footnote>
is not informative as it is not aligned with the 3-rd singular subspace of $\mathcal{P}$, spanned by $\mathbf{X}^{(3)}$. The second row of Figure 3 depicts the alignments of the spiked eigenvectors $\hat{\mathbf{u}}_{q\ell}^{(\ell)}$ with the corresponding singular subspaces of $\mathcal{P}$ as well as the asymptotic alignment given by Theorem 9 as a function of the position of the associated eigenvalue. It appears that the higher is the SNR $\rho_{q\ell}^{(\ell)}$, the farther is the isolated eigenvalue from the bulk and the more is the corresponding eigenvector aligned with the span of $\mathbf{X}^{(\ell)}$. This assertion can be intuitively understood in terms of “energy” $s_{q\ell}^2(\mathbf{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 $\mathcal{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 $\mathcal{T}$ following the general spiked tensor model (2), we consider its low-rank approximation $\hat{\mathcal{T}} = [[\hat{\mathcal{G}}; \hat{\mathcal{U}}^{(1)}, \dots, \hat{\mathcal{U}}^{(d)}]]$ where $\hat{\mathcal{U}}^{(\ell)}$ is the $n_\ell \times r_\ell$ matrix whose columns are the $r_\ell$ dominant singular vectors of $\mathcal{T}$ and $\hat{\mathcal{G}} = \mathcal{T}(\hat{\mathcal{U}}^{(1)}, \dots, \hat{\mathcal{U}}^{(d)})$. $\hat{\mathcal{T}}$ is the projection of $\mathcal{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{\mathcal{T}}$ is $\mathcal{P} = [[\mathcal{H}; \mathbf{X}^{(1)}, \dots, \mathbf{X}^{(d)}]]$. The reconstruction performance of $\hat{\mathcal{T}}$ hence depends on how well the subspace$^{15}$ spanned by $\hat{\mathcal{U}}^{(\ell)}$ estimates the one spanned by $\mathbf{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, $\mathrm{Span} \mathbf{X}^{(\ell)}$ and $\mathrm{Span} \hat{\mathcal{U}}^{(\ell)}$), one can define a set of mutual angles which are invariant under isometric transformation.

**Definition 11 (Principal angles)** The principal angles $\theta_{q\ell}^{(\ell)} \in [0, \frac{\pi}{2}]$ between the subspaces spanned by $\mathbf{X}^{(\ell)}$ and $\hat{\mathcal{U}}^{(\ell)}$ are recursively defined for $q_\ell = 1, \dots, r_\ell$ by

$$ \cos \theta_{q\ell}^{(\ell)} = x_{q\ell}^{\top} u_{q\ell} \quad \text{with} \quad (x_{q\ell}, u_{q\ell}) \in \underset{(x, u) \in \mathrm{Span} \mathbf{X}^{(\ell)} \times \mathrm{Span} \hat{\mathcal{U}}^{(\ell)}}{\arg \max} x^{\top} u. $$
$$ x^{\top} x_{q\ell} = 0, \quad u^{\top} u_{q\ell} = 0, \quad 1 \le q_{\ell} < q_{\ell-1} $$

Moreover, we have the following useful property.

**Proposition 12 (Björck and Golub, 1973)** The $q_\ell$-th singular value of $\mathbf{X}^{(\ell)\top}\hat{\mathcal{U}}^{(\ell)}$ in non-increasing order equals the cosine of the $q_\ell$-th principal angle,

$$ s_{q\ell}(\mathbf{X}^{(\ell)\top}\hat{\mathcal{U}}^{(\ell)}) = \cos \theta_{q\ell}^{(\ell)}, \quad \ell \in [d], \quad q\ell \in [r_\ell]. $$

In fact, information about the alignment between the subspaces induced by $\mathbf{X}^{(\ell)}$ and $\hat{\mathcal{U}}^{(\ell)}$ are contained entirely in the $r_\ell \times r_\ell$ matrix $\mathbf{X}^{(\ell)\top}\hat{\mathcal{U}}^{(\ell)}$. Following Definition 11, we know from Theorem 9 that, as $N \to +\infty$, $\cos^2 \theta_{q\ell}^{(\ell)} \to [\zeta_{q\ell}^{(\ell)}]^+$ almost surely$^{16}$, where we

<footnote n="15">15. The object of importance is indeed the subspace and not the matrix $\hat{\mathcal{U}}^{(\ell)}$ since any other matrix $\hat{\mathcal{U}}^{(\ell)}\mathcal{O}^{(\ell)}$, with $\mathcal{O}^{(\ell)}$ an $r_\ell \times r_\ell$ orthogonal matrix, would span the same subspace and therefore give the same approximation.</footnote>

<footnote n="16">16. Indeed, since $\zeta_1^{(\ell)} \ge \dots \ge \zeta_{r_\ell}^{(\ell)}$, observe that, in Definition 11, $|x_{q\ell}^{\top} u_{q\ell}|^2$ is asymptotically bounded by $[\zeta_{q\ell}^{(\ell)}]^+$.</footnote>
have used the handy notation $[\cdot]^{+} = \max(\cdot, 0)$ since $\rho_{q\ell}^{(\ell)} > 1 \iff \zeta_{q\ell}^{(\ell)} > 0$. Hence, using Proposition 12,

$$ \frac{1}{r_{\ell}} \| \mathbf{X}^{(\ell)\top} \hat{\mathbf{U}}^{(\ell)} \|_{\mathrm{F}}^{2} = \frac{1}{r_{\ell}} \sum_{q_{\ell}=1}^{r_{\ell}} \cos^2 \theta_{q_{\ell}}^{(\ell)} \xrightarrow[N \to +\infty]{\text{a.s.}} \frac{1}{r_{\ell}} \sum_{q_{\ell}=1}^{r_{\ell}} [\zeta_{q_{\ell}}^{(\ell)}]^{+}. \quad (7) $$

Therefore, the quantity $\frac{1}{r_{\ell}} \| \mathbf{X}^{(\ell)\top} \hat{\mathbf{U}}^{(\ell)} \|_{\mathrm{F}}^{2} \in [0, 1]$ appears as a relevant measure of alignment between the singular subspaces of $\mathcal{P}$ and $\hat{\mathcal{T}}$ and does not depend on the chosen orthonormal bases $\mathbf{X}^{(\ell)}$ and $\hat{\mathbf{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{\mathcal{T}}$ and $\mathcal{P} = \sqrt{\omega}\mathcal{P}_{\circ}$ with $\|\mathcal{P}_{\circ}\|^{2} = \sigma_{N}$ as a function of the signal-to-noise ratio $\omega$. The fact that $\|\mathcal{P}_{\circ}\|^{2} = \sigma_{N}$ ensures that the estimation problem is non-trivial (neither too easy nor too hard) as $\rho_{q\ell}^{(\ell)} = \Theta(1)$. Plain curves are the alignments given by Theorem 9 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, $\mathcal{P}_{\circ}$ has multilinear rank (3, 4, 5). Hence, its “energy” $\|\mathcal{P}_{\circ}\|^{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_{q\ell}^{(\ell)} = 1$, that is, $\omega = \sigma_{N}/s_{\ell}^{2}(\mathbf{P}_{\circ}^{(\ell)})$. In other words, there are $r_{\ell}$ phase transitions along mode $\ell$ whose positions depend on the singular values of $\mathbf{P}_{\circ}^{(\ell)}$. If $\omega$ is too small (here, $\omega \lesssim 2$), truncated MLSVD fails to recover any direction of the singular subspaces of $\mathcal{P}_{\circ}$. As $\omega$ passes the first phase transition (here, at $\omega = \sigma_{N}/s_{\ell}^{2}(\mathbf{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 $\mathcal{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 $\mathcal{T}$ with the HOOI algorithm, which is discussed in the next section.

**Remark 13 (Reconstruction of $\mathcal{H}$ from $\mathcal{T}$)** Guarantees on the recovery of the core tensor $\mathcal{H}$ can be deduced from Theorem 9 as well. Without loss of generality, we can assume that $\mathbf{X}^{(\ell)\top}\hat{\mathbf{U}}^{(\ell)}$ is, up to an almost-surely vanishing additive term, a diagonal matrix with entries $\zeta_{1}^{(\ell)}, \dots, \zeta_{r_{\ell}}^{(\ell)}$ (otherwise replace $\hat{\mathbf{U}}^{(\ell)}$ by $\hat{\mathbf{U}}^{(\ell)}\mathbf{O}^{(\ell)}$ for a well-chosen orthogonal matrix $\mathbf{O}^{(\ell)}$). Then, $\hat{\mathcal{G}} = \mathcal{T}(\hat{\mathbf{U}}^{(1)}, \dots, \hat{\mathbf{U}}^{(d)})$ is the corresponding estimator of $\mathcal{H}$ and

$$ \hat{\mathcal{G}} = \hat{\mathcal{H}} + \frac{1}{\sqrt{N}}\mathcal{N}(\hat{\mathbf{U}}^{(1)}, \dots, \hat{\mathbf{U}}^{(d)}) \quad \text{with} \quad \hat{\mathcal{H}} = \left[ \mathcal{H}; \hat{\mathbf{U}}^{(1)\top}\mathbf{X}^{(1)}, \dots, \hat{\mathbf{U}}^{(d)\top}\mathbf{X}^{(d)} \right]. $$

We will see in Lemma 14 below that $\|\frac{1}{\sqrt{N}}\mathcal{N}(\hat{\mathbf{U}}^{(1)}, \dots, \hat{\mathbf{U}}^{(d)})\|_{\mathrm{F}} = \mathcal{O}(1)$ almost surely as $N \to +\infty$. On the other hand, we know that $\|\mathcal{H}\|_{\mathrm{F}} = \|\mathcal{P}\|_{\mathrm{F}} = \Theta(N^{\frac{d-2}{4}}) \gg \mathcal{O}(1)$ as soon
as $d \ge 3$ and the entries of $\hat{\mathcal{H}}$ are proportional to those of $\mathcal{H}$:

$$ \hat{\mathcal{H}}_{q_1, \dots, q_d} = \mathcal{H}_{q_1, \dots, q_d} \prod_{\ell=1}^{d} \zeta_{q_\ell}^{(\ell)} + \epsilon_{q_1, \dots, q_d}, \quad \ell \in [d], \quad q_\ell \in [r_\ell], $$

up to an almost-surely vanishing additive term $\epsilon_{q_1, \dots, q_d}$ as $N \to +\infty$.

Moreover, regarding the reconstruction of $\mathcal{T}$, we know from De Lathauwer et al. (2000b, Property 10) that

$$ \|\mathcal{T} - \hat{\mathcal{T}}\|_{\mathrm{F}}^{2} \le \sum_{\ell=1}^{d} \sum_{i_{\ell}=r_{\ell}+1}^{n_{\ell}} s_{i_{\ell}}^{2}(\mathbf{T}^{(\ell)}) = \sum_{\ell=1}^{d} \left( \|\mathcal{T}\|_{\mathrm{F}}^{2} - \sum_{q_{\ell}=1}^{r_{\ell}} s_{q_{\ell}}^{2}(\mathbf{T}^{(\ell)}) \right) $$

and the asymptotic behavior of $s_{q_\ell}^2(\mathbf{T}^{(\ell)})$ is given by Theorem 9.

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

In search of an efficient estimator of the planted signal $\mathcal{P}$, one naturally considers the best low-multilinear-rank approximation of $\mathcal{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, 2018). 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

$$ (U_{\star}^{(1)}, \dots, U_{\star}^{(d)}) \in \underset{U^{(\ell)} \in V_{r_{\ell}}(\mathbb{R}^{n_{\ell}}), \ell \in [d]}{\arg \max} \frac{1}{2} \|\mathcal{T}(U^{(1)}, \dots, U^{(d)})\|_{\mathrm{F}}^{2} \quad (8) $$

where $V_{r_\ell}(\mathbb{R}^{n_\ell}) \stackrel{\text{def}}{=} \{U^{(\ell)} \in \mathbb{R}^{n_\ell \times r_\ell} | U^{(\ell)\top}U^{(\ell)} = 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 $\mathcal{T}(U^{(1)}, \dots, U^{(d)})$ is equal to the Frobenius norm of any of its unfoldings,

$$ \|\mathcal{T}(U^{(1)}, \dots, U^{(d)})\|_{\mathrm{F}} = \left\| U^{(\ell)\top} \mathbf{T}^{(\ell)} \boxtimes_{\ell' \neq \ell} U^{(\ell')} \right\|_{\mathrm{F}}, \quad \ell \in [d]. $$

And we see from Problem (8) that $U_{\star}^{(\ell)}$ is the matrix gathering the $r_\ell$ dominant left singular vectors of $\mathbf{T}^{(\ell)} \boxtimes_{\ell' \neq \ell} U_{\star}^{(\ell')}$. This is precisely what motivates the HOOI algorithm presented
<caption>**Algorithm 1:** Higher-Order Orthogonal Iteration (De Lathauwer et al., 2000a)</caption>

for $\ell = 1, \dots, d$ do $U_0^{(\ell)} \leftarrow r_\ell$ dominant left singular vectors of $T^{(\ell)}$
repeat
| for $\ell = 1, \dots, d$ do $U_{t+1}^{(\ell)} \leftarrow r_\ell$ dominant left singular vectors of $T^{(\ell)} \boxtimes_{\ell' \neq \ell} U_t^{(\ell')}$
until convergence at $t = T$
$\mathcal{G}_{\text{HOOI}} \leftarrow \mathcal{T}(U_T^{(1)}, \dots, U_T^{(d)})$

in Algorithm 1. It performs fixed-point iterations to compute a solution $U_{\star}^{(1)}, \dots, U_{\star}^{(d)}$ satisfying the previous property$^{17}$.

The HOOI algorithm is initialized with $U_0^{(1)}, \dots, U_0^{(d)}$, the truncated MLSVD$^{18}$ of $\mathcal{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 $U_{\star}^{(1)}, \dots, U_{\star}^{(d)}$ such that $U_{\star}^{(\ell)}$ spans the left $r_{\ell}$-dimensional dominant singular subspace of $T^{(\ell)} \boxtimes_{\ell' \neq \ell} U_{\star}^{(\ell')}$ 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, $\|\mathcal{T}(U_t^{(1)}, \dots, U_t^{(d)})\|_{\text{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 15 below, we show that an initialization with non-zero alignment with the signal $\mathcal{P}$ is sufficient to ensure that the HOOI algorithm perfectly reconstructs it *after a single iteration*.

Before introducing Theorem 15, we formulate an important preliminary result which essentially states that, given $A^{(\ell)} \in V_{r_\ell}(\mathbb{R}^{n_\ell})$, $\ell \in [d]$, the quantity $\frac{1}{\sqrt{N}}\|\mathcal{N}(A^{(1)}, \dots, A^{(d)})\|_{\text{F}}$ is almost surely bounded as $N \to +\infty$.

**Lemma 14** *With probability at least $1 - \delta$, *

$$ \sup_{A^{(\ell)} \in V_{r_\ell}(\mathbb{R}^{n_\ell}), \ell \in [d]} \left\| \mathcal{N}(A^{(1)}, \dots, A^{(d)}) \right\|_{\text{F}}^2 \\ \le 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 D.

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

<footnote n="17">17. In fact, this property corresponds to first-order optimality conditions of Problem (8), with the squared singular values as Lagrange multipliers.</footnote>

<footnote n="18">18. Consistently with the notations of Section 3, this is $\hat{U}^{(1)}, \dots, \hat{U}^{(d)}$.</footnote>
**Theorem 15 (Asymptotic optimality of HOOI)** As $N \to +\infty$, if $\|\mathcal{P}\|_F \gg 1$ and

$$ \min_{\ell \in [d], q_\ell \in [r_\ell]} \|\mathcal{P}(U_0^{(1)}, \dots, x_{q_\ell}^{(\ell)}, \dots, U_0^{(d)})\|_F \stackrel{\text{def}}{=} L_N \gg \|\mathcal{P}\|_F^{1/2} $$

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

$$ \frac{1}{r_\ell} \|X^{(\ell)\top} U_1^{(\ell)}\|_F^2 = 1 + \mathcal{O}\left(\frac{\|\mathcal{P}\|_F}{L_N^2}\right) \quad \text{almost surely.} $$

**Proof** See Appendix E.

<figure><img src="image_4.png" /></figure>

It is important to carefully understand the assumptions of Theorem 15. Firstly, it assumes that $\|\mathcal{P}\|_F \gg 1$ as $N \to +\infty$, that is, the signal is not necessarily in the non-trivial regime $\Theta(N^{d-2}/4)$ but can be smaller or bigger as long as $\|\mathcal{P}\|_F \to +\infty$, regardless its speed. Then, the second assumption $\|\mathcal{P}(U_0^{(1)}, \dots, x_{q_\ell}^{(\ell)}, \dots, U_0^{(d)})\|_F \gg \|\mathcal{P}\|_F^{1/2}$ means that each principal directions of the $\ell$-th singular subspace are sufficiently preserved after contraction on $\{U_0^{(\ell')}\}_{\ell' \neq \ell}$. When these assumptions are verified, Theorem 15 states that the matrices $U_1^{(1)}, \dots, U_1^{(d)}$ computed after the first iteration of HOOI perfectly reconstruct the singular subspaces of the sought signal $\mathcal{P}$ as the dimensions of the tensor, $n_1, \dots, n_d$, grow large. More formally, as $N \to +\infty$, the alignment $\frac{1}{r_\ell}\|X^{(\ell)\top} U_1^{(\ell)}\|_F^2$ approaches 1 almost surely. Furthermore, the speed of this convergence behaves like $\|\mathcal{P}\|_F/L_N^2$.

Theorem 15 does not assume a particular choice of initialization $U_0^{(1)}, \dots, U_0^{(d)}$ 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 = \Theta(\|\mathcal{P}\|)$ is verified as soon as enough principal directions are recovered. According to Theorem 9, this is only possible if $\|\mathcal{P}\|_F \gg \sqrt{\sigma_N} = \Theta(N^{d-2}/4)$ since a necessary condition is $\rho_1^{(\ell)} > 1$ for all $\ell \in [d]$, while a sufficient condition is $\rho_{q_\ell}^{(\ell)} > 1$ for all $\ell \in [d]$ and $q_\ell \in [d]$. In other words, convergence at speed $\|\mathcal{P}\|_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 $\mathcal{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 15 can already be verified as soon as $\min_\ell \rho_1^{(\ell)} > 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_2^{(\ell)} < 1$. If the singular subspaces of $\mathcal{P}$ are sufficiently preserved with the initialization $U_0^{(1)}, \dots, U_0^{(d)}$ — i.e., if $L_N \gg \|\mathcal{P}\|_F^{1/2}$ —, then the other principal directions still emerge after the first iteration.

In the simpler rank-1 case, these technical considerations vanish and we recover the result of Feldman and Donoho (2023, Theorem 4.2): if $\mathcal{P} = \beta_N \otimes_{\ell \in [d]} x^{(\ell)}$ then a necessary and sufficient condition for $L_N = \Theta(\|\mathcal{P}\|_F)$ is simply $\beta_N^2 > \sigma_N$. Indeed, $\rho_1^{(\ell)} = \beta_N^2/\sigma_N$ for all $\ell \in [d]$. Moreover, when this assumption is verified, Theorem 15 ensures the asymptotic exact reconstruction of $x^{(1)}, \dots, x^{(d)}$ in a single iteration with a $\|\mathcal{P}\|_F^{-1} = \beta_N^{-1}$ speed of convergence.
<figure><img src="image_5.png" /><figcaption>Figure 4: Alignments between singular subspaces of the observation $\mathcal{T} = \mathcal{P} + \frac{1}{\sqrt{N}}\mathcal{N}$ and of the signal $\mathcal{P}$, with $\|\mathcal{P}\|_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}\|\mathbf{X}^{(\ell)\top}\mathbf{U}_0^{(\ell)}\|_F^2$. Middle: $\frac{1}{r_\ell}\|\mathbf{X}^{(\ell)\top}\mathbf{U}_1^{(\ell)}\|_F^2$. Right: $(1 - \frac{1}{r_\ell}\|\mathbf{X}^{(\ell)\top}\mathbf{U}_1^{(\ell)}\|_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)$.</figcaption></figure>

**Remark 16 (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 15 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$.

Theorem 15 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 $\|\mathcal{P}\|_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 9. 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 15. This is specified in the right panel where the value $(1 - \frac{1}{r_\ell}\|\mathbf{X}^{(\ell)\top}\mathbf{U}_1^{(\ell)}\|_F^2) \times \sqrt{\sigma_N}$ is plotted. According to Theorem 15, this value should be $O(1)$ since $L_N = \Theta(\|\mathcal{P}\|_F) = \Theta(\sqrt{\sigma_N})$ here. The observed behavior confirms the $\|\mathcal{P}\|_F^{-1}$ speed of convergence asserted in Theorem 15.

## 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 $\|\mathcal{P}\|_F$ is at least $\Theta(N^{\frac{d-2}{4}})$ and its reconstruction performance has been discussed in Section 3.2. However, Theorem 15 suggests that it is possible to perfectly reconstruct the signal $\mathcal{P}$ from the observation $\mathcal{T}$ as long as $\|\mathcal{P}\|_F \gg 1$ and HOOI is accurately initialized. Yet, it is known that, without prior information on $\mathcal{P}$, maximum likelihood estimation (1) is NP-hard below the $\Theta(N^{\frac{d-2}{4}})$ computational threshold (Zhang and Xia, 2018), which lies precisely in the non-trivial regime of truncated MLSVD.

In fact, what can be understood from Theorem 15 is that it suffices to have an initialization $U_0^{(\ell)}$, $\ell \in [d]$, slightly aligned with the underlying signal $\mathcal{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 $\mathcal{P}$ as soon as $\|\mathcal{P}\|_F = \Theta(1)$ (Ben Arous et al., 2019; Jagannath et al., 2020; Zhang and Xia, 2018). Hence, as the HOOI algorithm is meant to compute a maximum likelihood estimator, with the assumption $\|\mathcal{P}\|_F \gg 1$ made in Theorem 15, it is expected that these iterations allow to perfectly recover the signal asymptotically. Maximum likelihood estimation is indeed (theoretically) trivial if $\|\mathcal{P}\|_F \to +\infty$. It is more surprising however that this already happens at the first iteration.

As said previously, the choice of initialization $U_0^{(\ell)}$, $\ell \in [d]$, does not matter in Theorem 15. In fact, without prior information, a truncated MLSVD is the best choice as it allows to partially reconstruct the signal at the $\Theta(N^{\frac{d-2}{4}})$ computational threshold. Nevertheless, had one prior information allowing such a reconstruction in the regime $1 \ll \mathcal{P} \ll N^{\frac{d-2}{4}}$ — where truncated MLSVD would not be fruitful —, HOOI would still be able to perfectly reconstruct the signal $\mathcal{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 $U_0^{(\ell)}$, $\ell \in [d]$, in the right basin of attraction, which, without prior information, is only possible above the $\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 $\Theta(1)$ or $\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 $\mathcal{T} = \sqrt{\omega}\mathcal{P}_o + \frac{1}{\sqrt{N}}\mathcal{N}$. Although $\|\mathcal{P}_o\|_F = \sqrt{\sigma_N} = \Theta(N^{\frac{d-2}{4}})$, HOOI does not achieve perfect recovery of $\mathcal{P}$ as one might expect from Theorem 15 (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 yields theoretical and practical insights into the estimation of a low-rank signal from an observation $\mathcal{T} = \mathcal{P} + \frac{1}{\sqrt{N}}\mathcal{N}$ following the most general spiked tensor model. While Zhang and Xia (2018) gave a general overview of the different regimes governing the estimation of $\mathcal{P}$ with a low-multilinear-rank approximation of $\mathcal{T}$ — thereby confirming the existence of a computational-to-statistical gap —, our results
shed light on the non-trivial aspects at stake around the $\Theta(N^{\frac{d-2}{4}})$ computational thresh-
old. 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 iter-
ation, which are efficient above a $\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 $\mathcal{T}$ in the large $N$ limit. Specifi-
cally, our first main result shows that, when properly centered and scaled, the eigenvalues of
$\mathcal{T}^{(\ell)}\mathcal{T}^{(\ell)\top}$ asymptotically follow a semicircle distribution. The rescaling exhibits their mean
$\mu_N^{(\ell)} = \Theta(N^{d-2})$ and a quantity $\sigma_N = \Theta(N^{\frac{d-2}{2}})$ governing their spread. From our denois-
ing 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 $\mathcal{P}$. Yet,
the existence of these outlier eigenvalues hinges on the values of the signal-to-noise ratios
$\rho_{q_\ell}^{(\ell)} = s_{q_\ell}^2(\mathcal{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_{q_\ell}^{(\ell)}$. These results justify the use of a truncated MLSVD
to estimate $\mathcal{P}$ from the observation $\mathcal{T}$ and allow the precise characterization of the achiev-
able reconstruction performances in the non-trivial regime, i.e., close to the computational
threshold. In particular, we have seen that each singular value of $\mathcal{P}^{(\ell)}$ determines the posi-
tion 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 $\mathcal{P}$ in a sense precised in Theorem 15, 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_{q_\ell}^{(\ell)}(\mathcal{P}^{(\ell)})$ of the unfoldings of $\mathcal{P}$. Yet, given
prior information, HOOI can still reconstruct the maximum likelihood solution below the
computational threshold $\|\mathcal{P}\|_F = \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 approxima-
tion problem near the computational threshold, which has both practical and theoretical
implications. Besides, from a theoretical perspective, the behavior of the maximum likeli-
hood estimator is still unclear near the statistical threshold — that is, in the regime where
$\|\mathcal{P}\|_F = \Theta(1)$. Several works have studied the rank-1 symmetric case (Ben Arous et al.,
2019; Jagannath et al., 2020) 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 to low-rank tensor estimation from spiked models.

## Acknowledgments and Disclosure of Funding

This work is supported by the MIAI LargeDATA Chair at Université Grenoble Alpes. The authors acknowledge support of the Institut Henri Poincaré (UAR 839 CNRS-Sorbonne Université), and LabEx CARMIN (ANR-10-LABX-59-01)

## Appendix A. More Random Matrix Tools

In this appendix we state a few useful results for our proofs in Appendices B and C. This completes the basic tools given in Section 2.3. Let us also recall valuable references for the interested reader: Bai and Silverstein (2010); Pastur and Shcherbina (2011); Tao (2012); Potters and Bouchaud (2020); Couillet and Liao (2022).

### A.1 Stieltjes Transforms

Following the definition of Stieltjes transform (Definition 2), we state a few of its properties. They can be found in Pastur and Shcherbina (2011, Proposition 2.1.2).

**Proposition 17 (Properties of Stieltjes transforms)** Let $m$ be the Stieltjes transform of a real probability measure $\mu$.

1. $m$ is analytic on $\mathbb{C} \setminus \mathbb{R}$ and $m(\bar{z}) = \overline{m(z)}$;

2. if $\Im z \neq 0$ then $\Im[m(z)]\Im z > 0$;

3. for all $z \in \mathbb{C} \setminus \mathbb{R}$, $|m(z)| \le 1/|\Im z|$;

4. if $\mu$ has a density at $x \in \mathbb{R}$ then $\frac{d\mu}{dx}(x) = \lim_{y \downarrow 0} \Im[m(x + iy)]$.

Note that the resolvent $Q_S(z) \stackrel{\text{def}}{=} (S - zI_n)^{-1}$ of a symmetric matrix $S$ also satisfies $\|Q_S(z)\| \le 1/|\Im z|$.

In order to take into account the potential divergence of $m(z)$ or $\|Q_S(z)\|$ near the real axis, we introduce the following notation.

**Definition 18** We write $u_n(z) = \mathcal{O}_z(v_n)$ whenever there exist two polynomials $P, Q$ with nonnegative coefficients and an integer $n_0$ such that $|u_n(z)| \le \frac{P(|z|)}{|\Im z|Q(|\Im z|)}|v_n|$ as soon as $n \ge n_0$.

In particular, $m(z) = \mathcal{O}_z(1)$ and $\|Q_S(z)\| = \mathcal{O}_z(1)$.
## A.2 Concentration Tools

In order to prove the concentration around its expectation of a scalar quantity, we will rely on the two following lemmas.

**Lemma 19 (Poincaré-Nash inequality)** Let $z \sim \mathcal{N}(0_p, I_p)$ and $f: \mathbb{R}^p \to \mathbb{C}$ be a differentiable function with polynomially bounded partial derivatives $\partial_1 f, \dots, \partial_p f$. Then,

$$ \mathrm{Var}(f(z)) \leq \mathbb{E}[\|\nabla f(z)\|^2] $$

where $\nabla = [\partial_1 \ \dots \ \partial_p]^T$.

**Proof** See Pastur and Shcherbina (2011, Proposition 2.1.6). ■

**Lemma 20** Let $(X_n)_{n \geq 0}$ be a sequence of random variables on $\mathbb{C}$. If $\mathbb{E}[X_n] \to \bar{X} \in \mathbb{C}$ as $n \to +\infty$ and there exists an integer $k \geq 1$ such that $\sum_{n \geq 0} \mathbb{E}[|X_n - \mathbb{E}[X_n]|^k] < +\infty$, then $X_n \to \bar{X}$ almost surely as $n \to +\infty$.

**Proof** Let $\varepsilon > 0$. By Markov's inequality (Billingsley, 2012, equation 5.31),

$$ \sum_{n \geq 0} \mathbb{P}(|X_n - \mathbb{E}[X_n]| \geq \varepsilon) \leq \frac{1}{\varepsilon^{\kappa}} \sum_{n \geq 0} \mathbb{E}[|X_n - \mathbb{E}[X_n]|^{\kappa}] < +\infty. $$

Thus, by the first Borel-Cantelli lemma<fnref n="19" /> (Billingsley, 2012, Theorem 4.3),

$$ \mathbb{P}\left(\limsup_{n \to +\infty} |X_n - \mathbb{E}[X_n]| \geq \varepsilon\right) = 0. $$

This implies the almost sure convergence of $(X_n)_{n \geq 0}$ to $\bar{X}$:

$$ \mathbb{P}\left(\limsup_{n \to +\infty} |X_n - \bar{X}| \geq \varepsilon\right) \leq \mathbb{P}\left(\limsup_{n \to +\infty} \left[|X_n - \mathbb{E}[X_n]| + |\mathbb{E}[X_n] - \bar{X}|\right] \geq \varepsilon\right) = 0. $$

## A.3 The Hellffer-Sjöstrand Formula

A basic property of Stieltjes transforms is that, if $f: \mathbb{R} \to \mathbb{R}$ is a bounded continuous function and $m$ is the Stieltjes transform of a real probability measure, then $\int f d\mu = \lim_{y \downarrow 0} \frac{1}{\pi} \int f(x)m(x + iy)dx$. Yet, this formula is not always satisfying because $m(x + iy)$ may diverge as $y \downarrow 0$. The following powerful result gives a more pleasing formula when $f$ is smooth and compactly supported.

<footnote n="19">19. Given a sequence of events $(A_n)_{n \geq 0}$ in a probability space $(\Omega, \mathcal{A}, \mathbb{P})$, if $\sum_{n \geq 0} \mathbb{P}(A_n)$ converges then $\mathbb{P}(\limsup_{n \geq 0} A_n) = 0$.</footnote>
**Proposition 21 (Helffer-Sjöstrand formula)** Let $\mu$ be a probability measure on $\mathbb{R}$ and $m$ be its Stieltjes transform. Let $f: \mathbb{R} \to \mathbb{R}$ be a compactly supported function which has a continuous $(k+1)$-th derivative $(k \ge 1)$. We define $\Phi_k[f]$, the quasi-analytic extension of $f$ on $\mathbb{C}^+ \stackrel{\text{def}}{=} \{z \in \mathbb{C} \mid \Im z > 0\}$ as

$$ \Phi_k[f](z) = \sum_{l=0}^{k} \frac{(i\Im z)^l}{l!} f^{(l)}(\Re z) \chi(\Im z) $$

where $\chi: \mathbb{R} \to [0,1]$ is an infinitely differentiable even function such that$^{20}$ $\chi(y) = 0$ if $|y| \ge 1$ and $\chi(y) = 1$ if $|y| \le \delta$ for some $\delta \in [0,1]$. Then,

$$ \int_{\mathbb{R}} f \, d\mu = \frac{2}{\pi} \Re \int_{\mathbb{C}^{+}} \frac{\partial \Phi_k[f]}{\partial \bar{z}}(z) m(z) \, dz $$

where $\frac{\partial}{\partial \bar{z}} = \frac{1}{2} \left( \frac{\partial}{\partial \Re z} + i \frac{\partial}{\partial \Im z} \right)$ is the Wirtinger derivative.

**Proof** The support of $\Phi_k[f]$ is compact therefore an integration by parts gives

$$ \begin{aligned} & \frac{2}{\pi} \Re \int_{\mathbb{C}^+} \frac{\partial \Phi_k[f]}{\partial \bar{z}}(z) m(z) \, dz \\ &= \frac{2}{\pi} \Re \left[ \frac{1}{2} \int_0^{+\infty} \left( \int_{-\infty}^{+\infty} \frac{\partial \Phi_k[f]}{\partial x}(x + iy) m(x + iy) \, dx \right) dy \right. \\ & \qquad \left. + \frac{i}{2} \int_{-\infty}^{+\infty} \left( \int_0^{+\infty} \frac{\partial \Phi_k[f]}{\partial y}(x + iy) m(x + iy) \, dy \right) dx \right] \\ &= \frac{2}{\pi} \Re \left[ \frac{-i}{2} \lim_{y \downarrow 0} \{ \Phi_k[f](x + iy) m(x + iy) \} \, dx - \int_{\mathbb{C}^+} \Phi_k[f](z) \frac{\partial m}{\partial \bar{z}}(z) \, dz \right]. \end{aligned} $$

Since $m$ is an analytic function, the Cauchy-Riemann equations give $\partial m / \partial \bar{z} = 0$ and the second integral of the right-hand side is zero. Then, as the limit of a product is the product of the limits and $\lim_{y \downarrow 0} \Phi[f](x + iy) = f(x)$, we have

$$ \frac{2}{\pi} \Re \int_{\mathbb{C}^+} \frac{\partial \Phi_k[f]}{\partial \bar{z}}(z) m(z) \, dz = \Re \left[ \frac{-i}{\pi} \int_{\mathbb{R}} f(x) \left( \lim_{y \downarrow 0} m(x + iy) \right) \, dx \right]. $$

Finally, observing that $\Re[-iz] = \Im z$ and $\frac{dx}{\pi} \lim_{y \downarrow 0} \Im[m(x + iy)] = d\mu(x)$ (Proposition 17) concludes the proof. ■

Since $\chi$ is constant on $[0, \delta]$, $0 < \Im z \le \delta \implies \frac{\partial \Phi_k[f]}{\partial \bar{z}}(z) = \frac{1}{2} \frac{(i \Im z)^k}{k!} f^{(k+1)}(\Re z)$. The greater $k$, the faster the convergence to zero of $\frac{\partial \Phi_k[f]}{\partial \bar{z}}$ as $\Im z \downarrow 0$. This is particularly useful to evaluate quantities such as $\int f d\mu_1 - \int f d\mu_2$ from the knowledge that the Stieltjes transforms $m_1, m_2$ of $\mu_1, \mu_2$ (respectively) satisfy $|m_1(z) - m_2(z)| = \mathcal{O}(|\Im z|^{-k})$ as $\Im z \downarrow 0$:

$$ \int f \, d\mu_1 - \int f \, d\mu_2 = \frac{2}{\pi} \Re \int_{\mathbb{C}^+} \underbrace{\frac{\partial \Phi_k[f]}{\partial \bar{z}}(z)(m_1(z) - m_2(z))}_{=\mathcal{O}(1) \text{ as } \Im z \downarrow 0} \, dz. $$

<footnote n="20">20. An example of such function is $\chi: y \mapsto g\left(\frac{1+y}{1-\delta}\right)g\left(\frac{1-y}{1-\delta}\right)$ with $g(u) = \frac{f(u)}{f(u)+f(1-u)}$ and $f(u) = \begin{cases} \exp(-1/u) & \text{if } u > 0 \\ 0 & \text{if } u \le 0. \end{cases}$</footnote>
The smoothness of $f$ compensates for the divergence of $m_1(z) - m_2(z)$ near the real axis.

The next lemma is a handy result to evaluate some integrals which appear after using the Helffer-Sjöstrand formula.

**Lemma 22** Let $K$ be a compact subset of $\mathbb{R}$ and $h$ be an analytic function on $\mathbb{C} \setminus K$ such that

1. $\lim_{|z| \to +\infty} h(z) = 0$,

2. $h(\bar{z}) = \overline{h(z)}$ for all $z \in \mathbb{C} \setminus K$,

3. there exist an integer $n_0 \ge 1$ and a constant $C_K > 0$ such that

$$|h(z)| \le C_K \max(\mathrm{Dist}(z, K)^{-n_0}, 1) \quad \text{for all } z \in \mathbb{C} \setminus K.$$

$$\text{Then, } \lim_{y \downarrow 0} \frac{1}{\pi} \int_{\mathbb{R}} \mathfrak{S}[h(x + iy)] \, dx = \lim_{y \to +\infty} -iyh(iy).$$

**Proof** According to Theorem 5.4 of Schultz (2005) (see also Capitaine et al., 2009, Theorem 4.3; Loubaton, 2016, Lemma 9.1), $h$ is the Stieltjes transform of a compactly-supported Schwartz distribution (Rudin, 1991, Chapter 6) and therefore satisfies the stated result. ■

## Appendix B. Proof of Theorem 5

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

$$Q^{(\ell)}(z) \stackrel{\text{def}}{=} (T^{(\ell)}T^{(\ell)\top} - zI_{n_\ell})^{-1}.$$

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

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

$$P^{(\ell)}T^{(\ell)\top}Q^{(\ell)} + \frac{1}{\sqrt{N}}N^{(\ell)}T^{(\ell)\top}Q^{(\ell)} - zQ^{(\ell)} = I_{n_\ell}. \quad (9)$$

### B.1 Expressions with Stein's Lemma

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

$$\mathbb{E}[\boldsymbol{P}^{(\ell)}\boldsymbol{T}^{(\ell)\top}\boldsymbol{Q}^{(\ell)}] = \mathbb{E}[\boldsymbol{P}^{(\ell)}\boldsymbol{P}^{(\ell)\top}\boldsymbol{Q}^{(\ell)}] - \frac{1}{N}\mathbb{E}[\boldsymbol{P}^{(\ell)}\boldsymbol{T}^{(\ell)\top}\boldsymbol{Q}^{(\ell)}\mathrm{Tr}\boldsymbol{Q}^{(\ell)} + \boldsymbol{P}^{(\ell)}\boldsymbol{T}^{(\ell)\top}\boldsymbol{Q}^{(\ell)2}], \quad (10)$$

$$\mathbb{E}[\boldsymbol{N}^{(\ell)}\boldsymbol{T}^{(\ell)\top}\boldsymbol{Q}^{(\ell)}] = \frac{\prod_{\ell' \ne \ell} n_{\ell'}}{\sqrt{N}} \mathbb{E}[\boldsymbol{Q}^{(\ell)}] - \frac{1}{\sqrt{N}} \mathbb{E}[(n_{\ell} + 1)\boldsymbol{Q}^{(\ell)} + z(\boldsymbol{Q}^{(\ell)2} + \boldsymbol{Q}^{(\ell)}\mathrm{Tr}\boldsymbol{Q}^{(\ell)})]. \quad (11)$$
**Derivatives of $Q^{(\ell)}$** Firstly, we need to show that

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

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

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

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

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

**Proof of Equation (10)** Since $\mathbf{T}^{(\ell)} = \mathbf{P}^{(\ell)} + \frac{1}{\sqrt{N}} \mathbf{N}^{(\ell)}$, we have,

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

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

$$ \begin{aligned} \mathbb{E}[\mathbf{P}^{(\ell)} \mathbf{N}^{(\ell)\top} \mathbf{Q}^{(\ell)}]_{i,j} &= \sum_{k=1}^{\prod_{\ell' \neq \ell} n_{\ell'}} \sum_{l=1}^{n_{\ell}} \mathbb{E}[P_{i,k}^{(\ell)} N_{l,k}^{(\ell)} Q_{l,j}^{(\ell)}] \\ &= \sum_{k=1}^{\prod_{\ell' \neq \ell} n_{\ell'}} \sum_{l=1}^{n_{\ell}} \mathbb{E}\left[P_{i,k}^{(\ell)} \frac{\partial Q_{l,j}^{(\ell)}}{\partial N_{l,k}^{(\ell)}}\right] \\ &= -\frac{1}{\sqrt{N}} \sum_{k=1}^{\prod_{\ell' \neq \ell} n_{\ell'}} \sum_{l=1}^{n_{\ell}} \mathbb{E}\left[P_{i,k}^{(\ell)} Q_{l,l}^{(\ell)} [\mathbf{T}^{(\ell)\top} \mathbf{Q}^{(\ell)}]_{k,j}\right] \\ &\quad - \frac{1}{\sqrt{N}} \sum_{k=1}^{\prod_{\ell' \neq \ell} n_{\ell'}} \sum_{l=1}^{n_{\ell}} \mathbb{E}\left[P_{i,k}^{(\ell)} [\mathbf{Q}^{(\ell)} \mathbf{T}^{(\ell)}]_{l,k} Q_{l,j}^{(\ell)}\right] \\ &= -\frac{1}{\sqrt{N}} \mathbb{E}\left[\mathbf{P}^{(\ell)} \mathbf{T}^{(\ell)\top} \mathbf{Q}^{(\ell)} \operatorname{Tr} \mathbf{Q}^{(\ell)} + \mathbf{P}^{(\ell)} \mathbf{T}^{(\ell)\top} \mathbf{Q}^{(\ell)2}\right]_{i,j}. \end{aligned} $$
**Proof of Equation (11)** We proceed similarly with Stein's lemma (Lemma 4) and Equation (12).

$$
\begin{align*}
& \mathbb{E}[\mathbf{N}^{(\ell)} \mathbf{T}^{(\ell)\top} \mathbf{Q}^{(\ell)}]_{i,j} \\
&= \sum_{k=1}^{\prod_{\ell' \neq \ell} n_{\ell'}} \sum_{l=1}^{n_{\ell}} \mathbb{E}[N_{i,k}^{(\ell)} T_{l,k}^{(\ell)} Q_{l,j}^{(\ell)}] \\
&= \sum_{k=1}^{\prod_{\ell' \neq \ell} n_{\ell'}} \sum_{l=1}^{n_{\ell}} \mathbb{E}\left[\frac{\partial T_{l,k}^{(\ell)}}{\partial N_{i,k}^{(\ell)}} Q_{l,j}^{(\ell)} + T_{l,k}^{(\ell)} \frac{\partial Q_{l,j}^{(\ell)}}{\partial N_{i,k}^{(\ell)}}\right] \\
&= \frac{\prod_{\ell' \neq \ell} n_{\ell'}}{\sqrt{N}} \mathbb{E}[\mathbf{Q}^{(\ell)}]_{i,j} - \frac{1}{\sqrt{N}} \sum_{k=1}^{\prod_{\ell' \neq \ell} n_{\ell'}} \sum_{l=1}^{n_{\ell}} \mathbb{E}\left[T_{l,k}^{(\ell)} \mathbf{Q}_{l,i}^{(\ell)} [\mathbf{T}^{(\ell)\top} \mathbf{Q}^{(\ell)}]_{k,j}\right] \\
&\quad - \frac{1}{\sqrt{N}} \sum_{k=1}^{\prod_{\ell' \neq \ell} n_{\ell'}} \sum_{l=1}^{n_{\ell}} \mathbb{E}\left[T_{l,k}^{(\ell)} [\mathbf{Q}^{(\ell)} \mathbf{T}^{(\ell)}]_{l,k} \mathbf{Q}_{i,j}^{(\ell)}\right] \\
&= \frac{\prod_{\ell' \neq \ell} n_{\ell'}}{\sqrt{N}} \mathbb{E}[\mathbf{Q}^{(\ell)}]_{i,j} - \frac{1}{\sqrt{N}} \mathbb{E}[\mathbf{Q}^{(\ell)} \mathbf{T}^{(\ell)} \mathbf{T}^{(\ell)\top} \mathbf{Q}^{(\ell)} + \mathbf{Q}^{(\ell)} \mathrm{Tr} \mathbf{T}^{(\ell)} \mathbf{T}^{(\ell)\top} \mathbf{Q}^{(\ell)}]_{i,j}.
\end{align*}
$$

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

## B.2 Asymptotic Behavior of the Resolvent

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

$$
\mathbb{E}[\mathbf{P}^{(\ell)} \mathbf{T}^{(\ell)\top} \mathbf{Q}^{(\ell)}] + \frac{\prod_{\ell' \neq \ell} n_{\ell'}}{N} \mathbb{E}[\mathbf{Q}^{(\ell)}] - \frac{1}{N} \mathbb{E}[(n_{\ell} + 1)\mathbf{Q}^{(\ell)} + z(\mathbf{Q}^{(\ell)2} + \mathbf{Q}^{(\ell)} \mathrm{Tr} \mathbf{Q}^{(\ell)})] - z \mathbb{E}[\mathbf{Q}^{(\ell)}] = \mathbf{I}_{n_{\ell}}.
$$

We rearrange this expression into the more convenient following form,

$$
\begin{align*}
& z \frac{n_{\ell}}{N} \mathbb{E} \left[ \mathbf{Q}^{(\ell)} \frac{\mathrm{Tr} \mathbf{Q}^{(\ell)}}{n_{\ell}} \right] + \left( z + \frac{n_{\ell} - \prod_{\ell' \neq \ell} n_{\ell'}}{N} \right) \mathbb{E} \left[ \mathbf{Q}^{(\ell)} \right] + \mathbf{I}_{n_{\ell}} \\
&= - \frac{1}{N} \mathbb{E} \left[ \mathbf{Q}^{(\ell)} + z \mathbf{Q}^{(\ell)2} \right] + \mathbb{E} \left[ \mathbf{P}^{(\ell)} \mathbf{T}^{(\ell)\top} \mathbf{Q}^{(\ell)} \right].
\end{align*}
$$

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

$$
\begin{align*}
\tilde{z} &\stackrel{\text{def}}{=} \frac{z - \mu_N^{(\ell)}}{\sigma_N}, & \tilde{\mathbf{Q}}^{(\ell)}(\tilde{z}) &\stackrel{\text{def}}{=} \left( \frac{\mathbf{T}^{(\ell)}\mathbf{T}^{(\ell)\top} - \mu_N^{(\ell)}\mathbf{I}_{n_{\ell}}}{\sigma_N} - \tilde{z}\mathbf{I}_{n_{\ell}} \right)^{-1} = \sigma_N \mathbf{Q}^{(\ell)}(z), \\
\text{with} \quad & \mu_N^{(\ell)} = \frac{1}{N} \prod_{\ell' \neq \ell} n_{\ell'}, & \sigma_N &= \frac{1}{N} \sqrt{\prod_{\ell \in [d]} n_{\ell}}.
\end{align*}
$$
This changes our equation into

$$ \frac{\mu_N^{(\ell)} + \sigma_N \tilde{z} n_\ell}{\sigma_N^2 N} \mathbb{E} \left[ \tilde{Q}^{(\ell)} \frac{\text{Tr} \tilde{Q}^{(\ell)}}{n_\ell} \right] + \left( \frac{\mu_N^{(\ell)}}{\sigma_N} + \tilde{z} + \frac{n_\ell - \prod_{\ell' \neq \ell} n_{\ell'}}{\sigma_N N} \right) \mathbb{E} \left[ \tilde{Q}^{(\ell)} \right] + I_{n_\ell} \\ = - \frac{1}{\sigma_N N} \mathbb{E} \left[ \tilde{Q}^{(\ell)} + \left( \frac{\mu_N^{(\ell)}}{\sigma_N} + \tilde{z} \right) \tilde{Q}^{(\ell)2} \right] + \frac{1}{\sigma_N} \mathbb{E} \left[ \boldsymbol{P}^{(\ell)} \boldsymbol{T}^{(\ell)\top} \tilde{Q}^{(\ell)} \right]. $$

Let us define $\tilde{m}_N^{(\ell)} : \tilde{z} \mapsto \frac{1}{n_\ell} \text{Tr} \tilde{Q}^{(\ell)}(\tilde{z})$, the Stieltjes transform of the empirical spectral distribution of $\frac{1}{\sigma_N} [\boldsymbol{T}^{(\ell)} \boldsymbol{T}^{(\ell)\top} - \mu_N^{(\ell)} \boldsymbol{I}_{n_\ell}]$. With this definition, and using the fact that $\frac{\mu_N^{(\ell)}}{\sigma_N^2 N} = 1$, we can rewrite the previous relation as

$$ \mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z}) \tilde{Q}^{(\ell)}] + \tilde{z} \mathbb{E}[\tilde{Q}^{(\ell)}] + I_{n_\ell} - \frac{1}{\sigma_N} \mathbb{E}[\boldsymbol{P}^{(\ell)} \boldsymbol{T}^{(\ell)\top} \tilde{Q}^{(\ell)}] \\ = -\frac{1}{\sigma_N} \frac{n_\ell}{N} \mathbb{E}[(\tilde{z} \tilde{m}_N^{(\ell)}(\tilde{z}) + 1) \tilde{Q}^{(\ell)}] - \frac{1}{N} \left( \frac{\mu_N^{(\ell)}}{\sigma_N^2} + \frac{\tilde{z}}{\sigma_N} \right) \mathbb{E}[\tilde{Q}^{(\ell)2}] - \frac{1}{\sigma_N} \frac{1}{N} \mathbb{E}[\tilde{Q}^{(\ell)}] \quad (13) $$

where we have kept only the non-vanishing terms on the left-hand side.

## B.3 Concentration of Bilinear Forms and Traces

### B.3.1 OVERVIEW

Let us show that $\mathbb{E}[\tilde{Q}^{(\ell)}]$ is a deterministic equivalent of $\tilde{Q}^{(\ell)}$ (Definition 3). That is, for all bounded (sequences of) vectors $\boldsymbol{a}, \boldsymbol{b} \in \mathbb{R}^{n_\ell}$ and matrices $\boldsymbol{A} \in \mathbb{R}^{n_\ell \times n_\ell}$,

$$ \boldsymbol{a}^{\top} (\tilde{Q}^{(\ell)} - \mathbb{E}[\tilde{Q}^{(\ell)}]) \boldsymbol{b} \xrightarrow{N \to +\infty} 0 \quad \text{and} \quad \frac{1}{n_\ell} \text{Tr} \boldsymbol{A} (\tilde{Q}^{(\ell)} - \mathbb{E}[\tilde{Q}^{(\ell)}]) \xrightarrow{N \to +\infty} 0. $$

This result is derived with the Poincaré-Nash inequality (Lemma 19) and Lemma 20. Indeed, relying on Lemma 19 and Equation (12), we find below that

$$ \mathbb{E} \left[ \left| \text{Tr} \boldsymbol{A} (\tilde{Q}^{(\ell)}(\tilde{z}) - \mathbb{E}[\tilde{Q}^{(\ell)}(\tilde{z})]) \right|^2 \right] = \mathcal{O}_{\tilde{z}}(N^{-2}), $$

$$ \mathbb{E} \left[ \left| \boldsymbol{a}^{\top} (\tilde{Q}^{(\ell)}(\tilde{z}) - \mathbb{E}[\tilde{Q}^{(\ell)}(\tilde{z})]) \boldsymbol{b} \right|^4 \right] = \mathcal{O}_{\tilde{z}}(N^{-2}) $$

for all deterministic matrix $\boldsymbol{A} \in \mathbb{R}^{n_\ell \times n_\ell}$ and vectors $\boldsymbol{a}, \boldsymbol{b} \in \mathbb{R}^{n_\ell}$ of bounded norms (spectral and Euclidean norms respectively). By Lemma 20, this implies that

$$ \text{Tr} \boldsymbol{A} (\tilde{Q}^{(\ell)}(\tilde{z}) - \mathbb{E}[\tilde{Q}^{(\ell)}(\tilde{z})]) \xrightarrow{N \to +\infty} 0 \quad \text{and} \quad \boldsymbol{a}^{\top} (\tilde{Q}^{(\ell)}(\tilde{z}) - \mathbb{E}[\tilde{Q}^{(\ell)}(\tilde{z})]) \boldsymbol{b} \xrightarrow{N \to +\infty} 0. $$

Hence, $\mathbb{E}[\tilde{Q}^{(\ell)}(\tilde{z})]$ is a deterministic equivalent of $\tilde{Q}^{(\ell)}(\tilde{z})$ according to Definition 3.
### B.3.2 PROOF OF THE CONCENTRATION

**Concentration of Traces** Let $\mathbf{A} \in \mathbb{R}^{n_{\ell} \times n_{\ell}}$ be a deterministic matrix with bounded norm. The Poincaré-Nash inequality (Lemma 19) yields

$$ \mathbb{E}\left[\left(\frac{1}{n_{\ell}} \operatorname{Tr} \mathbf{A}(\tilde{\mathbf{Q}}^{(\ell)} - \mathbb{E}[\tilde{\mathbf{Q}}^{(\ell)}])\right)^2\right] = \operatorname{Var}\left(\frac{1}{n_{\ell}} \operatorname{Tr} \mathbf{A}\tilde{\mathbf{Q}}^{(\ell)}\right) \\ \leq \frac{1}{n_{\ell}^2} \sum_{i=1}^{n_{\ell}} \prod_{\ell' \neq \ell} \sum_{j=1}^{n_{\ell'}} \mathbb{E}\left[\left(\frac{\partial \operatorname{Tr} \mathbf{A}\tilde{\mathbf{Q}}^{(\ell)}}{\partial N_{i,j}^{(\ell)}}\right)^2\right]. $$

With Equation (12), we have,

$$ \begin{aligned} \frac{\partial \operatorname{Tr} \mathbf{A}\tilde{\mathbf{Q}}^{(\ell)}}{\partial N_{i,j}^{(\ell)}} &= \sum_{k=1}^{n_{\ell}} \sum_{l=1}^{n_{\ell}} A_{k,l} \frac{\partial \tilde{Q}_{l,k}^{(\ell)}}{\partial N_{i,j}^{(\ell)}} \\ &= -\frac{1}{\sigma_N \sqrt{N}} \sum_{k=1}^{n_{\ell}} \sum_{l=1}^{n_{\ell}} A_{k,l} [\tilde{Q}_{l,i}^{(\ell)} [\mathbf{T}^{(\ell)\top} \tilde{\mathbf{Q}}^{(\ell)}]_{j,k} + [\tilde{\mathbf{Q}}^{(\ell)} \mathbf{T}^{(\ell)}]_{l,j} \tilde{Q}_{i,k}^{(\ell)}] \\ &= -\frac{1}{\sigma_N \sqrt{N}} [\tilde{\mathbf{Q}}^{(\ell)} \mathbf{A}^{\top} \tilde{\mathbf{Q}}^{(\ell)} \mathbf{T}^{(\ell)} + \tilde{\mathbf{Q}}^{(\ell)} \mathbf{A} \tilde{\mathbf{Q}}^{(\ell)} \mathbf{T}^{(\ell)}]_{i,j}. \end{aligned} $$

Thus, using the fact that $(a+b)^2 \leq 2(a^2+b^2)$, we have,

$$ \begin{aligned} & \frac{1}{n_{\ell}^2} \sum_{i=1}^{n_{\ell}} \prod_{\ell' \neq \ell} \sum_{j=1}^{n_{\ell'}} \mathbb{E}\left[\left(\frac{\partial \operatorname{Tr} \mathbf{A}\tilde{\mathbf{Q}}^{(\ell)}}{\partial N_{i,j}^{(\ell)}}\right)^2\right] \\ & \leq \frac{2}{\sigma_N^2 N n_{\ell}^2} \sum_{i=1}^{n_{\ell}} \prod_{\ell' \neq \ell} \sum_{j=1}^{n_{\ell'}} \mathbb{E}\left[\left[\tilde{\mathbf{Q}}^{(\ell)} \mathbf{A}^{\top} \tilde{\mathbf{Q}}^{(\ell)} \mathbf{T}^{(\ell)}\right]_{i,j}^2 + \left[\tilde{\mathbf{Q}}^{(\ell)} \mathbf{A} \tilde{\mathbf{Q}}^{(\ell)} \mathbf{T}^{(\ell)}\right]_{i,j}^2\right] \\ & = \frac{2}{\sigma_N^2 N n_{\ell}^2} \mathbb{E}\left[\left\|\tilde{\mathbf{Q}}^{(\ell)} \mathbf{A}^{\top} \tilde{\mathbf{Q}}^{(\ell)} \mathbf{T}^{(\ell)}\right\|_{\mathrm{F}}^2 + \left\|\tilde{\mathbf{Q}}^{(\ell)} \mathbf{A} \tilde{\mathbf{Q}}^{(\ell)} \mathbf{T}^{(\ell)}\right\|_{\mathrm{F}}^2\right] \\ & \leq \frac{4}{\sigma_N^2 N n_{\ell}^2} \mathbb{E}\left[\left\|\tilde{\mathbf{Q}}^{(\ell)}\right\|^4 \left\|\mathbf{T}^{(\ell)}\right\|^2\right] \end{aligned} $$

where the last inequality is obtained using the property $\|\mathbf{A}\mathbf{B}\|_{\mathrm{F}} \leq \|\mathbf{A}\|\|\mathbf{B}\|_{\mathrm{F}}$ and the sub-multiplicative property $\|\mathbf{A}\mathbf{B}\| \leq \|\mathbf{A}\|\|\mathbf{B}\|$. Since $\|\mathbf{A}\|_{\mathrm{F}}^2 \leq n_{\ell}\|\mathbf{A}\|^2 = \mathcal{O}(N)$ (because $n_{\ell} = \Theta(N)$ and $\|\mathbf{A}\|$ is bounded), $\|\tilde{\mathbf{Q}}^{(\ell)}\| = \mathcal{O}_{\tilde{z}}(1)$ and $\frac{1}{\sigma_N^2}\|\mathbf{T}^{(\ell)}\|^2 = \mathcal{O}(1)$, we find that $\operatorname{Var}(\frac{1}{n_{\ell}} \operatorname{Tr} \mathbf{A}\tilde{\mathbf{Q}}^{(\ell)}) = \mathcal{O}_{\tilde{z}}(N^{-2})$; Q.E.D.

**Concentration of Bilinear Forms** Let $\mathbf{a}, \mathbf{b} \in \mathbb{R}^{n_{\ell}}$ be two deterministic vectors with bounded norms. The Poincaré-Nash inequality (Lemma 19) yields

$$ \mathbb{E}\left[\left(\mathbf{a}^{\top}(\tilde{\mathbf{Q}}^{(\ell)} - \mathbb{E}[\tilde{\mathbf{Q}}^{(\ell)}])\mathbf{b}\right)^2\right] = \operatorname{Var}(\mathbf{a}^{\top}\tilde{\mathbf{Q}}^{(\ell)}\mathbf{b}) \leq \sum_{i=1}^{n_{\ell}} \prod_{\ell' \neq \ell} \sum_{j=1}^{n_{\ell'}} \mathbb{E}\left[\left(\frac{\partial \mathbf{a}^{\top}\tilde{\mathbf{Q}}^{(\ell)}\mathbf{b}}{\partial N_{i,j}^{(\ell)}}\right)^2\right]. $$
With Equation (12), we have,

$$
\begin{align*}
\frac{\partial \mathbf{a}^{\top} \tilde{\mathbf{Q}}^{(\ell)} \mathbf{b}}{\partial N_{i,j}^{(\ell)}} &= \sum_{k=1}^{n_{\ell}} \sum_{l=1}^{n_{\ell}} a_k \frac{\partial \tilde{Q}_{k,l}^{(\ell)}}{\partial N_{i,j}^{(\ell)}} b_l \\
&= -\frac{1}{\sigma_N \sqrt{N}} \sum_{k=1}^{n_{\ell}} \sum_{l=1}^{n_{\ell}} a_k [\tilde{Q}_{k,i}^{(\ell)} [\mathbf{T}^{(\ell)\top} \tilde{\mathbf{Q}}^{(\ell)}]_{j,l} + [\tilde{\mathbf{Q}}^{(\ell)} \mathbf{T}^{(\ell)}]_{k,j} \tilde{Q}_{i,l}^{(\ell)}] b_l \\
&= -\frac{1}{\sigma_N \sqrt{N}} [\tilde{\mathbf{Q}}^{(\ell)} \mathbf{a} \mathbf{b}^{\top} \tilde{\mathbf{Q}}^{(\ell)} \mathbf{T}^{(\ell)} + \tilde{\mathbf{Q}}^{(\ell)} \mathbf{b} \mathbf{a}^{\top} \tilde{\mathbf{Q}}^{(\ell)} \mathbf{T}^{(\ell)}]_{i,j}.
\end{align*}
$$

Thus, using the fact that $(a+b)^2 \le 2(a^2+b^2)$, we have,

$$
\begin{align*}
& \sum_{i=1}^{n_{\ell}} \sum_{j=1}^{\prod_{\ell' \neq \ell} n_{\ell'}} \mathbb{E} \left[ \left( \frac{\partial \mathbf{a}^{\top} \tilde{\mathbf{Q}}^{(\ell)} \mathbf{b}}{\partial N_{i,j}^{(\ell)}} \right)^2 \right] \\
& \le \frac{2}{\sigma_N^2 N} \sum_{i=1}^{n_{\ell}} \sum_{j=1}^{\prod_{\ell' \neq \ell} n_{\ell'}} \mathbb{E} \left[ [\tilde{\mathbf{Q}}^{(\ell)} \mathbf{a} \mathbf{b}^{\top} \tilde{\mathbf{Q}}^{(\ell)} \mathbf{T}^{(\ell)}]_{i,j}^2 + [\tilde{\mathbf{Q}}^{(\ell)} \mathbf{b} \mathbf{a}^{\top} \tilde{\mathbf{Q}}^{(\ell)} \mathbf{T}^{(\ell)}]_{i,j}^2 \right] \\
& = \frac{2}{\sigma_N^2 N} \mathbb{E} \left[ \left\| \tilde{\mathbf{Q}}^{(\ell)} \mathbf{a} \mathbf{b}^{\top} \tilde{\mathbf{Q}}^{(\ell)} \mathbf{T}^{(\ell)} \right\|_{\mathrm{F}}^2 + \left\| \tilde{\mathbf{Q}}^{(\ell)} \mathbf{b} \mathbf{a}^{\top} \tilde{\mathbf{Q}}^{(\ell)} \mathbf{T}^{(\ell)} \right\|_{\mathrm{F}}^2 \right] \\
& \le \frac{2}{\sigma_N^2 N} \|\mathbf{a}\|^2 \|\mathbf{b}\|^2 \mathbb{E} \left[ \left\| \tilde{\mathbf{Q}}^{(\ell)} \right\|^2 \left\| \mathbf{T}^{(\ell)} \right\|^2 \right] = \mathcal{O}_{\tilde{z}}(N^{-1})
\end{align*}
$$

where the last line is obtained similarly to the previous case. Hence, $\mathrm{Var}(\mathbf{a}^{\top}\tilde{\mathbf{Q}}^{(\ell)}\mathbf{b}) = \mathcal{O}_{\tilde{z}}(N^{-1})$, which is not enough to apply Lemma 20. We must therefore evaluate the moment of order 4,

$$
\mathbb{E} \left[ \left( \mathbf{a}^{\top} \left( \tilde{\mathbf{Q}}^{(\ell)} - \mathbb{E}[\tilde{\mathbf{Q}}^{(\ell)}] \right) \mathbf{b} \right)^4 \right] = \mathrm{Var} \left( \left( \mathbf{a}^{\top} \left( \tilde{\mathbf{Q}}^{(\ell)} - \mathbb{E}[\tilde{\mathbf{Q}}^{(\ell)}] \right) \mathbf{b} \right)^2 \right) \\
+ \mathbb{E} \left[ \left( \mathbf{a}^{\top} \left( \tilde{\mathbf{Q}}^{(\ell)} - \mathbb{E}[\tilde{\mathbf{Q}}^{(\ell)}] \right) \mathbf{b} \right)^2 \right]^2.
$$

Since we already known that the rightmost term is $O(N^{-2})$, we use the Poincaré-Nash inequality (Lemma 19) on the variance of $(\mathbf{a}^{\top}(\tilde{\mathbf{Q}}^{(\ell)} - \mathbb{E}[\tilde{\mathbf{Q}}^{(\ell)}])\mathbf{b})^2$,

$$
\mathrm{Var} \left( \left( \mathbf{a}^{\top} \left( \tilde{\mathbf{Q}}^{(\ell)} - \mathbb{E}[\tilde{\mathbf{Q}}^{(\ell)}] \right) \mathbf{b} \right)^2 \right) \\
\le \sum_{i,j} \mathbb{E} \left[ \left( \mathbf{a}^{\top} \frac{\partial}{\partial N_{i,j}^{(\ell)}} [\tilde{\mathbf{Q}}^{(\ell)} \mathbf{b} \mathbf{a}^{\top} \tilde{\mathbf{Q}}^{(\ell)} - \tilde{\mathbf{Q}}^{(\ell)} \mathbf{b} \mathbf{a}^{\top} \mathbb{E}[\tilde{\mathbf{Q}}^{(\ell)}] - \mathbb{E}[\tilde{\mathbf{Q}}^{(\ell)}] \mathbf{b} \mathbf{a}^{\top} \tilde{\mathbf{Q}}^{(\ell)}] \mathbf{b} \right)^2 \right].
$$
Similarly, injecting the expression of $\partial \tilde{Q}^{(\ell)}/\partial N_{i,j}^{(\ell)}$, straightforward computations yield

$$
\begin{align*}
\mathrm{Var}((\boldsymbol{a}^{\top}(\tilde{\boldsymbol{Q}}^{(\ell)} - \mathbb{E}[\tilde{\boldsymbol{Q}}^{(\ell)}])\boldsymbol{b})^2) & \\
&\le \sum_{i,j} \mathbb{E}\left[\left(-\frac{2}{\sigma_N \sqrt{N}}(\boldsymbol{a}^{\top}(\tilde{\boldsymbol{Q}}^{(\ell)} - \mathbb{E}[\tilde{\boldsymbol{Q}}^{(\ell)}])\boldsymbol{b})[\tilde{\boldsymbol{Q}}^{(\ell)}(\boldsymbol{a}\boldsymbol{b}^{\top} + \boldsymbol{b}\boldsymbol{a}^{\top})\tilde{\boldsymbol{Q}}^{(\ell)}\boldsymbol{T}^{(\ell)}]_{i,j}\right)^2\right] \\
&\le \frac{8}{\sigma_N^2 N} \sum_{i,j} \mathbb{E}\left[(\boldsymbol{a}^{\top}(\tilde{\boldsymbol{Q}}^{(\ell)} - \mathbb{E}[\tilde{\boldsymbol{Q}}^{(\ell)}])\boldsymbol{b})^2 ([\tilde{\boldsymbol{Q}}^{(\ell)}\boldsymbol{a}\boldsymbol{b}^{\top}\tilde{\boldsymbol{Q}}^{(\ell)}\boldsymbol{T}^{(\ell)}]_{i,j}^2 + [\tilde{\boldsymbol{Q}}^{(\ell)}\boldsymbol{b}\boldsymbol{a}^{\top}\tilde{\boldsymbol{Q}}^{(\ell)}\boldsymbol{T}^{(\ell)}]_{i,j}^2)\right] \\
&\le \frac{16}{\sigma_N^2 N} \|\boldsymbol{a}\|^2 \|\boldsymbol{b}\|^2 \mathbb{E}\left[(\boldsymbol{a}^{\top}(\tilde{\boldsymbol{Q}}^{(\ell)} - \mathbb{E}[\tilde{\boldsymbol{Q}}^{(\ell)}])\boldsymbol{b})^2 \|\tilde{\boldsymbol{Q}}^{(\ell)}\|^4 \|\boldsymbol{T}^{(\ell)}\|^2\right].
\end{align*}
$$

Since we have shown that $\mathbb{E}[(\boldsymbol{a}^{\top}(\tilde{\boldsymbol{Q}}^{(\ell)} - \mathbb{E}[\tilde{\boldsymbol{Q}}^{(\ell)}])\boldsymbol{b})^2] = \mathcal{O}_{\tilde{z}}(N^{-1})$, we have $\mathrm{Var}((\boldsymbol{a}^{\top}(\tilde{\boldsymbol{Q}}^{(\ell)} - \mathbb{E}[\tilde{\boldsymbol{Q}}^{(\ell)}])\boldsymbol{b})^2) = \mathcal{O}_{\tilde{z}}(N^{-2})$; Q.E.D.

## B.4 Expansion of the Mean Empirical Stieltjes Transform

Let $\tilde{\boldsymbol{Q}}_0^{(\ell)}$ denote the resolvent of $\frac{1}{\sigma_N} [\frac{1}{N} \boldsymbol{N}^{(\ell)} \boldsymbol{N}^{(\ell)\top} - \mu_N^{(\ell)} \boldsymbol{I}_{n_\ell}]$. Using the resolvent identity $\boldsymbol{A}^{-1} - \boldsymbol{B}^{-1} = \boldsymbol{A}^{-1}(\boldsymbol{B} - \boldsymbol{A})\boldsymbol{B}^{-1}$, we can see that $\frac{1}{n_\ell} \mathrm{Tr}(\tilde{\boldsymbol{Q}}_0^{(\ell)} - \tilde{\boldsymbol{Q}}^{(\ell)}) = \mathcal{O}_{\tilde{z}}(N^{-1})$ therefore the low-rank perturbation $\boldsymbol{P}^{(\ell)}$ does not change the limiting spectral distribution and we can consider $\boldsymbol{P}^{(\ell)} = \boldsymbol{0}_{n_\ell \times \prod_{\ell' \neq \ell} n_{\ell'}}$ from now on.

Applying $\frac{1}{n_\ell} \mathrm{Tr}$ to Equation (13) and using the fact that $\frac{1}{n_\ell} \mathrm{Tr} \tilde{\boldsymbol{Q}}^{(\ell)2} = \tilde{m}_N^{(\ell)\prime}$ as well as $\mathbb{E}[\tilde{m}_N^{(\ell)2}(\tilde{z})] = \mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})]^2 + \mathrm{Var}(\tilde{m}_N^{(\ell)}(\tilde{z}))$ with $\mathrm{Var}(\tilde{m}_N^{(\ell)}(\tilde{z})) = \mathcal{O}_{\tilde{z}}(N^{-2})$, we find

$$
\begin{align}
\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})]^2 + \tilde{z}\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})] + 1 & \nonumber \\
&= -\frac{1}{\sigma_N} \frac{n_\ell}{N} \mathbb{E}[\tilde{z}\tilde{m}_N^{(\ell)2}(\tilde{z}) + \tilde{m}_N^{(\ell)}(\tilde{z})] - \frac{1}{N} \frac{\mu_N^{(\ell)}}{\sigma_N^2} \mathbb{E}[\tilde{m}_N^{(\ell)\prime}(\tilde{z})] + \mathcal{O}_{\tilde{z}}(N^{-\min(\frac{d}{2}, 2)}) \tag{14}
\end{align}
$$

where the $\mathcal{O}_{\tilde{z}}(N^{-\min(\frac{d}{2}, 2)})$ stems from $\sigma_N = \Theta(N^{\frac{d-2}{2}})$. We know that the Stieltjes transform $m_{\mathrm{SC}}$ of the Wigner semicircle distribution on $[-2, 2]$ satisfies $m_{\mathrm{SC}}^2(\tilde{z}) + \tilde{z}m_{\mathrm{SC}}(\tilde{z}) + 1 = 0$ for all $\tilde{z} \in \mathbb{C} \setminus [-2, 2]$ (see, e.g., Pastur and Shcherbina, 2011, Corollary 2.2.8). Let us subtract this relation to Equation (14) and factorize by $\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})] - m_{\mathrm{SC}}(\tilde{z})$ using the relation $a^2 - b^2 = (a - b)(a + b)$.

$$
\begin{align*}
& (\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})] - m_{\mathrm{SC}}(\tilde{z})) (\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})] + m_{\mathrm{SC}}(\tilde{z}) + \tilde{z}) \\
&= -\frac{1}{\sigma_N} \frac{n_\ell}{N} \mathbb{E}[\tilde{z} \tilde{m}_N^{(\ell)2}(\tilde{z}) + \tilde{m}_N^{(\ell)}(\tilde{z})] - \frac{1}{N} \frac{\mu_N^{(\ell)}}{\sigma_N^2} \mathbb{E}[\tilde{m}_N^{(\ell)\prime}(\tilde{z})] + \mathcal{O}_{\tilde{z}}(N^{-\min(\frac{d}{2}, 2)}).
\end{align*}
$$
Let $g_N^{(\ell)}: \tilde{z} \mapsto \frac{-1}{\tilde{z}+m_{\text{SC}}(\tilde{z})+\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})]}$. Since, $m_{\text{SC}}$ and $\mathbb{E}[\tilde{m}_N^{(\ell)}]$ are Stieltjes transforms$^{21}$, we have $|g_N^{(\ell)}(\tilde{z})| \le |\mathfrak{F}[\tilde{z}+m_{\text{SC}}(\tilde{z})+\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})]]|^{-1} \le |\mathfrak{F}\tilde{z}|^{-1}$ therefore

$$ \mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})] - m_{\text{SC}}(\tilde{z}) = g_N^{(\ell)}(\tilde{z}) \left[ \frac{1}{\sigma_N} \frac{n_\ell}{N} \mathbb{E}[\tilde{z} \tilde{m}_N^{(\ell)2}(\tilde{z}) + \tilde{m}_N^{(\ell)}(\tilde{z})] + \frac{1}{N} \frac{\mu_N^{(\ell)}}{\sigma_N^2} \mathbb{E}[\tilde{m}_N^{(\ell)\prime}(\tilde{z})] \right] + \mathcal{O}_{\tilde{z}}(N^{-\min(\frac{d}{2},2)}) \quad (15) $$

$$ = g_N^{(\ell)}(\tilde{z}) \frac{n_\ell}{N} \mathbb{E}[\tilde{z} \tilde{m}_N^{(\ell)2}(\tilde{z}) + \tilde{m}_N^{(\ell)}(\tilde{z})] + \mathcal{O}_{\tilde{z}}(N^{-1}) \quad (16) $$

$$ = \mathcal{O}_{\tilde{z}}(N^{-\min(\frac{d-2}{2},1)}). \quad (17) $$

Notice that the dominant term in the difference $\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})] - m_{\text{SC}}(\tilde{z})$ differs depending on whether $d=3$, $d=4$ or $d \ge 5$ because $\sigma_N = \Theta(N^{\frac{d-2}{2}})$. The higher $d$, the faster the convergence of the empirical spectral distribution to the semicircle distribution. Indeed, at this point, Equation (17) already shows the pointwise convergence of $\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})]$ to $m_{\text{SC}}$ and therefore the weak convergence of the corresponding distributions (Geronimo and Hill, 2003). Yet, in order to show the confinement of the spectrum below, we need an explicit expansion of $\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})]$ with all the terms not dominated by $N^{-1}$, which we state in the following lemma.

**Lemma 23** $\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})] = m_{\text{SC}}(\tilde{z}) + h_N^{(\ell)}(\tilde{z}) + \mathcal{O}_{\tilde{z}}(N^{-\min(\frac{d}{2},2)})$ with

$$ \begin{gathered} h_N^{(\ell)}(\tilde{z}) = \tau_N^{(\ell)}(\tilde{z}) \left( \frac{1}{\sigma_N} \frac{n_\ell}{N} (\tilde{z} t_N^{(\ell)2}(\tilde{z}) + t_N^{(\ell)}(\tilde{z})) - \frac{1}{N} \frac{\mu_N^{(\ell)}}{\sigma_N^2} \frac{m_{\text{SC}}(\tilde{z})}{\tilde{z} + 2m_{\text{SC}}(\tilde{z})} \right) \\ \text{where } \tau_N^{(\ell)}(\tilde{z}) = \frac{-1}{\tilde{z} + 2m_{\text{SC}}(\tilde{z})} \left( 1 + \frac{1}{\sigma_N} \frac{n_\ell}{N} \frac{\tilde{z} m_{\text{SC}}^2(\tilde{z}) + m_{\text{SC}}(\tilde{z})}{(\tilde{z} + 2m_{\text{SC}}(\tilde{z}))^2} \right) \\ \text{and } t_N^{(\ell)}(\tilde{z}) = m_{\text{SC}}(\tilde{z}) - \frac{1}{\sigma_N} \frac{n_\ell}{N} \frac{\tilde{z} m_{\text{SC}}^2(\tilde{z}) + m_{\text{SC}}(\tilde{z})}{\tilde{z} + 2m_{\text{SC}}(\tilde{z})}. \end{gathered} $$

**Proof** Firstly, with the Poincaré-Nash inequality (Lemma 19), we can show that

$$ \mathrm{Var}(\tilde{m}_N^{(\ell)\prime}(\tilde{z})) \le \sum_{i=1}^{n_\ell} \sum_{j=1}^{\prod_{\ell' \ne \ell} n_{\ell'}} \mathbb{E} \left[ \left| \frac{1}{n_\ell} \frac{\partial \mathrm{Tr} \, \tilde{Q}^{(\ell)2}}{\partial N_{i,j}^{(\ell)}} \right|^2 \right] \le \frac{16 \mathbb{E} \left[ \left\| \tilde{Q}^{(\ell)} \right\|^6 \left\| \mathbf{N}^{(\ell)} \right\|^2 \right]}{n_\ell \sigma_N^2 N^2} = \mathcal{O}_{\tilde{z}}(N^{-2}). $$

Taking $\frac{1}{n_\ell}$ Tr of Equation (13) with $\mathbf{P}^{(\ell)} = \mathbf{0}_{n_\ell \times \prod_{\ell' \ne \ell} n_{\ell'}}$ and differentiating with respect to the complex variable $\tilde{z}$ under $\mathbb{E}$ (which is possible because the integrand can be upper

<footnote n="21">21. $\tilde{m}_N^{(\ell)}$ is obviously a Stieltjes transform and it can be shown that its expectation is also a Stieltjes transform with Herglotz theorem (Weidmann, 1980, Theorem B3).</footnote>
bounded on every compact subset of $C \setminus \mathbb{R})$, we find

$$ 
\begin{aligned}
& \mathbb{E}[(2\tilde{m}_N^{(\ell)}(\tilde{z}) + \tilde{z})\tilde{m}_N^{(\ell)'}(\tilde{z}) + \tilde{m}_N^{(\ell)}(\tilde{z})] \\
& = \mathbb{E}\left[-\frac{1}{\sigma_N}\frac{n_\ell}{N}(\tilde{m}_N^{(\ell)2}(\tilde{z}) + 2\tilde{z}\tilde{m}_N^{(\ell)}(\tilde{z})\tilde{m}_N^{(\ell)'}(\tilde{z}) + \tilde{m}_N^{(\ell)'}(\tilde{z})\right. \\
& \qquad \left. -\frac{1}{N}\left(\frac{\mu_N^{(\ell)}}{\sigma_N^2} + \frac{\tilde{z}}{\sigma_N}\right)\tilde{m}_N^{(\ell)''}(\tilde{z}) - \frac{2\tilde{m}^{(\ell)'}(\tilde{z})}{\sigma_N N}\right]. \quad (18)
\end{aligned}
$$

Then, since $\mathrm{Var} \tilde{m}_N^{(\ell)}(\tilde{z}) = \mathcal{O}_{\tilde{z}}(N^{-2})$ and $\mathrm{Var} \tilde{m}_N^{(\ell)'}(\tilde{z}) = \mathcal{O}_{\tilde{z}}(N^{-2})$, we have $\mathbb{E}[\tilde{m}_N^{(\ell)2}(\tilde{z})] = \mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})]^2 + \mathcal{O}_{\tilde{z}}(N^{-2})$ and $\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})\tilde{m}_N^{(\ell)'}(\tilde{z})] = \mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})]\mathbb{E}[\tilde{m}_N^{(\ell)'}(\tilde{z})] + \mathcal{O}_{\tilde{z}}(N^{-2})$ (using the Cauchy-Schwarz inequality to upper bound the covariance). And, with the fact that $\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})] = m_{\mathrm{SC}}(\tilde{z}) + \mathcal{O}_{\tilde{z}}(N^{-\min(\frac{d-2}{2},1)})$ (Equation (17)), we obtain from Equation (18) that

$$ \mathbb{E}[\tilde{m}_N^{(\ell)'}(\tilde{z})] = \frac{-m_{\mathrm{SC}}(\tilde{z})}{\tilde{z} + 2m_{\mathrm{SC}}(\tilde{z})} + \mathcal{O}_{\tilde{z}}(N^{-\min(\frac{d-2}{2},1)}) $$

since $|\tilde{z} + 2m_{\mathrm{SC}}(\tilde{z})|^{-1} \le |\tilde{z}|^{-1} \le |\Im \tilde{z}|^{-1}$. Moreover, using Equation (17) in Equation (16), we find

$$ \mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})] = m_{\mathrm{SC}}(\tilde{z}) + g_N^{(\ell)}(\tilde{z}) \frac{n_\ell}{N} \left( \tilde{z} m_{\mathrm{SC}}^2(\tilde{z}) + m_{\mathrm{SC}}(\tilde{z}) \right) + \mathcal{O}_{\tilde{z}}(N^{-1}). \quad (19) $$

We can now inject the last two relations into Equation (15):

$$ 
\begin{aligned}
\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})] = m_{\mathrm{SC}}(\tilde{z}) &+ g_N^{(\ell)}(\tilde{z}) \frac{\tilde{z}}{\sigma_N} \frac{n_\ell}{N} \left( m_{\mathrm{SC}}(\tilde{z}) + g_N^{(\ell)}(\tilde{z}) \frac{1}{\sigma_N} \frac{n_\ell}{N} \left( \tilde{z} m_{\mathrm{SC}}^2(\tilde{z}) + m_{\mathrm{SC}}(\tilde{z}) \right) \right)^2 \\
&+ g_N^{(\ell)}(\tilde{z}) \frac{1}{\sigma_N} \frac{n_\ell}{N} \left( m_{\mathrm{SC}}(\tilde{z}) + g_N^{(\ell)}(\tilde{z}) \frac{1}{\sigma_N} \frac{n_\ell}{N} \left( \tilde{z} m_{\mathrm{SC}}^2(\tilde{z}) + m_{\mathrm{SC}}(\tilde{z}) \right) \right) \\
&- g_N^{(\ell)}(\tilde{z}) \frac{1}{N} \frac{\mu_N^{(\ell)}}{\sigma_N^2} \frac{m_{\mathrm{SC}}(\tilde{z})}{\tilde{z} + 2m_{\mathrm{SC}}(\tilde{z})} + \mathcal{O}_{\tilde{z}}(N^{-\min(\frac{d}{2},1)}).
\end{aligned}
\quad (20) $$

We only need to handle the asymptotic behavior of $g_N^{(\ell)}(\tilde{z})$ to conclude the proof. With Equation (17) and the fact that $|\tilde{z} + 2m_{\mathrm{SC}}(\tilde{z})|^{-1} \le |\Im \tilde{z}|^{-1}$, we have $g_N^{(\ell)}(\tilde{z}) = \frac{-1}{\tilde{z} + 2m_{\mathrm{SC}}(\tilde{z})} + \mathcal{O}_{\tilde{z}}(N^{-\min(\frac{d-2}{2},1)})$. We can then use this relation in Equation (19):

$$ \mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})] - m_{\mathrm{SC}}(\tilde{z}) = -\frac{1}{\sigma_N} \frac{n_\ell}{N} \frac{\tilde{z} m_{\mathrm{SC}}^2(\tilde{z}) + m_{\mathrm{SC}}(\tilde{z})}{\tilde{z} + 2m_{\mathrm{SC}}(\tilde{z})} + \mathcal{O}_{\tilde{z}}(N^{-1}). $$
Therefore, we have

$$
\begin{align*}
g_N^{(\ell)}(\tilde{z}) &= \frac{-1}{\tilde{z} + m_{\text{SC}}(\tilde{z}) + \mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})]} \\
&= \frac{-1}{\tilde{z} + 2m_{\text{SC}}(\tilde{z})} \left( 1 + \frac{\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})] - m_{\text{SC}}(\tilde{z})}{\tilde{z} + 2m_{\text{SC}}(\tilde{z})} \right)^{-1} \\
&= \frac{-1}{\tilde{z} + 2m_{\text{SC}}(\tilde{z})} \left( 1 - \frac{\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})] - m_{\text{SC}}(\tilde{z})}{\tilde{z} + 2m_{\text{SC}}(\tilde{z})} + \mathcal{O}_{\tilde{z}}(N^{-\min(d-2,2)}) \right) \\
&= \frac{-1}{\tilde{z} + 2m_{\text{SC}}(\tilde{z})} \left( 1 + \frac{1}{\sigma_N} \frac{n_\ell}{N} \frac{\tilde{z}m_{\text{SC}}^2(\tilde{z}) + m_{\text{SC}}(\tilde{z})}{(\tilde{z} + 2m_{\text{SC}}(\tilde{z}))^2} \right) + \mathcal{O}_{\tilde{z}}(N^{-1}).
\end{align*}
$$

Eventually, Equation (20) becomes

$$
\begin{align*}
\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})] &= m_{\text{SC}}(\tilde{z}) + \tau_N^{(\ell)}(\tilde{z}) \frac{\tilde{z}}{\sigma_N} \frac{n_\ell}{N} \left( m_{\text{SC}}(\tilde{z}) - \frac{1}{\sigma_N} \frac{n_\ell}{N} \frac{\tilde{z} m_{\text{SC}}^2(\tilde{z}) + m_{\text{SC}}(\tilde{z})}{\tilde{z} + 2m_{\text{SC}}(\tilde{z})} \right)^2 \\
&\quad + \tau_N^{(\ell)}(\tilde{z}) \frac{1}{\sigma_N} \frac{n_\ell}{N} \left( m_{\text{SC}}(\tilde{z}) - \frac{1}{\sigma_N} \frac{n_\ell}{N} \frac{\tilde{z} m_{\text{SC}}^2(\tilde{z}) + m_{\text{SC}}(\tilde{z})}{\tilde{z} + 2m_{\text{SC}}(\tilde{z})} \right) \\
&\quad - \tau_N^{(\ell)}(\tilde{z}) \frac{1}{N} \frac{\mu_N^{(\ell)}}{\sigma_N^2} \frac{m_{\text{SC}}(\tilde{z})}{\tilde{z} + 2m_{\text{SC}}(\tilde{z})} + \mathcal{O}_{\tilde{z}}(N^{-\min(\frac{d}{2}, 2)}).
\end{align*}
$$

$$
\text{with } \tau_N^{(\ell)} : \tilde{z} \mapsto \frac{-1}{\tilde{z} + 2m_{\text{SC}}(\tilde{z})} \left( 1 + \frac{1}{\sigma_N} \frac{n_\ell}{N} \frac{\tilde{z}m_{\text{SC}}^2(\tilde{z}) + m_{\text{SC}}(\tilde{z})}{(\tilde{z} + 2m_{\text{SC}}(\tilde{z}))^2} \right).
$$

<figure><img src="image_6.png" /></figure>

B.5 Confinement of the Spectrum

We have just found, in Lemma 23, the following expansion: $\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})] = m_{\text{SC}}(\tilde{z}) + h_N^{(\ell)}(\tilde{z}) + \mathcal{O}_{\tilde{z}}(N^{-\min(\frac{d}{2},2)}).$ Note that this shows that the limiting spectral distribution of $\frac{1}{\sigma_N}[\frac{1}{N}N^{(\ell)}N^{(\ell)\top} - \mu_N^{(\ell)}I_{n_\ell}]$ (and therefore that of $\frac{1}{\sigma_N}[T^{(\ell)}T^{(\ell)\top} - \mu_N^{(\ell)}I_{n_\ell}]$) is the semicircle distribution on $[-2, 2]$. We will now use this expansion to prove that, almost surely, no eigenvalue of $\frac{1}{\sigma_N}[\frac{1}{N}N^{(\ell)}N^{(\ell)\top} - \mu_N^{(\ell)}I_{n_\ell}]$ stays outside $[-2, 2]$ as $N \to +\infty$. That is, for all $\varepsilon > 0$, there exists an integer $N_0$ such that, for all $\lambda \in \text{Sp}(\frac{1}{\sigma_N}[\frac{1}{N}N^{(\ell)}N^{(\ell)\top} - \mu_N^{(\ell)}I_{n_\ell}])$, $\text{Dist}(\lambda, \text{Supp } \mu_{\text{MP}}) \le \varepsilon$ as soon as $N \ge N_0$.

Let $\varepsilon > 0$, $\varphi : \mathbb{R} \mapsto [0, 1]$ be an infinitely differentiable function which equals 1 on $[-2, 2]$ and 0 on $\mathbb{R} \setminus [-2 - \varepsilon, 2 + \varepsilon]$ and $\psi = 1 - \varphi$. We want to show that

$$
\mathrm{Tr} \left( \psi \left( \frac{1}{\sigma_N} \left[ \frac{1}{N} \mathbf{N}^{(\ell)} \mathbf{N}^{(\ell)\top} - \mu_N^{(\ell)} \mathbf{I}_{n_\ell} \right] \right) \right) \to 0 \quad \text{almost surely as } N \to +\infty.
$$
First of all, we show the convergence in mean with the Hellfer-Sjöstrand formula (Proposition 21).

$$ 
\begin{aligned}
& \mathbb{E} \left[ \frac{1}{n_{\ell}} \operatorname{Tr} \left( \varphi \left( \frac{1}{\sigma_N} \left[ \frac{1}{N} \mathbf{N}^{(\ell)} \mathbf{N}^{(\ell)\top} - \mu_N^{(\ell)} \mathbf{I}_{n_{\ell}} \right] \right) \right) \right] \\
&= \frac{2}{\pi} \Re \int_{\mathbb{C}^+} \frac{\partial \Phi_q[\varphi]}{\partial \tilde{z}}(\tilde{z}) \mathbb{E} \left[ \tilde{m}_N^{(\ell)}(\tilde{z}) \right] \mathrm{d}\tilde{z} \\
&= \frac{2}{\pi} \Re \int_{\mathbb{C}^+} \frac{\partial \Phi_q[\varphi]}{\partial \tilde{z}}(\tilde{z}) m_{\mathrm{SC}}(\tilde{z}) \mathrm{d}\tilde{z} + \frac{2}{\pi} \Re \int_{\mathbb{C}^+} \frac{\partial \Phi_q[\varphi]}{\partial \tilde{z}}(\tilde{z}) h_N^{(\ell)}(\tilde{z}) \mathrm{d}\tilde{z} \\
&\quad + \frac{2}{\pi} \Re \int_{\mathbb{C}^+} \frac{\partial \Phi_q[\varphi]}{\partial \tilde{z}}(\tilde{z}) \times \mathcal{O}_{\tilde{z}}(N^{-\min(\frac{d}{2}, 2)}) \mathrm{d}\tilde{z}.
\end{aligned}
$$

where we have used the expression of $\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})]$ given by Lemma 23. The first integral is $\int_{\mathbb{R}} \varphi \, d\mu_{\mathrm{SC}} = 1$ while the last one is $\mathcal{O}(N^{-\min(\frac{d}{2}, 2)})$ with $q$ chosen sufficiently large so that $\frac{\partial \Phi_q[\varphi]}{\partial \tilde{z}}(\tilde{z})$ cancels the divergence of $\mathcal{O}_{\tilde{z}}(N^{-\min(\frac{d}{2}, 2)})$ near the real axis. In order to evaluate the second integral, we perform an integration by parts:

$$ 
\begin{aligned}
& \frac{2}{\pi} \Re \int_{\mathbb{C}^+} \frac{\partial \Phi_q[\varphi]}{\partial \tilde{z}}(\tilde{z}) h_N^{(\ell)}(\tilde{z}) \mathrm{d}\tilde{z} \\
&= \frac{2}{\pi} \Re \left[ \frac{1}{2} \int_0^{+\infty} \left( \int_{-\infty}^{+\infty} \frac{\partial \Phi_q[\varphi]}{\partial x}(x + iy) h_N^{(\ell)}(x + iy) \mathrm{d}x \right) \mathrm{d}y \right. \\
&\qquad \left. + \frac{i}{2} \int_{-\infty}^{+\infty} \left( \int_0^{+\infty} \frac{\partial \Phi_q[\varphi]}{\partial y}(x + iy) h_N^{(\ell)}(x + iy) \mathrm{d}y \right) \mathrm{d}x \right] \\
&= \frac{2}{\pi} \Re \left[ \frac{-i}{2} \int_{\mathbb{R}} \lim_{y \downarrow 0} \left\{ \Phi_q[\varphi](x + iy) h_N^{(\ell)}(x + iy) \right\} \mathrm{d}x - \int_{\mathbb{C}^+} \Phi_q[\varphi](\tilde{z}) \frac{\partial h}{\partial \tilde{z}}(\tilde{z}) \mathrm{d}\tilde{z} \right].
\end{aligned}
$$

Since $h_N^{(\ell)}$ is an analytic function (as sums and products of Stieltjes transforms), we have $\partial h/\partial \tilde{z} = 0$ by the Cauchy-Riemann equations. Moreover, we have $\lim_{y \downarrow 0} \Phi_q[\varphi](x + iy) = \varphi(x)$ and, from the definition of $h_N^{(\ell)}$, $\lim_{y \downarrow 0} \Im[h_N^{(\ell)}(x + iy)] = 0$ for all $x \in \mathbb{R} \setminus [-2, 2]$. Therefore,

$$ \frac{2}{\pi} \Re \int_{\mathbb{C}^+} \frac{\partial \Phi_q[\varphi]}{\partial \tilde{z}}(\tilde{z}) h_N^{(\ell)}(\tilde{z}) \mathrm{d}\tilde{z} = \lim_{y \downarrow 0} \frac{1}{\pi} \int_{-2}^{2} \Im[h_N^{(\ell)}(x + iy)] \mathrm{d}x. $$

We use Lemma 22 to show that this integral equals $\lim_{y \to +\infty} -iyh_N^{(\ell)}(iy) = 0$. The function $h_N^{(\ell)}$ is analytic on $\mathbb{C} \setminus [-2, 2]$, $\lim_{|\tilde{z}| \to +\infty} h(\tilde{z}) = 0$ and $h_N^{(\ell)}(\tilde{z}) = \overline{h(\tilde{z})}$ for all $\tilde{z} \in \mathbb{C} \setminus [-2, 2]$. Thus, we just need to show that there exist an integer $n_0$ and a constant $C > 0$ such that $|h_N^{(\ell)}(\tilde{z})| \le C \max(\mathrm{Dist}(\tilde{z}, [-2, 2])^{-n_0}, 1)$ for all $\tilde{z} \in \mathbb{C} \setminus [-2, 2]$. Firstly, we find an upper bound for $t_N^{(\ell)}(\tilde{z})$. Since $\tilde{z} \mapsto \frac{-1}{\tilde{z}+2m_{\mathrm{SC}}(\tilde{z})}$ is the Stieltjes transform of a probability measure on $[-2, 2]$ and $-(\tilde{z}m_{\mathrm{SC}}(\tilde{z}) + 1) = m_{\mathrm{SC}}^2(\tilde{z})$, we have

$$ 
\begin{aligned}
|t_N^{(\ell)}(\tilde{z})| &= |m_{\mathrm{SC}}(\tilde{z})| \left| 1 - \frac{1}{\sigma_N} \frac{n_\ell}{N} \frac{\tilde{z} m_{\mathrm{SC}}(\tilde{z}) + 1}{\tilde{z} + 2m_{\mathrm{SC}}(\tilde{z})} \right| \\
&\le \frac{1}{\mathrm{Dist}(\tilde{z}, [-2, 2])} \left( 1 + \frac{1}{\sigma_N} \frac{n_\ell}{N} \frac{1}{\mathrm{Dist}(\tilde{z}, [-2, 2])^3} \right) \\
&\le \left( 1 + \frac{1}{\sigma_N} \frac{n_\ell}{N} \right) \max(\mathrm{Dist}(\tilde{z}, [-2, 2])^{-4}, 1).
\end{aligned}
$$
Moreover, since $|\tilde{z}| \le \text{Dist}(\tilde{z}, [-2, 2]) + 2$, we also have

$$ \begin{aligned} |\tilde{z}t_N^{(\ell)}(\tilde{z})| &\le \left(1 + \frac{2}{\text{Dist}(\tilde{z}, [-2, 2])}\right) \left(1 + \frac{1}{\sigma_N} \frac{n_\ell}{N} \frac{1}{\text{Dist}(\tilde{z}, [-2, 2])^3}\right) \\ &\le 3 \left(1 + \frac{1}{\sigma_N} \frac{n_\ell}{N}\right) \max(\text{Dist}(\tilde{z}, [-2, 2])^{-4}, 1). \end{aligned} $$

Similarly,

$$ \begin{aligned} |\tau_N^{(\ell)}(\tilde{z})| &= \left| \frac{-1}{\tilde{z} + 2m_{\text{SC}}(\tilde{z})} \right| \left| 1 + \frac{1}{\sigma_N} \frac{n_\ell}{N} \frac{m_{\text{SC}}(\tilde{z})(\tilde{z}m_{\text{SC}}(\tilde{z}) + 1)}{(\tilde{z} + 2m_{\text{SC}}(\tilde{z}))^2} \right| \\ &\le \left(1 + \frac{1}{\sigma_N} \frac{n_\ell}{N}\right) \max(\text{Dist}(\tilde{z}, [-2, 2])^{-6}, 1). \end{aligned} $$

Hence, we can upper bound $|h_N^{(\ell)}(\tilde{z})|$:

$$ \begin{aligned} |h_N^{(\ell)}(\tilde{z})| &\le |\tau_N^{(\ell)}(\tilde{z})| \left( \frac{1}{\sigma_N} \frac{n_\ell}{N} |t_N^{(\ell)}(\tilde{z})| \left( |\tilde{z}t_N^{(\ell)}(\tilde{z})| + 1 \right) + \frac{1}{N} \frac{\mu_N^{(\ell)}}{\sigma_N^2} \frac{|m_{\text{SC}}(\tilde{z})|}{|\tilde{z} + 2m_{\text{SC}}(\tilde{z})|} \right) \\ &\le \left(1 + \frac{1}{\sigma_N} \frac{n_\ell}{N}\right) \left[ \frac{1}{\sigma_N} \frac{n_\ell}{N} \left(1 + \frac{1}{\sigma_N} \frac{n_\ell}{N}\right) \left[ 3\left(1 + \frac{1}{\sigma_N} \frac{n_\ell}{N}\right) + 1 \right] + \frac{1}{N} \frac{\mu_N^{(\ell)}}{\sigma_N^2} \right] \\ &\qquad \times \max(\text{Dist}(\tilde{z}, [-2, 2])^{-14}, 1). \end{aligned} $$

Therefore we can conclude that $\frac{2}{\pi} \Re \int_{\mathbb{C}+} \frac{\partial \Phi_q[\varphi]}{\partial \tilde{z}}(\tilde{z}) h_N^{(\ell)}(\tilde{z}) \, d\tilde{z} = 0$ (Lemma 22) and

$$ \begin{aligned} \mathbb{E} \left[ \frac{1}{n_\ell} \text{Tr} \left( \varphi \left( \frac{1}{\sigma_N} \left[ \frac{1}{N} \mathbf{N}^{(\ell)} \mathbf{N}^{(\ell)\top} - \mu_N^{(\ell)} \mathbf{I}_{n_\ell} \right] \right) \right) \right] &= 1 + \mathcal{O}(N^{-\min(\frac{d}{2}, 2)}), \\ \text{i.e.,} \quad \mathbb{E} \left[ \text{Tr} \left( \psi \left( \frac{1}{\sigma_N} \left[ \frac{1}{N} \mathbf{N}^{(\ell)} \mathbf{N}^{(\ell)\top} - \mu_N^{(\ell)} \mathbf{I}_{n_\ell} \right] \right) \right) \right] &= \mathcal{O}(N^{-\min(\frac{d-2}{2}, 1)}). \end{aligned} $$

Secondly, we prove the almost sure convergence of $\text{Tr}(\psi(\frac{1}{\sigma_N}[\frac{1}{N}\mathbf{N}^{(\ell)}\mathbf{N}^{(\ell)\top} - \mu_N^{(\ell)}\mathbf{I}_{n_\ell}]))$ to 0 by showing that its variance is $\mathcal{O}(N^{-\min(\frac{d}{2}, 2)})$ (and Lemma 20 implies the result). With
the Poincaré-Nash inequality (Lemma 19), we have

$$
\begin{align*}
& \mathrm{Var}\left(\mathrm{Tr}\left(\psi\left(\frac{1}{\sigma_N}\left[\frac{1}{N}\mathbf{N}^{(\ell)}\mathbf{N}^{(\ell)\top} - \mu_N^{(\ell)}\mathbf{I}_{n_\ell}\right]\right)\right)\right) \\
&= \mathrm{Var}\left(\mathrm{Tr}\left(\varphi\left(\frac{1}{\sigma_N}\left[\frac{1}{N}\mathbf{N}^{(\ell)}\mathbf{N}^{(\ell)\top} - \mu_N^{(\ell)}\mathbf{I}_{n_\ell}\right]\right)\right)\right) \\
&\le \sum_{i=1}^{n_\ell} \sum_{j=1}^{\prod_{\ell' \ne \ell} n_{\ell'}} \mathbb{E}\left[\left|\frac{\partial \mathrm{Tr}\left(\varphi\left(\frac{1}{\sigma_N}\left[\frac{1}{N}\mathbf{N}^{(\ell)}\mathbf{N}^{(\ell)\top} - \mu_N^{(\ell)}\mathbf{I}_{n_\ell}\right]\right)\right)}{\partial N_{i,j}^{(\ell)}}\right|^2\right] \\
&= \sum_{i=1}^{n_\ell} \sum_{j=1}^{\prod_{\ell' \ne \ell} n_{\ell'}} \mathbb{E}\left[\left|\mathrm{Tr}\left(\varphi'\left(\frac{1}{\sigma_N}\left[\frac{1}{N}\mathbf{N}^{(\ell)}\mathbf{N}^{(\ell)\top} - \mu_N^{(\ell)}\mathbf{I}_{n_\ell}\right]\right)\frac{\partial}{\partial N_{i,j}^{(\ell)}}\left[\frac{1}{\sigma_N}\frac{1}{N}\mathbf{N}^{(\ell)}\mathbf{N}^{(\ell)\top}\right]\right)\right|^2\right] \\
&= \frac{1}{\sigma_N^2}\frac{1}{N^2}\sum_{i=1}^{n_\ell} \sum_{j=1}^{\prod_{\ell' \ne \ell} n_{\ell'}} \mathbb{E}\left[\left|\left[2\varphi'\left(\frac{1}{\sigma_N}\left[\frac{1}{N}\mathbf{N}^{(\ell)}\mathbf{N}^{(\ell)\top} - \mu_N^{(\ell)}\mathbf{I}_{n_\ell}\right]\right)\mathbf{N}^{(\ell)}\right]_{i,j}\right|^2\right] \\
&= \frac{1}{\sigma_N}\frac{4}{N}\mathbb{E}\left[\mathrm{Tr}\, u\left(\frac{1}{\sigma_N}\left[\frac{1}{N}\mathbf{N}^{(\ell)}\mathbf{N}^{(\ell)\top} - \mu_N^{(\ell)}\mathbf{I}_{n_\ell}\right]\right)\right] \\
&\quad + \frac{4}{N}\frac{\mu_N^{(\ell)}}{\sigma_N^2}\mathbb{E}\left[\mathrm{Tr}\, \varphi'^2\left(\frac{1}{\sigma_N}\left[\frac{1}{N}\mathbf{N}^{(\ell)}\mathbf{N}^{(\ell)\top} - \mu_N^{(\ell)}\mathbf{I}_{n_\ell}\right]\right)\right]
\end{align*}
$$

where $u: x \mapsto x\varphi'^2(x)$ and $\varphi'^2$ are infinitely differentiable functions with compact support which equal 0 on $[-2, 2]$. Hence, applying similarly the Helffer-Sjöstrand formula (Proposition 21), we find

$$
\begin{align*}
\frac{1}{\sigma_N} \frac{4}{N} \mathbb{E} \left[ \mathrm{Tr} \, u \left( \frac{1}{\sigma_N} \left[ \frac{1}{N} \mathbf{N}^{(\ell)} \mathbf{N}^{(\ell)\top} - \mu_N^{(\ell)} \mathbf{I}_{n_\ell} \right] \right) \right]
&= \frac{4}{\sigma_N} \frac{n_\ell}{N} \frac{2}{\pi} \Re \int_{\mathbb{C}^+} \frac{\partial \Phi_q[u]}{\partial \tilde{z}}(\tilde{z}) \mathbb{E} \left[ \tilde{m}_N^{(\ell)}(\tilde{z}) \right] \, d\tilde{z} \\
&= \mathcal{O} \left( N^{-\min(d-1, \frac{d+2}{2})} \right)
\end{align*}
$$

and

$$
\begin{align*}
& \frac{4}{N} \frac{\mu_N^{(\ell)}}{\sigma_N^2} \mathbb{E} \left[ \mathrm{Tr} \, \varphi'^2 \left( \frac{1}{\sigma_N} \left[ \frac{1}{N} \mathbf{N}^{(\ell)} \mathbf{N}^{(\ell)\top} - \mu_N^{(\ell)} \mathbf{I}_{n_\ell} \right] \right) \right] \\
&= 4 \frac{n_\ell}{N} \frac{\mu_N^{(\ell)}}{\sigma_N^2} \frac{2}{\pi} \Re \int_{\mathbb{C}^+} \frac{\partial \Phi_q[\varphi']}{\partial \tilde{z}}(\tilde{z}) \mathbb{E} \left[ \tilde{m}_N^{(\ell)}(\tilde{z}) \right] \, d\tilde{z} = \mathcal{O} \left( N^{-\min(\frac{d}{2}, 2)} \right)
\end{align*}
$$

for $q$ chosen sufficiently large. Thus,

$$
\mathrm{Var}\left(\mathrm{Tr}\left(\psi\left(\frac{1}{\sigma_N}\left[\frac{1}{N}\mathbf{N}^{(\ell)}\mathbf{N}^{(\ell)\top} - \mu_N^{(\ell)}\mathbf{I}_{n_\ell}\right]\right)\right)\right) = \mathcal{O}\left(N^{-\min\left(\frac{d}{2}, 2\right)}\right)
$$

and we can conclude on the almost sure convergence with Lemma 20.
## B.6 Deterministic Equivalent

With the rescaling $(z, Q^{(\ell)}(z)) \sim (\tilde{z}, \tilde{Q}^{(\ell)}(\tilde{z}))$, Equation (10) becomes

$$ \frac{1}{\sigma_N} \mathbb{E} [\mathbf{P}^{(\ell)} \mathbf{T}^{(\ell)\top} \tilde{\mathbf{Q}}^{(\ell)}] = \frac{1}{\sigma_N} \mathbb{E} [\mathbf{P}^{(\ell)} \mathbf{P}^{(\ell)\top} \tilde{\mathbf{Q}}^{(\ell)}] \\ - \frac{1}{\sigma_N^2} \mathbb{E} \left[ \frac{n_\ell}{N} \tilde{m}_N^{(\ell)}(\tilde{z}) \mathbf{P}^{(\ell)} \mathbf{T}^{(\ell)\top} \tilde{\mathbf{Q}}^{(\ell)} + \frac{1}{N} \mathbf{P}^{(\ell)} \mathbf{T}^{(\ell)\top} \tilde{\mathbf{Q}}^{(\ell)2} \right] $$

where $\left\| \frac{1}{\sigma_N^2} \mathbb{E} \left[ \frac{n_\ell}{N} \tilde{m}_N^{(\ell)}(\tilde{z}) \mathbf{P}^{(\ell)} \mathbf{T}^{(\ell)\top} \tilde{\mathbf{Q}}^{(\ell)} + \frac{1}{N} \mathbf{P}^{(\ell)} \mathbf{T}^{(\ell)\top} \tilde{\mathbf{Q}}^{(\ell)2} \right] \right\| \to 0$ as $N \to +\infty$ since $\| \mathbf{P}^{(\ell)} \| = \mathcal{O}(N^{\frac{d-2}{4}})$ and $\| \mathbf{T}^{(\ell)} \| = \mathcal{O}(N^{\frac{d-2}{2}})$. Hence, with Equation (13) and Equation (17), we have

$$ \left\| m_{\mathrm{SC}}(\tilde{z}) \mathbb{E}[\tilde{\mathbf{Q}}^{(\ell)}] + \tilde{z} \mathbb{E}[\tilde{\mathbf{Q}}^{(\ell)}] + \mathbf{I}_{n_\ell} - \frac{1}{\sigma_N} \mathbf{P}^{(\ell)} \mathbf{P}^{(\ell)\top} \mathbb{E}[\tilde{\mathbf{Q}}^{(\ell)}] \right\| \xrightarrow[N \to +\infty]{} 0 $$

and, since $m_{\mathrm{SC}}(\tilde{z}) + \tilde{z} = \frac{-1}{m_{\mathrm{SC}}(\tilde{z})}$, we can define the following deterministic equivalent (Definition 3):

$$ \tilde{\mathbf{Q}}^{(\ell)}(\tilde{z}) \longleftrightarrow \bar{\mathbf{Q}}^{(\ell)}(\tilde{z}) \stackrel{\text{def}}{=} \left( \frac{1}{\sigma_N} \mathbf{P}^{(\ell)} \mathbf{P}^{(\ell)\top} + \frac{1}{m_{\mathrm{SC}}(\tilde{z})} \mathbf{I}_{n_\ell} \right)^{-1}. $$

## Appendix C. Proof of Theorem 9

Recall that $\tilde{\mathbf{Q}}^{(\ell)}$ is the resolvent of $\frac{1}{\sigma_N} [\mathbf{T}^{(\ell)} \mathbf{T}^{(\ell)\top} - \mu_N^{(\ell)} \mathbf{I}_{n_\ell}]$ while $\tilde{\mathbf{Q}}_0^{(\ell)}$ denotes the resolvent of the same model without signal, $\frac{1}{\sigma_N} [\frac{1}{N} \mathbf{N}^{(\ell)} \mathbf{N}^{(\ell)\top} - \mu_N^{(\ell)} \mathbf{I}_{n_\ell}]$.

### C.1 Convergence of Bilinear Forms

First of all, we must show the almost sure convergence $\mathbf{a}^\top \tilde{\mathbf{Q}}_0^{(\ell)} \mathbf{b} - m_{\mathrm{SC}}(\tilde{z}) \langle \mathbf{a}, \mathbf{b} \rangle \to 0$ for all bounded (sequences of) vectors $\mathbf{a}, \mathbf{b} \in \mathbb{R}^{n_\ell}$. Given the concentration result proven in Section B.3, we just need to show that $\mathbf{a}^\top \mathbb{E}[\tilde{\mathbf{Q}}_0^{(\ell)}] \mathbf{b} - m_{\mathrm{SC}}(\tilde{z}) \langle \mathbf{a}, \mathbf{b} \rangle \to 0$ as $N \to +\infty$.

Let us multiply Equation (13) when $\mathbf{P}^{(\ell)} = \mathbf{0}_{n_\ell} \times \prod_{\ell' \neq \ell} n_{\ell'}$ by $\mathbf{a}^\top$ on the left and $\mathbf{b}$ on the right.

$$ \mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z}) \mathbf{a}^\top \tilde{\mathbf{Q}}_0^{(\ell)} \mathbf{b}] + \tilde{z} \mathbb{E}[\mathbf{a}^\top \tilde{\mathbf{Q}}_0^{(\ell)} \mathbf{b}] + \langle \mathbf{a}, \mathbf{b} \rangle = \mathcal{O}_z(N^{-\min(\frac{d-2}{2},1)}). $$

Then, using the fact that $\mathbb{E}[\tilde{m}_N^{(\ell)}(\tilde{z})] = m_{\mathrm{SC}}(\tilde{z}) + \mathcal{O}_{\tilde{z}}(N^{-\min(\frac{d-2}{2},1)})$ and $m_{\mathrm{SC}}(\tilde{z}) + \tilde{z} = \frac{-1}{m_{\mathrm{SC}}(\tilde{z})}$, we obtain the desired result: $\mathbb{E}[\mathbf{a}^\top \tilde{\mathbf{Q}}_0^{(\ell)} \mathbf{b}] = m_{\mathrm{SC}}(\tilde{z}) \langle \mathbf{a}, \mathbf{b} \rangle + \mathcal{O}_{\tilde{z}}(N^{-\min(\frac{d-2}{2},1)})$.

### C.2 Isolated Eigenvalues

We seek eigenvalues of $\frac{1}{\sigma_N} [\mathbf{T}^{(\ell)} \mathbf{T}^{(\ell)\top} - \mu_N^{(\ell)} \mathbf{I}_{n_\ell}]$ which stay outside the support of the semi-circle distribution $[-2, 2]$. That is, we seek $\tilde{\xi} \in \mathbb{R} \setminus [-2, 2]$ such that

$$ \det\left(\frac{1}{\sigma_N} [\mathbf{T}^{(\ell)} \mathbf{T}^{(\ell)\top} - \mu_N^{(\ell)} \mathbf{I}_{n_\ell}] - \tilde{\xi} \mathbf{I}_{n_\ell}\right) = 0. $$
Using the expansion $T^{(\ell)} = P^{(\ell)} + \frac{1}{\sqrt{N}}N^{(\ell)}$, this is equivalent to

$$ \det\left(\frac{1}{\sigma_N}\left(P^{(\ell)}P^{(\ell)\top} + \frac{1}{\sqrt{N}}P^{(\ell)}N^{(\ell)\top} + \frac{1}{\sqrt{N}}N^{(\ell)}P^{(\ell)\top}\right)\tilde{Q}_0^{(\ell)}(\tilde{\xi}) + I_{n_\ell}\right) \times \det\left(\frac{1}{\sigma_N}\left[\frac{1}{N}N^{(\ell)}N^{(\ell)\top} - \mu_N^{(\ell)}I_{n_\ell}\right] - \tilde{\xi}I_{n_\ell}\right) = 0 $$

where the second determinant is non-zero for $N$ large enough from the confinement of the spectrum proven in Section B.5. Then, we know that $\mathcal{P} = [\mathcal{H}; X^{(1)}, \dots, X^{(d)}]$ therefore $P^{(\ell)} = X^{(\ell)}L^{(\ell)}$ with $L^{(\ell)} = H^{(\ell)}(\otimes_{\ell' \neq \ell} X^{(\ell)\top})$ and we can write

$$ \frac{1}{\sigma_N}\left(P^{(\ell)}P^{(\ell)\top} + \frac{1}{\sqrt{N}}P^{(\ell)}N^{(\ell)\top} + \frac{1}{\sqrt{N}}N^{(\ell)}P^{(\ell)\top}\right) $$

as the matrix product $[X^{(\ell)} X^{(\ell)} \frac{1}{\sigma_N \sqrt{N}} N^{(\ell)} L^{(\ell)\top}] \left[ \frac{1}{\sigma_N H^{(\ell)} H^{(\ell)\top} X^{(\ell)\top}} \frac{1}{\sigma_N \sqrt{N}} L^{(\ell)} N^{(\ell)\top} \frac{1}{X^{(\ell)\top}} \right]$. Thus, with Sylvester's identity $(\det(I_n + AB) = \det(I_k + BA))$, we are left to evaluate a $3r_\ell \times 3r_\ell$ determinant:

$$ \det(I_{3r_\ell} + M) = 0 $$

where $M$ is the following matrix

$$ \left[ \begin{array}{ccc} \frac{1}{\sigma_N H^{(\ell)} H^{(\ell)\top} X^{(\ell)\top} \tilde{Q}_0^{(\ell)}(\tilde{\xi}) X^{(\ell)}} & \frac{1}{\sigma_N H^{(\ell)} H^{(\ell)\top} X^{(\ell)\top} \tilde{Q}_0^{(\ell)}(\tilde{\xi}) X^{(\ell)}} & \frac{1}{\sigma_N^2 \sqrt{N}} H^{(\ell)} H^{(\ell)\top} X^{(\ell)\top} \tilde{Q}_0^{(\ell)}(\tilde{\xi}) N^{(\ell)} L^{(\ell)\top} \\ \frac{1}{\sigma_N \sqrt{N}} L^{(\ell)} N^{(\ell)\top} \tilde{Q}_0^{(\ell)}(\tilde{\xi}) X^{(\ell)} & \frac{1}{\sigma_N \sqrt{N}} L^{(\ell)} N^{(\ell)\top} \tilde{Q}_0^{(\ell)}(\tilde{\xi}) X^{(\ell)} & \frac{1}{\sigma_N^2 \sqrt{N}} L^{(\ell)} N^{(\ell)\top} \tilde{Q}_0^{(\ell)}(\tilde{\xi}) N^{(\ell)} L^{(\ell)\top} \\ X^{(\ell)\top} \tilde{Q}_0^{(\ell)}(\tilde{\xi}) X^{(\ell)} & X^{(\ell)\top} \tilde{Q}_0^{(\ell)}(\tilde{\xi}) X^{(\ell)} & \frac{1}{\sigma_N \sqrt{N}} X^{(\ell)\top} \tilde{Q}_0^{(\ell)}(\tilde{\xi}) N^{(\ell)} L^{(\ell)\top} \end{array} \right] $$

From the convergence of bilinear forms and the orthonormality of the columns of $X^{(\ell)}$, we have $X^{(\ell)\top} \tilde{Q}_0^{(\ell)}(\tilde{\xi}) X^{(\ell)} \to m_{\text{SC}}(\tilde{\xi}) I_{r_\ell}$ almost surely. Moreover, we can see that the matrix $\frac{1}{\sigma_N \sqrt{N}} L^{(\ell)} N^{(\ell)\top} \tilde{Q}_0^{(\ell)}(\tilde{\xi}) X^{(\ell)}$ vanishes almost surely as $N \to +\infty$ since $\|L^{(\ell)}\| = \mathcal{O}(N^{d-2})$ and $\|N^{(\ell)}a\| = \mathcal{O}(\sqrt{N})$ almost surely$^{22}$ for all bounded (sequences of) vectors

<footnote n="22">22. This fact is not so easy to see (note that it does not depend on $d$!). If we naively upper bound $\|N^{(\ell)}a\|^2$ by $\|N^{(\ell)}\|^2 \|a\|^2 = \mathcal{O}(N^{d-2})$, we do not find the desired upper bound. Instead, we can remark that $\|N^{(\ell)}a\|^2 = a^\top N^{(\ell)\top} N^{(\ell)} a = a^\top V D V^\top a$ where $D = \text{Diag}(\lambda_1(N^{(\ell)\top} N^{(\ell)}), \dots, \lambda_{n_\ell}(N^{(\ell)\top} N^{(\ell)}))$ (we assume $N$ large enough so that $\prod_{\ell' \neq \ell} n_{\ell'} > n_\ell$ thus $\lambda_i(N^{(\ell)\top} N^{(\ell)}) = 0$ for $i > n_\ell$) and $V$ follows a uniform distribution on the Stiefel manifold $V_{n_\ell}(\mathbb{R} \prod_{\ell' \neq \ell} n_{\ell'})$ (Chikuse, 2003, Theorem 2.2.1). Therefore $\frac{1}{N}\|N^{(\ell)}a\|^2 \le (\mu_N^{(\ell)} + 2\sigma_N) a^\top V V^\top a$. Without loss of generality, we can assume that $a_i = \|a\|\delta_{1,i}$ (replace $a$ and $V$ by $Oa$ and $OV$ for a well-chosen orthogonal matrix $O$). From Mardia and Khatri (1977), we know that $[VV^\top]_{1,1}$ follows a beta distribution with parameters $\frac{n_\ell}{2}, \frac{\prod_{\ell' \neq \ell} n_{\ell'} - n_\ell}{2}$ so its moments are given by $\mathbb{E}[[VV^\top]_{1,1}] = \prod_{r=0}^{k-1} \frac{n_\ell + 2r}{\prod_{\ell' \neq \ell} n_{\ell'} + 2r}$ for all $k \ge 1$. This is enough to see that $(\mu_N^{(\ell)} + 2\sigma_N)\mathbb{E}[[VV^\top]_{1,1}] = (\mu_N^{(\ell)} + 2\sigma_N) \frac{n_\ell}{\prod_{\ell' \neq \ell} n_{\ell'}} = \mathcal{O}(1)$ and $(\mu_N^{(\ell)} + 2\sigma_N)^4 \mathbb{E}[[VV^\top]_{1,1} - \mathbb{E}[[VV^\top]_{1,1}]]^4] = \mathcal{O}(N^{-2})$, whence the almost sure statement $\frac{1}{N}\|N^{(\ell)}a\|^2 = \mathcal{O}(1)$.</footnote>
$a \in \mathbb{R}^{\prod_{l' \neq l} n_{l'}}$. The only term which remains to evaluate is the block (2,3):

$$ \frac{1}{\sigma_N^2 N} \| \mathbf{L}^{(\ell)} \mathbf{N}^{(\ell)\top} \tilde{\mathbf{Q}}_0^{(\ell)}(\tilde{\xi}) \mathbf{N}^{(\ell)} \mathbf{L}^{(\ell)\top} \| \\ \le \frac{1}{\sigma_N} \frac{\| \mathbf{H}^{(\ell)} \|^2}{\sigma_N} \frac{\| \mathbf{N}^{(\ell)} (\otimes_{l' \neq l} \mathbf{X}^{(l')}) \|^2}{N} \| \tilde{\mathbf{Q}}_0^{(\ell)}(\tilde{\xi}) \| \xrightarrow[N \to +\infty]{\text{a.s.}} 0. $$

Eventually, as the determinant is a continuous function in the entries of the matrix, we have, in the large $N$ limit,

$$ \det \left[ \begin{array}{ccc} \frac{m_{\mathrm{SC}}(\tilde{\xi})}{\sigma_N} \mathbf{H}^{(\ell)} \mathbf{H}^{(\ell)\top} + \mathbf{I}_{r_\ell} & \frac{m_{\mathrm{SC}}(\tilde{\xi})}{\sigma_N} \mathbf{H}^{(\ell)} \mathbf{H}^{(\ell)\top} & \mathbf{0}_{r_\ell \times r_\ell} \\ \mathbf{0}_{r_\ell \times r_\ell} & \mathbf{I}_{r_\ell} & \mathbf{0}_{r_\ell \times r_\ell} \\ m_{\mathrm{SC}}(\tilde{\xi}) \mathbf{I}_{r_\ell} & m_{\mathrm{SC}}(\tilde{\xi}) \mathbf{I}_{r_\ell} & \mathbf{I}_{r_\ell} \end{array} \right] = 0. $$

Using twice the relation $\det[\frac{A}{C} \frac{B}{D}] = \det(AD^{-1} - BD^{-1}CD)$ when $D$ is invertible, this simplifies into

$$ \det \left( \frac{m_{\mathrm{SC}}(\tilde{\xi})}{\sigma_N} \mathbf{H}^{(\ell)} \mathbf{H}^{(\ell)\top} + \mathbf{I}_{r_\ell} \right) = 0 $$

Thus, we seek $\tilde{\xi}_{q_\ell}^{(\ell)} \in \mathbb{R} \setminus [-2,2]$ such that

$$ \frac{m_{\mathrm{SC}}(\tilde{\xi}_{q_\ell}^{(\ell)})}{\sigma_N} s_{q_\ell}^2(\mathbf{P}^{(\ell)}) + 1 = 0, \quad q_\ell \in [r_\ell]. $$

Injecting the expression $m_{\mathrm{SC}}(\tilde{\xi}_{q_\ell}^{(\ell)}) = -\frac{\sigma_N}{s_{q_\ell}^2(\mathbf{P}^{(\ell)})}$ into the equation $m_{\mathrm{SC}}^2(\tilde{\xi}_{q_\ell}^{(\ell)}) + \tilde{\xi}_{q_\ell}^{(\ell)} m_{\mathrm{SC}}(\tilde{\xi}_{q_\ell}^{(\ell)}) + 1 = 0$ yields

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

As $\tilde{\xi}_{q_\ell}^{(\ell)} > 0$ by definition, it must be strictly greater than 2 (the right edge of the semicircle). This is true only if $\rho_{q_\ell}^{(\ell)} \stackrel{\text{def}}{=} \frac{s_{q_\ell}^2(\mathbf{P}^{(\ell)})}{\sigma_N} > 1$.

### C.3 Eigenvector Alignments

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

$$ \tilde{\mathbf{Q}}^{(\ell)}(\tilde{z}) = \sum_{i_\ell=1}^{n_\ell} \frac{\hat{\mathbf{u}}_{q_\ell}^{(\ell)} \hat{\mathbf{u}}_{q_\ell}^{(\ell)\top}}{\frac{1}{\sigma_N} [s_{i_\ell}^2(\mathbf{T}^{(\ell)}) - \mu_N^{(\ell)}] - \tilde{z}}. $$

If $\rho_{q_\ell}^{(\ell)} > 1$, $q_\ell \in [r_\ell]$, then $s_{q_\ell}^2(\mathbf{T}^{(\ell)})$ is an isolated eigenvalue in the spectrum of $\mathbf{T}^{(\ell)}\mathbf{T}^{(\ell)\top}$ and $\frac{1}{\sigma_N} [s_{q_\ell}^2(\mathbf{T}^{(\ell)}) - \mu_N^{(\ell)}] \xrightarrow[N \to +\infty]{\text{a.s.}} \tilde{\xi}_{q_\ell}^{(\ell)}$. Hence, for any positively-oriented simple closed complex
contour $\gamma_{q_l}^{(\ell)}$ circling around $\tilde{\xi}_{q_l}^{(\ell)}$, leaving all the other $\tilde{\xi}_{q_l}^{(\ell)}$, $q_l' \neq q_l$, outside and not crossing $[-2, 2]$, Cauchy's integral formula yields, for $N$ large enough and any $\mathbf{a} \in \mathbb{R}^{n_l}$,

$$ \left| \mathbf{a}^{\top} \hat{\mathbf{u}}_{q_l}^{(\ell)} \right|^2 = -\frac{1}{2i\pi} \oint_{\gamma_{q_l}^{(\ell)}} \mathbf{a}^{\top} \tilde{\mathbf{Q}}^{(\ell)}(\tilde{z}) \mathbf{a} \, \mathrm{d}\tilde{z} \xrightarrow[N \to +\infty]{\text{a.s.}} -\frac{1}{2i\pi} \oint_{\gamma_{q_l}^{(\ell)}} \mathbf{a}^{\top} \bar{\mathbf{Q}}^{(\ell)}(\tilde{z}) \mathbf{a} \, \mathrm{d}\tilde{z} $$

by the dominated convergence theorem since, for all $\tilde{z} \in \gamma_{q_l}^{(\ell)}$, $\mathbf{a}^{\top} \tilde{\mathbf{Q}}^{(\ell)}(\tilde{z}) \mathbf{a} \to \mathbf{a}^{\top} \bar{\mathbf{Q}}^{(\ell)}(\tilde{z}) \mathbf{a}$ almost surely as $N \to +\infty$ by definition of the deterministic equivalent (Definition 3) and $\tilde{z} \in \gamma_{q_l}^{(\ell)} \mapsto |\mathbf{a}^{\top} \tilde{\mathbf{Q}}^{(\ell)}(\tilde{z}) \mathbf{a}|$ is almost surely bounded for $N$ large enough because we can choose $\gamma_{q_l}^{(\ell)}$ such that $\mathrm{Dist}(\tilde{\xi}_{q_l}^{(\ell)}, \gamma_{q_l}^{(\ell)}) \geq \varepsilon > 0$ and therefore $|\mathbf{a}^{\top} \tilde{\mathbf{Q}}^{(\ell)}(\tilde{z}) \mathbf{a}| \leq \|\mathbf{a}\|^2 \|\tilde{\mathbf{Q}}^{(\ell)}(\tilde{z})\| \leq \|\mathbf{a}\|^2/\varepsilon$ almost surely.

Using residue calculus, we can compute,

$$ -\frac{1}{2i\pi} \oint_{\gamma_{q_l}^{(\ell)}} \mathbf{a}^{\top} \bar{\mathbf{Q}}^{(\ell)}(\tilde{z}) \mathbf{a} \, \mathrm{d}\tilde{z} = - \lim_{\tilde{z} \to \tilde{\xi}_{q_l}^{(\ell)}} (\tilde{z} - \tilde{\xi}_{q_l}^{(\ell)}) \mathbf{a}^{\top} \left( \frac{1}{\sigma_N} \mathbf{P}^{(\ell)} \mathbf{P}^{(\ell)\top} + \frac{1}{m_{\mathrm{SC}}(\tilde{z})} \mathbf{I}_{n_l} \right)^{-1} \mathbf{a}. $$

Note that $\mathbf{P}^{(\ell)}\mathbf{P}^{(\ell)\top} = \mathbf{X}^{(\ell)}\mathbf{H}^{(\ell)}\mathbf{H}^{(\ell)\top}\mathbf{X}^{(\ell)\top}$ and there exist an $r_l \times r_l$ orthogonal matrix $\mathbf{O}^{(\ell)}$ such that $\mathbf{H}^{(\ell)}\mathbf{H}^{(\ell)\top} = \mathbf{O}^{(\ell)}\mathbf{\Lambda}^{(\ell)}\mathbf{O}^{(\ell)\top}$ with $\mathbf{\Lambda}^{(\ell)} = \mathrm{Diag}(s_1^2(\mathbf{P}^{(\ell)}), \dots, s_{r_l}^2(\mathbf{P}^{(\ell)}))$. Hence,

$$ \begin{aligned} & \mathbf{a}^{\top} \left( \frac{1}{\sigma_N} \mathbf{P}^{(\ell)} \mathbf{P}^{(\ell)\top} + \frac{1}{m_{\mathrm{SC}}(\tilde{z})} \mathbf{I}_{n_l} \right)^{-1} \mathbf{a} \\ & = \mathbf{a}^{\top} \mathbf{X}^{(\ell)} \mathbf{O}^{(\ell)} \left( \frac{1}{\sigma_N} \mathbf{\Lambda}^{(\ell)} + \frac{1}{m_{\mathrm{SC}}(\tilde{z})} \mathbf{I}_{n_l} \right)^{-1} \mathbf{O}^{(\ell)\top} \mathbf{X}^{(\ell)\top} \mathbf{a}. \end{aligned} $$

Let us therefore compute the following quantity,

$$ - \lim_{\tilde{z} \to \tilde{\xi}_{q_l}^{(\ell)}} (\tilde{z} - \tilde{\xi}_{q_l}^{(\ell)}) \left( \frac{s_{q_l'}^2(\mathbf{P}^{(\ell)})}{\sigma_N} + \frac{1}{m_{\mathrm{SC}}(\tilde{z})} \right)^{-1} = \begin{cases} 0 & \text{if } q_l' \neq q_l, \\ \zeta_{q_l}^{(\ell)} & \text{if } q_l' = q_l, \end{cases} \quad q_l' \in [r_l], $$

where we have used the fact that $m_{\mathrm{SC}}(\tilde{\xi}_{q_l}^{(\ell)}) = -\frac{\sigma_N}{s_{q_l}^2(\mathbf{P}^{(\ell)})}$. In order to handle the case $q_l' = q_l$, we use L'Hôpital's rule,

$$ \begin{aligned} \zeta_{q_l}^{(\ell)} &= - \left( \frac{\mathrm{d}}{\mathrm{d}\tilde{z}} \left[ \frac{s_{q_l}^2(\mathbf{P}^{(\ell)})}{\sigma_N} + \frac{1}{m_{\mathrm{SC}}(\tilde{z})} \right] \right)_{\tilde{z}=\tilde{\xi}_{q_l}^{(\ell)}}^{-1} \\ &= \frac{m_{\mathrm{SC}}^2(\tilde{\xi}_{q_l}^{(\ell)})}{m'_{\mathrm{SC}}(\tilde{\xi}_{q_l}^{(\ell)})} \\ &= \frac{\sigma_N^2}{s_{q_l}^4(\mathbf{P}^{(\ell)}) m'_{\mathrm{SC}}(\tilde{\xi}_{q_l}^{(\ell)})}. \end{aligned} $$

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

$$ m'_{\mathrm{SC}}(\tilde{z}) = - \frac{m_{\mathrm{SC}}(\tilde{z})}{2m_{\mathrm{SC}}(\tilde{z}) + \tilde{z}}. $$
Hence,

$$ 
\begin{aligned}
m'_{\text{SC}}(\tilde{\xi}_{q_{\ell}}^{(\ell)}) &= -\frac{-\frac{\sigma_N}{s_{q_{\ell}}^2(\mathbf{P}^{(\ell)})}}{-2\frac{\sigma_N}{s_{q_{\ell}}^2(\mathbf{P}^{(\ell)})} + \frac{s_{q_{\ell}}^2(\mathbf{P}^{(\ell)})}{\sigma_N} + \frac{\sigma_N}{s_{q_{\ell}}^2(\mathbf{P}^{(\ell)})}} \\
&= \frac{1}{\frac{s_{q_{\ell}}^4(\mathbf{P}^{(\ell)})}{\sigma_N^2} - 1}.
\end{aligned}
$$

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

$$ \zeta_{q_{\ell}}^{(\ell)} = 1 - \frac{\sigma_N^2}{s_{q_{\ell}}^4(\mathbf{P}^{(\ell)})}. $$

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

$$ -\frac{1}{2i\pi} \oint_{\gamma_{q_{\ell}}^{(\ell)}} \mathbf{a}^{\top} \bar{\mathbf{Q}}^{(\ell)}(\tilde{z}) \mathbf{a} \, d\tilde{z} = \mathbf{a}^{\top} \mathbf{X}^{(\ell)} \mathbf{O}^{(\ell)} \mathbf{Z}_{q_{\ell}}^{(\ell)} \mathbf{O}^{(\ell)\top} \mathbf{X}^{(\ell)\top} \mathbf{a} $$

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

$$ \| \mathbf{X}^{(\ell)\top} \hat{\mathbf{u}}_{q_{\ell}}^{(\ell)} \|^2 \xrightarrow[N \to +\infty]{\text{a.s.}} \zeta_{q_{\ell}}^{(\ell)} \sum_{q'_{\ell}=1}^{r_{\ell}} O_{q'_{\ell},q_{\ell}}^{(\ell)2} = \zeta_{q_{\ell}}^{(\ell)}. $$

## Appendix D. Proof of Lemma 14

Our proof of Lemma 14 uses the notion of $\varepsilon$-covering. An $\varepsilon$-covering of a compact set $\mathcal{K}$ for the norm $\|\cdot\|$ is a finite set $\mathcal{C} \subset \mathcal{K}$ such that for all $x \in \mathcal{K}$, there exists $\bar{x} \in \mathcal{C}$ such that $\|x - \bar{x}\| \le \varepsilon$. We also define the covering number $N(\varepsilon, \mathcal{K}, \|\cdot\|)$ 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}dt$ and the (upper) incomplete Gamma function $\Gamma(s, x) = \int_{x}^{+\infty} t^{s-1}e^{-t}dt$ for $s > 0$ and $x \ge 0$.

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

**Lemma 24** For $\ell \in [d]$ and $\varepsilon > 0$, let $\mathbf{\Delta}^{(\ell)} \in \mathbb{R}^{n_{\ell} \times r_{\ell}}$ be such that $\|\mathbf{\Delta}^{(\ell)}\| \le \varepsilon$ and $\mathbf{V}^{(\ell)} \in V_{r_{\ell}}(\mathbb{R}^{n_{\ell}})$ be the matrix of its left singular vectors. For all $\mathbf{A}^{(\ell')} \in \mathbb{R}^{n_{\ell'} \times r_{\ell'}}$, $\ell' \ne \ell$,

$$ \| \mathcal{N}(\mathbf{A}^{(1)}, \dots, \mathbf{\Delta}^{(\ell)}, \dots, \mathbf{A}^{(d)}) \|_{\mathrm{F}} \le \varepsilon \| \mathcal{N}(\mathbf{A}^{(1)}, \dots, \mathbf{V}^{(\ell)}, \dots, \mathbf{A}^{(d)}) \|_{\mathrm{F}}. $$

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

$$ 
\begin{aligned}
\| \mathcal{N}(\mathbf{A}^{(1)}, \dots, \mathbf{\Delta}^{(\ell)}, \dots, \mathbf{A}^{(d)}) \|_{\mathrm{F}} &= \left\| \mathbf{\Delta}^{(\ell)\top} \mathbf{N}^{(\ell)} \boxtimes_{\ell' \neq \ell} \mathbf{A}^{(\ell')} \right\|_{\mathrm{F}} \\
&= \left\| \mathbf{W}^{(\ell)\top} \mathbf{\Sigma}^{(\ell)} \mathbf{V}^{(\ell)\top} \mathbf{N}^{(\ell)} \boxtimes_{\ell' \neq \ell} \mathbf{A}^{(\ell')} \right\|_{\mathrm{F}} \\
&\le \underbrace{\| \mathbf{W}^{(\ell)\top} \mathbf{\Sigma}^{(\ell)} \|}_{=\varepsilon} \underbrace{\| \mathbf{V}^{(\ell)\top} \mathbf{N}^{(\ell)} \boxtimes_{\ell' \neq \ell} \mathbf{A}^{(\ell')} \|}_{=\varepsilon} \underbrace{\| \mathcal{N}(\mathbf{A}^{(1)}, \dots, \mathbf{V}^{(\ell)}, \dots, \mathbf{A}^{(d)}) \|_{\mathrm{F}}}_{=\varepsilon}.
\end{aligned}
$$
using the fact that $||\mathbf{A}\mathbf{B}||_F \le ||\mathbf{A}||||\mathbf{B}||_F$.

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

$$ \|\mathbf{N}(\mathbf{A}^{(1)}, \dots, \mathbf{A}^{(d)})\|_{F}^{2} \sim \chi^{2}\left(\prod_{\ell \in [d]} r_{\ell}\right). $$

**Proof** Firstly, observe that, for all $(q_1, \dots, q_d) \in \mathcal{X}_{\ell \in [d]}[r_\ell]$,

$$ [\mathcal{N}(\mathbf{A}^{(1)}, \dots, \mathbf{A}^{(d)})]_{q_1, \dots, q_d} = \sum_{i_1, \dots, i_d=1}^{n_1, \dots, n_d} \mathcal{N}_{i_1, \dots, i_d} A_{i_1, q_1}^{(1)} \dots A_{i_d, q_d}^{(d)} \sim \mathcal{N}(0, 1) $$

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

$$ \begin{aligned} & \mathbb{E}[[\mathcal{N}(\mathbf{A}^{(1)}, \dots, \mathbf{A}^{(d)})]_{q_1, \dots, q_d}[\mathcal{N}(\mathbf{A}^{(1)}, \dots, \mathbf{A}^{(d)})]_{q'_1, \dots, q'_d}] \\ &= \sum_{i_1, \dots, i_d=1}^{n_1, \dots, n_d} \sum_{i'_1, \dots, i'_d=1}^{n_1, \dots, n_d} \mathbb{E}[\mathcal{N}_{i_1, \dots, i_d} \mathcal{N}_{i'_1, \dots, i'_d}] A_{i_1, q_1}^{(1)} A_{i'_1, q'_1}^{(1)} \dots A_{i_d, q_d}^{(d)} A_{i'_d, q'_d}^{(d)} \\ &= \sum_{i_1, \dots, i_d=1}^{n_1, \dots, n_d} A_{i_1, q_1}^{(1)} A_{i'_1, q'_1}^{(1)} \dots A_{i_d, q_d}^{(d)} A_{i'_d, q'_d}^{(d)} \\ &= \begin{cases} 1 & \text{if } (q_1, \dots, q_d) = (q'_1, \dots, q'_d) \\ 0 & \text{otherwise} \end{cases}. \end{aligned} $$

Hence, the result follows from the fact that $||\mathcal{N}(\mathbf{A}^{(1)}, \dots, \mathbf{A}^{(d)})\|_F^2$ is the sum of $\prod_{\ell \in [d]} r_\ell$ squared independent $\mathcal{N}(0, 1)$ variables.

**Lemma 26 (Hinrichs et al., 2017, 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 $||\cdot||$,

$$ N(\varepsilon, V_r(\mathbb{R}^n), ||\cdot||) \le \left[ \frac{C}{\varepsilon} \right]^{r\left(n - \frac{r+1}{2}\right)}. $$

where $C > 0$ is a universal constant.

**Lemma 27** $\Gamma(s, x) \le \max(1, e^{s-1})\Gamma(s)e^{-x/2}$ for all $x \ge 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 \le 1$ when $C$ is well chosen. $f$ is continuously differentiable on $[0, +\infty[$ and

$$ f'(x) = \frac{1}{Ce^{-x/2}} \left( -x^{s-1}e^{-x} + \frac{1}{2}\Gamma(s, x) \right) \ge 0 \iff \Gamma(s, x) - 2x^{s-1}e^{-x} \ge 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'(x) = x^{s-1}e^{-x} - 2(s-1)x^{s-2}e^{-x} \ge 0 \iff x \ge 2(s-1).$$

We distinguish two cases.

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

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

<table><thead><tr><th>x</th><th>0</th><th>2(s - 1)</th><th>+∞</th></tr></thead><tbody><tr><td>g'(x)</td><td>-</td><td>0</td><td>+</td></tr><tr><td>g(x)</td><td colspan="3">Γ(s) <br> g(2(s - 1)) <br> 0</td></tr></tbody></table>

Since $g$ is strictly increasing on $[2(s - 1), +\infty[$ and $\lim_{x \to +\infty} g(x) = 0$, we necessarily have $g(2(s - 1)) < 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(s - 1)$. Let $x_0(s) \in [0, 2(s - 1)]$ 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)) \le \Gamma(s)$ and $e^{x_0(s)/2} \le e^{s-1}$.

<figure><img src="image_7.png" /></figure>

We are now ready to prove Lemma 14.

Let $\varepsilon > 0$ and $C_\ell$ be an $\varepsilon$-covering of $V_{r_\ell}(\mathbb{R}^{n_\ell})$ for the spectral norm $\|\cdot\|$, $\ell \in [d]$. Since $\times_{\ell \in [d]} V_{r_\ell}(\mathbb{R}^{n_\ell})$ is compact, it contains an element $(A_\star^{(1)}, \dots, A_\star^{(d)})$ such that

$$\sup_{A^{(\ell)} \in V_{r_\ell}(\mathbb{R}^{n_\ell}), \ell \in [d]} \left\| \mathcal{N}(A^{(1)}, \dots, A^{(d)}) \right\|_{\mathrm{F}} = \left\| \mathcal{N}(A_\star^{(1)}, \dots, A_\star^{(d)}) \right\|_{\mathrm{F}}.$$

Let $\bar{A}^{(\ell)} \in C_\ell$ be such that $A_\star^{(\ell)} = \bar{A}^{(\ell)} + \Delta^{(\ell)}$ with $\|\Delta^{(\ell)}\| \le \varepsilon$. Then, using the triangle inequality, Lemma 24 and the optimality of $(A_\star^{(1)}, \dots, A_\star^{(d)})$, we have,

$$\| \mathcal{N}(\bar{A}^{(1)} + \Delta^{(1)}, \dots, \bar{A}^{(d)} + \Delta^{(d)}) \|_{\mathrm{F}} \le \| \mathcal{N}(\bar{A}^{(1)}, \dots, \bar{A}^{(d)}) \|_{\mathrm{F}} + S \| \mathcal{N}(A_\star^{(1)}, \dots, A_\star^{(d)}) \|_{\mathrm{F}}$$

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

$$\| \mathcal{N}(A_\star^{(1)}, \dots, A_\star^{(d)}) \|_{\mathrm{F}} \le 2 \| \mathcal{N}(\bar{A}^{(1)}, \dots, \bar{A}^{(d)}) \|_{\mathrm{F}}$$
and, from the union bound, for any $t \ge 0$,

$$ \begin{aligned} \mathbb{P}(\|\mathcal{N}(\mathbf{A}^{(1)}, \dots, \mathbf{A}^{(d)})\|_{\mathrm{F}} \ge t) &\le \mathbb{P}\left(\bigcup_{\mathbf{A}^{(\ell)} \in \mathcal{C}_{\ell}, \ell \in [d]} \left\{\|\mathcal{N}(\mathbf{A}^{(1)}, \dots, \mathbf{A}^{(d)})\|_{\mathrm{F}} \ge \frac{t}{2}\right\}\right) \\ &\le \sum_{\mathbf{A}^{(\ell)} \in \mathcal{C}_{\ell}, \ell \in [d]} \mathbb{P}\left(\|\mathcal{N}(\mathbf{A}^{(1)}, \dots, \mathbf{A}^{(d)})\|_{\mathrm{F}} \ge \frac{t}{2}\right). \end{aligned} $$

Thus, combining Lemma 25 and 26, we have,

$$ \mathbb{P}(\|\mathcal{N}(\mathbf{A}^{(1)}, \dots, \mathbf{A}^{(d)})\|_{\mathrm{F}} \ge t) \le \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 \ge \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 27,

$$ \mathbb{P}\left(X \ge \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)} \le \max\left(1, e^{\frac{1}{2} \prod_{\ell \in [d]} r_{\ell} - 1}\right) e^{-t^2/16}. $$

We get the result stated in Lemma 14 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 E. Proof of Theorem 15

Recall the decomposition $\mathcal{P} = [\mathcal{H}; \mathbf{X}^{(1)}, \dots, \mathbf{X}^{(d)}]$. We use the following lemma whenever we state that $\|\mathcal{P}(\mathbf{U}^{(1)}, \dots, \mathbf{U}^{(d)})\|_{\mathrm{F}} = \mathcal{O}(\|\mathcal{P}\|_{\mathrm{F}})$ as $N \to +\infty$.

**Lemma 28** For all $(\mathbf{A}^{(1)}, \dots, \mathbf{A}^{(d)}) \in \mathcal{X}_{\ell \in [d]} V_{r_{\ell}}(\mathbb{R}^{n_{\ell}})$,

$$ \|\mathcal{P}(\mathbf{A}^{(1)}, \dots, \mathbf{A}^{(d)})\|_{\mathrm{F}} \le \|\mathcal{P}\|_{\mathrm{F}} \prod_{\ell \in [d]} \|\mathbf{X}^{(\ell)\top} \mathbf{A}^{(\ell)}\|. $$
**Proof** The proof relies on the property $||A B||_F \le ||A||_F ||B||_F$.

$$
\begin{align*}
\| \mathcal{P}(A^{(1)}, \dots, A^{(d)}) \|_{\mathrm{F}} &= \left\| A^{(1)\top} X^{(1)} H^{(1)} \boxtimes_{\ell=2}^{d} X^{(\ell)\top} A^{(\ell)} \right\|_{\mathrm{F}} \\
&\le \left\| A^{(1)\top} X^{(1)} \right\| \left\| H^{(1)} \boxtimes_{\ell=2}^{d} X^{(\ell)\top} A^{(\ell)} \right\|_{\mathrm{F}} \\
&= \left\| A^{(1)\top} X^{(1)} \right\| \left\| A^{(2)\top} X^{(2)} H^{(2)} \left( I_{r_1} \boxtimes \boxtimes_{\ell=3}^{d} X^{(\ell)\top} A^{(\ell)} \right) \right\|_{\mathrm{F}} \\
&\le \left\| A^{(1)\top} X^{(1)} \right\| \left\| A^{(2)\top} X^{(2)} \right\| \left\| H^{(2)} \left( I_{r_1} \boxtimes \boxtimes_{\ell=3}^{d} X^{(\ell)\top} A^{(\ell)} \right) \right\|_{\mathrm{F}} \\
&\dots \\
&\le \left( \prod_{\ell=1}^{d} \left\| A^{(\ell)\top} X^{(\ell)} \right\| \right) \underbrace{\left\| H^{(d)} \boxtimes_{\ell=1}^{d} I_{r_\ell} \right\|}_{= \|\mathcal{H}\|_{\mathrm{F}} = \|\mathcal{P}\|_{\mathrm{F}}} .
\end{align*}
$$

<figure><img src="image_8.png" /></figure>

Given $\ell \in [d]$, $U_1^{(\ell)}$ gathers the $r_\ell$ dominant left singular vectors of $T^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')}$, i.e.,
it is solution to

$$
\max_{U^{(\ell)} \in V_{r_\ell}(\mathbb{R}^{n_\ell})} \left\| U^{(\ell)\top} T^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_{\mathrm{F}}^2. \qquad (21)
$$

Consider also a solution $\tilde{U}_1^{(\ell)}$ to the following related problem

$$
\max_{U^{(\ell)} \in V_{r_\ell}(\mathbb{R}^{n_\ell})} \left\| U^{(\ell)\top} P^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_{\mathrm{F}}^2 . \qquad (22)
$$

Observe that, using the property $||A B||_F \le ||A||_F ||B||_F$, we have,

$$
\begin{align*}
\left\| U^{(\ell)\top} P^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_{\mathrm{F}}^2 &= \left\| U^{(\ell)\top} X^{(\ell)} H^{(\ell)} \boxtimes_{\ell' \neq \ell} X^{(\ell')\top} U_0^{(\ell')} \right\|_{\mathrm{F}}^2 \\
&\le \left\| U^{(\ell)\top} X^{(\ell)} \right\|^2 \left\| H^{(\ell)} \boxtimes_{\ell' \neq \ell} X^{(\ell')\top} U_0^{(\ell')} \right\|_{\mathrm{F}}^2 \le \left\| H^{(\ell)} \boxtimes_{\ell' \neq \ell} X^{(\ell')\top} U_0^{(\ell')} \right\|_{\mathrm{F}}^2
\end{align*}
$$

and this upper bound is only reached with $U^{(\ell)} = X^{(\ell)}O^{(\ell)}$, for any $r_\ell \times r_\ell$ orthogonal matrix $O^{(\ell)}$. Hence, $\tilde{U}_1^{(\ell)} = X^{(\ell)}O^{(\ell)}$. The strategy of our proof is to show that, as $N \to +\infty$, Problem (21) has the same solutions as Problem (22), which are known.
With the decomposition $\mathcal{T} = \mathcal{P} + \frac{1}{\sqrt{N}}\mathcal{N}$, we have,

$$ \left\| U^{(\ell)\top} T^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_{\mathrm{F}}^2 = \left\| U^{(\ell)\top} P^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_{\mathrm{F}}^2 + \frac{1}{N} \left\| U^{(\ell)\top} N^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_{\mathrm{F}}^2 \\ + \frac{2}{\sqrt{N}} \left\langle U^{(\ell)\top} P^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')}, U^{(\ell)\top} N^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\rangle_{\mathrm{F}}. $$

From Lemma 14, $\frac{1}{N}\|U^{(\ell)\top}N^{(\ell)}\boxtimes_{\ell' \neq \ell} U_0^{(\ell')}\|_{\mathrm{F}}^2 = \mathcal{O}(1)$ almost surely and

$$ \frac{2}{\sqrt{N}} \left| \left\langle U^{(\ell)\top} P^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')}, U^{(\ell)\top} N^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\rangle_{\mathrm{F}} \right| \\ \le \frac{2}{\sqrt{N}} \underbrace{\left\| U^{(\ell)\top} P^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_{\mathrm{F}}}_{=\mathcal{O}(\|\mathcal{P}\|_{\mathrm{F}})} \underbrace{\left\| U^{(\ell)\top} N^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_{\mathrm{F}}}_{=\mathcal{O}(\sqrt{N})}. $$

Therefore, for all $U^{(\ell)} \in V_{r_\ell}(\mathbb{R}^{n_\ell})$,

$$ \left\| U^{(\ell)\top} T^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_{\mathrm{F}}^2 = \left\| U^{(\ell)\top} P^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_{\mathrm{F}}^2 + \mathcal{O}(\|\mathcal{P}\|_{\mathrm{F}}) \quad \text{almost surely.} \quad (23) $$

In particular,

$$ \left\| U_1^{(\ell)\top} T^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_{\mathrm{F}}^2 = \left\| U_1^{(\ell)\top} P^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_{\mathrm{F}}^2 + \mathcal{O}(\|\mathcal{P}\|_{\mathrm{F}}) \quad \text{almost surely,} $$

$$ \left\| U_1^{(\ell)\top} T^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_{\mathrm{F}}^2 = \left\| \tilde{U}_1^{(\ell)\top} P^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_{\mathrm{F}}^2 + \mathcal{O}(\|\mathcal{P}\|_{\mathrm{F}}) \quad \text{almost surely,} $$

where the first equation is simply Equation (23) with $U^{(\ell)} = U_1^{(\ell)}$ and the second equation stems from the maximum over $U^{(\ell)} \in V_{r_\ell}(\mathbb{R}^{n_\ell})$ of both sides of Equation (23)$^{23}$. Hence,

$$ \left\| U_1^{(\ell)\top} P^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_{\mathrm{F}}^2 = \left\| \tilde{U}_1^{(\ell)\top} P^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_{\mathrm{F}}^2 + \mathcal{O}(\|\mathcal{P}\|_{\mathrm{F}}) \quad \text{almost surely.} \quad (24) $$

<footnote n="23">23. $\max\{\|U^{(\ell)\top}T^{(\ell)}\boxtimes_{\ell' \neq \ell} U_0^{(\ell')}\|_{\mathrm{F}}^2\} \le \max\{\|U^{(\ell)\top}P^{(\ell)}\boxtimes_{\ell' \neq \ell} U_0^{(\ell')}\|_{\mathrm{F}}^2\} + \max\{\mathcal{O}(\|\mathcal{P}\|_{\mathrm{F}})\}$ where each max is over $U^{(\ell)} \in V_{r_\ell}(\mathbb{R}^{n_\ell})$. Thus, $\|U_1^{(\ell)\top}T^{(\ell)}\boxtimes_{\ell' \neq \ell} U_0^{(\ell')}\|_{\mathrm{F}}^2 - \|\tilde{U}_1^{(\ell)\top}P^{(\ell)}\boxtimes_{\ell' \neq \ell} U_0^{(\ell')}\|_{\mathrm{F}}^2 = \mathcal{O}(\|\mathcal{P}\|_{\mathrm{F}})$ almost surely.</footnote>
Then, consider the singular value decomposition $X^{(\ell)\top} U_1^{(\ell)} = \sum_{q_\ell=1}^{r_\ell} s_{q_\ell}^{(\ell)} v_{q_\ell}^{(\ell)} w_{q_\ell}^{(\ell)\top}$.

$$
\begin{aligned}
& \left\| U_1^{(\ell)\top} P^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')} \right\|_F^2 = \left\| U_1^{(\ell)\top} X^{(\ell)} H^{(\ell)} \boxtimes_{\ell' \neq \ell} X^{(\ell')\top} U_0^{(\ell')} \right\|_F^2 \\
&= \left\| \sum_{q_\ell=1}^{r_\ell} s_{q_\ell}^{(\ell)} w_{q_\ell}^{(\ell)} v_{q_\ell}^{(\ell)\top} H^{(\ell)} \boxtimes_{\ell' \neq \ell} X^{(\ell')\top} U_0^{(\ell')} \right\|_F^2 \\
&= \sum_{q_\ell=1}^{r_\ell} s_{q_\ell}^{(\ell)2} \left\| v_{q_\ell}^{(\ell)\top} H^{(\ell)} \boxtimes_{\ell' \neq \ell} X^{(\ell')\top} U_0^{(\ell')} \right\|_F^2 \\
&= \sum_{q_\ell=1}^{r_\ell} \left\| v_{q_\ell}^{(\ell)\top} H^{(\ell)} \boxtimes_{\ell' \neq \ell} X^{(\ell')\top} U_0^{(\ell')} \right\|_F^2 - \sum_{q_\ell=1}^{r_\ell} (1 - s_{q_\ell}^{(\ell)2}) \left\| v_{q_\ell}^{(\ell)\top} H^{(\ell)} \boxtimes_{\ell' \neq \ell} X^{(\ell')\top} U_0^{(\ell')} \right\|_F^2 \\
&= \left\| H^{(\ell)} \boxtimes_{\ell' \neq \ell} X^{(\ell')\top} U_0^{(\ell')} \right\|_F^2 - \sum_{q_\ell=1}^{r_\ell} (1 - s_{q_\ell}^{(\ell)2}) \left\| v_{q_\ell}^{(\ell)\top} H^{(\ell)} \boxtimes_{\ell' \neq \ell} X^{(\ell')\top} U_0^{(\ell')} \right\|_F^2.
\end{aligned}
$$

Therefore, because $\|H^{(\ell)} \boxtimes_{\ell' \neq \ell} X^{(\ell)\top} U_0^{(\ell')}\|_F^2 = \|\tilde{U}_1^{(\ell)\top} P^{(\ell)} \boxtimes_{\ell' \neq \ell} U_0^{(\ell')}\|_F^2$, Equation (24) yields,

$$ \sum_{q_\ell=1}^{r_\ell} (1 - s_{q_\ell}^{(\ell)2}) \left\| v_{q_\ell}^{(\ell)\top} H^{(\ell)} \boxtimes_{\ell' \neq \ell} X^{(\ell')\top} U_0^{(\ell')} \right\|_F^2 = \mathcal{O}(\|\mathcal{P}\|_F) \quad \text{almost surely.} $$

Using the decomposition $v_{q_\ell}^{(\ell)} = \sum_{q_\ell=1}^{r_\ell} [v_{q_\ell}^{(\ell)}]_{q'_\ell} X^{(\ell)\top} x_{q'_\ell}^{(\ell)}$, we can see that

$$ \left\| v_{q_\ell}^{(\ell)\top} H^{(\ell)} \boxtimes_{\ell' \neq \ell} X^{(\ell')\top} U_0^{(\ell')} \right\|_F^2 = \sum_{q_\ell=1}^{r_\ell} [v_{q_\ell}^{(\ell)}]_{q'_\ell}^2 \left\| \mathcal{P}(U_0^{(1)}, \dots, x_{q_\ell}^{(\ell)}, \dots, U_0^{(d)}) \right\|_F^2 = \Theta(L_N^2). $$

Hence,

$$ \sum_{q_\ell=1}^{r_\ell} (1 - s_{q_\ell}^{(\ell)2}) = \mathcal{O} \left( \frac{\|\mathcal{P}\|_F}{L_N^2} \right) \quad \text{almost surely,} $$

which is the result stated in Theorem 15:

$$ \frac{1}{r_\ell} \| X^{(\ell)\top} U_1^{(\ell)} \|_F^2 = \frac{1}{r_\ell} \sum_{q_\ell=1}^{r_\ell} s_{q_\ell}^{(\ell)2} = 1 + \mathcal{O} \left( \frac{\|\mathcal{P}\|_F}{L_N^2} \right) \quad \text{almost surely.} $$

## References

Karim M. Abadir and Jan R. Magnus. *Matrix Algebra*. Econometric Exercises. Cambridge University Press, Cambridge, 2005. ISBN 978-0-521-53746-9. doi: 10.1017/CBO978051180800. URL https://www.cambridge.org/core/books/matrix-algebra/BCE8FD2D62006D4061F88E02615B5622.
Pierre-Antoine Absil, Robert Mahony, and Rodolphe Sepulchre. *Optimization Algorithms on Matrix Manifolds*. Princeton University Press, April 2009. ISBN 978-1-4008-3024-4. doi: 10.1515/9781400830244. URL https://www.degruyter.com/document/doi/10.1515/9781400830244/html. Publication Title: Optimization Algorithms on Matrix Manifolds.

Evrim Acar, Canan Aykut-Bingol, Haluk Bingol, Rasmus Bro, and Bülent Yener. Multiway analysis of epilepsy tensors. *Bioinformatics* (Oxford, England), 23(13):i10–18, July 2007. ISSN 1367-4811. doi: 10.1093/bioinformatics/btm210.

Animashree Anandkumar, Rong Ge, Daniel Hsu, and Sham Kakade. A Tensor Spectral Approach to Learning Mixed Membership Community Models. In *Proceedings of the 26th Annual Conference on Learning Theory*, pages 867–881. PMLR, June 2013. URL https://proceedings.mlr.press/v30/Anandkumar13.html. ISSN: 1938-7228.

Animashree Anandkumar, Rong Ge, Daniel Hsu, Sham M. Kakade, and Matus Telgarsky. Tensor Decompositions for Learning Latent Variable Models. *Journal of Machine Learning Research*, 15(80):2773–2832, 2014. ISSN 1533-7928. URL http://jmlr.org/papers/v15/anandkumar14b.html.

Zhidong Bai and Jack W. Silverstein. *Spectral analysis of large dimensional random matrices*, volume 20. Springer, 2010.

Zhidong Bai and Yanqing Yin. Convergence to the Semicircle Law. *The Annals of Probability*, 16(2):863–875, April 1988. ISSN 0091-1798, 2168-894X. doi: 10.1214/aop/1176991792. URL https://projecteuclid.org/journals/annals-of-probability/volume-16/issue-2/Convergence-to-Publisher: Institute of Mathematical Statistics.

Jinho Baik, Gérard Ben Arous, and Sandrine Péché. Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices. *The Annals of Probability*, 33(5):1643–1697, September 2005. ISSN 0091-1798, 2168-894X. doi: 10.1214/009117905000000233. URL https://projecteuclid.org/journals/annals-of-probability/volume-33/issue-5/Phase-transition-Publisher: Institute of Mathematical Statistics.

Afonso S. Bandeira, Amelia Perry, and Alexander S. Wein. Notes on computational-to-statistical gaps: predictions using statistical physics. *Portugaliae Mathematica*, 75 (2):159–186, December 2018. ISSN 0032-5155. doi: 10.4171/pm/2014. URL https://ems.press/journals/pm/articles/15908.

Marwa Banna, Florence Merlevède, and Magda Peligrad. On the limiting spectral distribution for a large class of symmetric random matrices with correlated entries. *Stochastic Processes and their Applications*, 125(7):2700–2726, July 2015. ISSN 0304-4149. doi: 10.1016/j.spa.2015.01.010. URL https://www.sciencedirect.com/science/article/pii/S0304414915000290.
Gérard Ben Arous, Song Mei, Andrea Montanari, and Mihai Nica. The Landscape of the Spiked Tensor Model. *Communications on Pure and Applied Mathematics*, 72(11):2282-2330, 2019. ISSN 1097-0312. doi: 10.1002/cpa.21861. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/cpa.21861.

Gérard Ben Arous, Reza Gheissari, and Aukosh Jagannath. Algorithmic thresholds for tensor PCA. *The Annals of Probability*, 48(4):2052-2087, July 2020. ISSN 0091-1798, 2168-894X. doi: 10.1214/19-AOP1415. URL https://projecteuclid.org/journals/annals-of-probability/volume-48/issue-4/Algorithmic-thresholds-for-tensor-pca. Publisher: Institute of Mathematical Statistics.

Gérard Ben Arous, Daniel Zhengyu Huang, and Jiaoyang Huang. Long random matrices and tensor unfolding. *The Annals of Applied Probability*, 33(6B):5753-5780, December 2023. ISSN 1050-5164, 2168-8737. doi: 10.1214/23-AAP1958. URL https://projecteuclid.org/journals/annals-of-applied-probability/volume-33/issue-6B/Long-rank-random-matrices-and-tensor-unfolding. Publisher: Institute of Mathematical Statistics.

Florent Benaych-Georges and Raj Rao Nadakuditi. The eigenvalues and eigenvectors of finite, low rank perturbations of large random matrices. *Advances in Mathematics*, 227(1):494-521, 2011. ISSN 0001-8708. doi: 10.1016/j.aim.2011.02.007. URL https://www.sciencedirect.com/science/article/pii/S0001870811000570.

Florent Benaych-Georges and Raj Rao Nadakuditi. The singular values and vectors of low rank perturbations of large rectangular random matrices. *Journal of Multivariate Analysis*, 111:120-135, October 2012. ISSN 0047-259X. doi: 10.1016/j.jmva.2012.04.019. URL https://www.sciencedirect.com/science/article/pii/S0047259X12001108.

Xuan Bi, Xiwei Tang, Yubai Yuan, Yanqing Zhang, and Annie Qu. Tensors in Statistics. *Annual Review of Statistics and Its Application*, 8(Volume 8, 2021):345-368, March 2021. ISSN 2326-8298, 2326-831X. doi: 10.1146/annurev-statistics-042720-020816. URL https://www.annualreviews.org/content/journals/10.1146/annurev-statistics-042720-020816. Publisher: Annual Reviews.

Patrick Billingsley. *Probability and Measure*. Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., 2012.

Åke Björck and Gene H. Golub. Numerical methods for computing angles between linear subspaces. *Mathematics of Computation*, 27(123):579-594, 1973. ISSN 0025-5718, 1088-6842. doi: 10.1090/S0025-5718-1973-0348991-3. URL https://www.ams.org/mcom/1973-27-123/S0025-5718-1973-0348991-3/.

Rasmus Bro and Claus A. Andersson. Improving the speed of multiway algorithms: Part II: Compression. *Chemometrics and Intelligent Laboratory Systems*, 42(1):105-113, August 1998. ISSN 0169-7439. doi: 10.1016/S0169-7439(98)00011-2. URL https://www.sciencedirect.com/science/article/pii/S0169743998000112.

Mireille Capitaine, Catherine Donati-Martin, and Delphine Féral. The largest eigenvalues of finite rank deformation of large Wigner matrices: convergence and nonuniversality of
the fluctuations. *The Annals of Probability*, 37(1), January 2009. ISSN 0091-1798. doi: 10.1214/08-AOP394. URL http://arxiv.org/abs/0706.0136. arXiv:0706.0136 [math].

Hongbin Chen, Jean-Christophe Mourrat, and Jiaming Xia. Statistical inference of finite-rank tensors. *Annales Henri Lebesgue*, 5:1161–1189, 2022. ISSN 2644-9463. doi: 10.5802/ahl.146. URL https://ahl.centre-mersenne.org/articles/10.5802/ahl.146/.

Wei-Kuo Chen, Madeline Handschy, and Gilad Lerman. Phase transition in random tensors with multiple independent spikes. *The Annals of Applied Probability*, 31(4): 1868–1913, August 2021. ISSN 1050-5164, 2168-8737. doi: 10.1214/20-AAP1636. URL https://projecteuclid.org/journals/annals-of-applied-probability/volume-31/issue-4/Phase-Publisher: Institute of Mathematical Statistics.

Antoine Chevreuil and Philippe Loubaton. On The Non-Detectability of Spiked Large Random Tensors. In *2018 IEEE Statistical Signal Processing Workshop (SSP)*, pages 443–447, June 2018. doi: 10.1109/SSP.2018.8450752. URL https://ieeexplore.ieee.org/document/8450752.

Yasuko Chikuse. *Statistics on Special Manifolds*, volume 174 of *Lecture Notes in Statistics*. Springer, New York, NY, 2003. ISBN 978-0-387-00160-9 978-0-387-21540-2. doi: 10.1007/978-0-387-21540-2. URL http://link.springer.com/10.1007/978-0-387-21540-2.

Andrzej Cichocki, Danilo P. Mandic, Anh Huy Phan, Cesar F. Caiafa, Guoxu Zhou, Qibin Zhao, and Lieven De Lathauwer. Tensor Decompositions for Signal Processing Applications: From two-way to multiway component analysis. *IEEE Signal Processing Magazine*, 32(2):145–163, March 2015. ISSN 1558-0792. doi: 10.1109/MSP.2013.2297439. URL https://ieeexplore.ieee.org/document/7038247. Conference Name: IEEE Signal Processing Magazine.

Pierre Comon. Tensors versus Matrices, usefulness and unexpected properties. In IEEE, editor, *IEEE Workshop on Statistical Signal Processing*, pages 780–788, Cardiff, United Kingdom, September 2009. IEEE. URL https://hal.archives-ouvertes.fr/hal-00417258.

Pierre Comon. Tensors: a Brief Introduction. *IEEE Signal Processing Magazine*, 31(3):44–53, May 2014. doi: 10.1109/MSP.2014.2298533. URL https://hal.archives-ouvertes.fr/hal-00923279. Publisher: Institute of Electrical and Electronics Engineers.

Romain Couillet and Zhenyu Liao. *Random Matrix Methods for Machine Learning*. Cambridge University Press, Cambridge, 2022. ISBN 978-1-00-912323-5. doi: 10.1017/9781009128490. URL https://www.cambridge.org/core/books/random-matrix-methods-for-machine-learning/6B681EB69

Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle. On the Best Rank-1 and Rank-(R1, R2, ..., RN) Approximation of Higher-Order Tensors. *SIAM Journal on Matrix Analysis and Applications*, 21(4):1324–1342, January 2000a. ISSN 0895-4798. doi: 10.1137/S0895479898346995. URL https://epubs.siam.org/doi/10.1137/S0895479898346995. Publisher: Society for Industrial and Applied Mathematics.
Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle. A Multilinear Singular Value Decomposition. *SIAM Journal on Matrix Analysis and Applications*, 21(4):1253-1278, January 2000b. ISSN 0895-4798. doi: 10.1137/S0895479896305696. URL https://epubs.siam.org/doi/10.1137/S0895479896305696. Publisher: Society for Industrial and Applied Mathematics.

Carl Eckart and Gale Young. The approximation of one matrix by another of lower rank. *Psychometrika*, 1(3):211–218, September 1936. ISSN 1860-0980. doi: 10.1007/BF02288367. URL https://doi.org/10.1007/BF02288367.

Samuel F. Edwards and Raymund C. Jones. The eigenvalue spectrum of a large symmetric random matrix. *Journal of Physics A: Mathematical and General*, 9(10):1595, October 1976. ISSN 0305-4470. doi: 10.1088/0305-4470/9/10/011. URL https://dx.doi.org/10.1088/0305-4470/9/10/011.

Hadi Fanaee-T and João Gama. EigenEvent: An Algorithm for Event Detection from Complex Data Streams in Syndromic Surveillance. *Intelligent Data Analysis*, 19(3):597–616, June 2015. ISSN 1088467X, 15714128. doi: 10.3233/IDA-150734. URL http://arxiv.org/abs/1406.3496. arXiv:1406.3496 [cs, stat].

Michael J. Feldman. Spiked singular values and vectors under extreme aspect ratios. *Journal of Multivariate Analysis*, 196:105187, July 2023. ISSN 0047-259X. doi: 10.1016/j.jmva.2023.105187. URL https://www.sciencedirect.com/science/article/pii/S0047259X23000337.

Michael J. Feldman and David L. Donoho. Sharp Recovery Thresholds of Tensor PCA Spectral Algorithms. *Advances in Neural Information Processing Systems*, 36:56628–56640, December 2023. URL https://proceedings.neurips.cc/paper_files/paper/2023/hash/b14d76c7266be21b338527cd25deac

Evgeny Frolov and Ivan Oseledets. Tensor methods and recommender systems. *WIREs Data Mining and Knowledge Discovery*, 7(3):e1201, 2017. ISSN 1942-4795. doi: 10.1002/widm.1201. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/widm.1201. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/widm.1201.

Delphine Féral and Sandrine Péché. The Largest Eigenvalue of Rank One Deformation of Large Wigner Matrices. *Communications in Mathematical Physics*, 272(1):185–228, May 2007. ISSN 1432-0916. doi: 10.1007/s00220-007-0209-3. URL https://doi.org/10.1007/s00220-007-0209-3.

Zoltán Füredi and János Komlós. The eigenvalues of random symmetric matrices. *Combinatorica*, 1(3):233–241, September 1981. ISSN 1439-6912. doi: 10.1007/BF02579329. URL https://doi.org/10.1007/BF02579329.

Jeffrey S. Geronimo and Theodore P. Hill. Necessary and sufficient condition that the limit of Stieltjes transforms is a Stieltjes transform. *Journal of Approximation Theory*, 121(1):54–60, March 2003. ISSN 0021-9045. doi: 10.1016/S0021-9045(02)00042-4. URL https://www.sciencedirect.com/science/article/pii/S0021904502000424.
José Henrique de M. Goulart, Romain Couillet, and Pierre Comon. A Random Matrix Perspective on Random Tensors. *Journal of Machine Learning Research*, 23(264):1–36, 2022. ISSN 1533-7928. URL http://jmlr.org/papers/v23/21-1038.html.

Lars Grasedyck, Daniel Kressner, and Christine Tobler. A literature survey of low-rank tensor approximation techniques. *GAMM-Mitteilungen*, 36(1):53–78, 2013. ISSN 1522-2608. doi: 10.1002/gamm.201310004. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/gamm.201310004. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/gamm.201310004.

Razvan Gurau. Universality for random tensors. *Annales de l'Institut Henri Poincaré, Probabilités et Statistiques*, 50(4):1474–1525, November 2014. ISSN 0246-0203. doi: 10.1214/13-AIHP567. URL https://projecteuclid.org/journals/annales-de-l-institut-henri-poincare-probabilites-et-sts. Publisher: Institut Henri Poincaré.

Wolfgang Hackbusch. *Tensor Spaces and Numerical Tensor Calculus*. Springer Series in Computational Mathematics. Springer, 2012. URL https://link.springer.com/book/10.1007/978-3-642-28027-6.

Christopher J. Hillar and Lek-Heng Lim. Most Tensor Problems Are NP-Hard. *Journal of the ACM*, 60(6):45:1–45:39, November 2013. ISSN 0004-5411. doi: 10.1145/2512329. URL https://doi.org/10.1145/2512329.

Aicke Hinrichs, Joscha Prochno, and Jan Vybíral. Entropy numbers of embeddings of Schatten classes. *Journal of Functional Analysis*, 273(10):3241–3261, November 2017. ISSN 0022-1236. doi: 10.1016/j.jfa.2017.08.008. URL https://www.sciencedirect.com/science/article/pii/S0022123617303221.

Frank L. Hitchcock. The Expression of a Tensor or a Polyadic as a Sum of Products. *Journal of Mathematics and Physics*, 6(1-4):164–189, 1927. ISSN 1467-9590. doi: 10.1002/sapm192761164. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/sapm192761164. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/sapm192761164.

Samuel B. Hopkins, Jonathan Shi, and David Steurer. Tensor principal component analysis via sum-of-square proofs. In *Proceedings of the 28th Conference on Learning Theory*, pages 956–1006. PMLR, June 2015. URL https://proceedings.mlr.press/v40/Hopkins15.html. ISSN: 1938-7228.

Samuel B. Hopkins, Tselil Schramm, Jonathan Shi, and David Steurer. Fast spectral algorithms from sum-of-squares proofs: tensor decomposition and planted sparse vectors. In *Proceedings of the forty-eighth annual ACM symposium on Theory of Computing*, STOC '16, pages 178–191, New York, NY, USA, 2016. Association for Computing Machinery. ISBN 978-1-4503-4132-5. doi: 10.1145/2897518.2897529. URL https://doi.org/10.1145/2897518.2897529.