# When Random Tensors meet Random Matrices

Mohamed El Amine Seddik, Maxime Guillaud, Romain Couillet

► To cite this version:

Mohamed El Amine Seddik, Maxime Guillaud, Romain Couillet. When Random Tensors meet Random Matrices. The Annals of Applied Probability, 2024, 34 (1A), 10.1214/23-AAP1962. hal-04102861

HAL Id: hal-04102861

https://inria.hal.science/hal-04102861v1

Submitted on 6 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
WHEN RANDOM TENSORS MEET RANDOM MATRICES

BY MOHAMED EL AMINE SEEDIK¹, MAXIME GUILLAUD¹ AND ROMAIN COUILLET²

¹Mathematical and Algorithmic Sciences Laboratory, Huawei Technologies France,
mohamed.seddik@huawei.com; maxime.guillaud@huawei.com

²Université Grenoble Alpes, CNRS, GIPSA-lab, Grenoble INP,
romain.couillet@univ-grenoble-alpes.fr

Relying on random matrix theory (RMT), this paper studies asymmetric order-d spiked tensor models with Gaussian noise. Using the variational definition of the singular vectors and values of (Lim, 2005) [15], we show that the analysis of the considered model boils down to the analysis of an equivalent spiked symmetric block-wise random matrix, that is constructed from contractions of the studied tensor with the singular vectors associated to its best rank-1 approximation. Our approach allows the exact characterization of the almost sure asymptotic singular value and alignments of the corresponding singular vectors with the true spike components, when $\frac{n_i}{\sum_{j=1}^d n_j} \to c_i \in (0, 1)$ with $n_i$'s the tensor dimensions. In contrast to other works that rely mostly on tools from statistical physics to study random tensors, our results rely solely on classical RMT tools such as Stein's lemma. Finally, classical RMT results concerning spiked random matrices are recovered as a particular case.

CONTENTS

<table><tbody><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Random matrix theory tools</td><td>6</td></tr><tr><td>2.1</td><td>The resolvent matrix</td><td>7</td></tr><tr><td>2.2</td><td>The Stieltjes transform</td><td>7</td></tr><tr><td>2.3</td><td>Gaussian calculations</td><td>8</td></tr><tr><td>3</td><td>The asymmetric rank-one spiked tensor model of order 3</td><td>9</td></tr><tr><td>3.1</td><td>Tensor singular value and vectors</td><td>9</td></tr><tr><td>3.2</td><td>Associated random matrix model</td><td>9</td></tr><tr><td>3.3</td><td>Limiting spectral measure of block-wise 3-order tensor contractions</td><td>12</td></tr><tr><td>3.4</td><td>Concentration of the singular value and the alignments</td><td>14</td></tr><tr><td>3.5</td><td>Asymptotic singular value and alignments</td><td>15</td></tr><tr><td>3.6</td><td>Cubic 3-order tensors: case $c_1 = c_2 = c_3 = \frac{1}{3}$</td><td>16</td></tr><tr><td>4</td><td>Random tensors meet random matrices</td><td>17</td></tr><tr><td>4.1</td><td>Limiting spectral measure</td><td>18</td></tr><tr><td>4.2</td><td>Limiting singular value and alignments</td><td>19</td></tr><tr><td>5</td><td>Generalization to arbitrary d-order tensors</td><td>20</td></tr><tr><td>5.1</td><td>Associated random matrix ensemble</td><td>20</td></tr><tr><td>5.2</td><td>Limiting spectral measure of block-wise d-order tensor contractions</td><td>21</td></tr><tr><td>5.3</td><td>Asymptotic singular value and alignments of hyper-rectangular tensors</td><td>23</td></tr><tr><td>6</td><td>Generalization to rank r tensor with orthogonal components</td><td>24</td></tr><tr><td>7</td><td>Discussion</td><td>25</td></tr><tr><td></td><td>Acknowledgments</td><td>26</td></tr></tbody></table>

MSC2020 subject classifications: Primary 60B20; secondary 15B52.

Keywords and phrases: random matrix theory, random tensor theory, spiked models.
<table><tr><td>A</td><td>Simulations</td><td>26</td></tr><tr><td>A.1</td><td>Matrix case</td><td>26</td></tr><tr><td>A.2</td><td>Order 3 tensors</td><td>26</td></tr><tr><td>B</td><td>Proofs</td><td>28</td></tr><tr><td>B.1</td><td>Derivative of tensor singular value and vectors</td><td>28</td></tr><tr><td>B.2</td><td>Proof of Theorem 3</td><td>30</td></tr><tr><td>B.3</td><td>Proof of Corollary 1</td><td>33</td></tr><tr><td>B.4</td><td>Proof of Theorem 4</td><td>33</td></tr><tr><td>B.5</td><td>Proof of Theorem 5</td><td>35</td></tr><tr><td>B.6</td><td>Proof of Corollary 3</td><td>36</td></tr><tr><td>B.7</td><td>Proof of Corollary 4</td><td>37</td></tr><tr><td>B.8</td><td>Proof of Corollary 5</td><td>38</td></tr><tr><td>B.9</td><td>Proof of Theorem 6</td><td>39</td></tr><tr><td>B.10</td><td>Proof of Corollary 6</td><td>40</td></tr><tr><td>B.11</td><td>Existence of g(z)</td><td>40</td></tr><tr><td>B.12</td><td>Proof of Theorem 7</td><td>40</td></tr><tr><td>B.13</td><td>Proof of Theorem 8</td><td>42</td></tr><tr><td>B.14</td><td>Additional lemmas</td><td>43</td></tr><tr><td>References</td><td></td><td>45</td></tr></table>

Notations: [n] denotes the set {1, ..., n}. The set of rectangular matrices of size m × n is denoted $M_{m,n}$. The set of square matrices of size n is denoted $M_n$. The set of d-order tensors of size $n_1 \times \cdots \times n_d$ is denoted $T_{n_1, \dots, n_d}$. The set of hyper-cubic tensors of size n and order d is denoted $T_n^d$. The notation $X \sim T_{n_1, \dots, n_d}(\mathcal{N}(0, 1))$ means that $X$ is a random tensor with i.i.d. Gaussian $\mathcal{N}(0, 1)$ entries. Scalars are denoted by lowercase letters as $a, b, c$. Vectors are denoted by bold lowercase letters as $a, b, c$. $e_i^d$ denotes the canonical vector in $\mathbb{R}^d$ with $[e_i^d]_j = \delta_{ij}$. Matrices are denoted by bold uppercase letters as $A, B, C$. Tensors are denoted as $A, B, C$. $T_{i_1, \dots, i_d}$ denotes the entry $(i_1, \dots, i_d)$ of the tensor $T$. $\langle u, v \rangle = \sum_i u_i v_i$ denotes the scalar product between $u$ and $v$, the $\ell_2$-norm of a vector $u$ is denoted as $\|u\|^2 = \langle u, u \rangle$. $\| \cdot \|$ denotes the spectral norm for tensors. $T(u^{(1)}, \dots, u^{(d)}) = \sum_{i_1, \dots, i_d} u_{i_1}^{(1)} \cdots u_{i_d}^{(d)} T_{i_1, \dots, i_d}$ denotes the contraction of tensor $T$ on the vectors given as arguments. Given some vectors $u^{(1)}, \dots, u^{(k)}$ with $k < d$, the contraction $T(u^{(1)}, \dots, u^{(k)}, \dots, \dots, \dots)$ denotes the resulting $(d-k)$-th order tensor. $\mathbb{S}^{N-1}$ denotes the $N$-dimensional unit sphere.

**1. Introduction.** The extraction of latent and low-dimensional structures from raw data is a key step in various machine learning and signal processing applications. Our present interest is in those modern techniques which rely on the extraction of such structures from a low-rank random tensor model [1] and which extends the ideas from matrix-type data to tensor structured data. We refer the reader to [23, 21, 25] and the references therein which introduce an extensive set of applications of tensor decomposition methods to machine learning, including dimensionality reduction, supervised and unsupervised learning, learning subspaces for feature extraction, low-rank tensor recovery etc. Although random matrix models have been extensively studied and well understood in the literature, the understanding of random tensor models is still in its infancy and the ideas from random matrix analysis do not easily extend to higher-order tensors. Indeed, the *resolvent* notion (see Definition 1) which is at the heart of random matrices does not generalize to tensors. In our present investigation, we consider the spiked tensor model, which consists in an observed d-order tensor $T \in T_{n_1, \dots, n_d}$ of the form

$$ (1) \qquad T = \beta \boldsymbol{x}^{(1)} \otimes \cdots \otimes \boldsymbol{x}^{(d)} + \frac{1}{\sqrt{N}} \mathbf{X}, $$
where $(\mathbf{x}^{(1)}, \dots, \mathbf{x}^{(d)}) \in \mathbb{S}^{n_1-1} \times \dots \times \mathbb{S}^{n_d-1}$, $\mathbf{X} \sim \mathbb{T}_{n_1, \dots, n_d}(\mathcal{N}(0,1))$ and $N = \sum_{i=1}^d n_i$.
Note that the tensor noise is normalized by $\sqrt{\sum_{i=1}^d n_i}$ with $n_i$ the tensor dimensions, since
the spectral norm of $\mathbf{X}$ is of order $\sqrt{\sum_{i=1}^d n_i}$ from Lemma 4. One aims at retrieving the
rank-1 component (or spike) $\beta \mathbf{x}^{(1)} \otimes \dots \otimes \mathbf{x}^{(d)}$ from the noisy tensor $\mathbf{T}$, where $\beta$ can be
seen as controlling the signal-to-noise ratio (SNR). The identification of the dominant rank-1 component is an important special case of the low-rank tensor approximation problem, a
noisy version of the classical canonical polyadic decomposition (CPD) [10, 13]. Extensive
efforts have been made to study the performance of low rank tensor approximation methods
in the large dimensional regime – when the tensor dimensions $n_i \to \infty$, however considering
symmetric tensor models where $n_1 = \dots = n_d$, $\mathbf{x}^{(1)} = \dots = \mathbf{x}^{(d)}$, and assuming that the
noise is symmetric [17, 20, 14, 9, 12, 8].

In particular, in the matrix case (i.e., $d=2$), the above spiked tensor model becomes a
so-called spiked matrix model. It is well-known that in the large dimensional regime, there
exists an order one critical value $\beta_c$ of the SNR below which it is information-theoretically
impossible to detect/recover the spike, while above $\beta_c$, it is possible to detect the spike and
recover the corresponding components in (at least) polynomial time using singular value
decomposition (SVD). This phenomenon is sometimes known as the BBP (Baik, Ben Arous,
and Péché) phase transition [2, 5, 7, 19].

In the (symmetric) spiked tensor model for $d \ge 3$, there also exists an order one critical value<fnref n="1" /> $\beta_c(d)$ (in the high-dimensional asymptotic) below which it is information-theoretically impossible to detect/recover the spike, while above $\beta_c(d)$ recovery is theoretically possible with the maximum likelihood (ML) estimator. Computing the maximum likelihood in the matrix case corresponds to the computation of the largest singular vectors of the considered matrix which has a polynomial time complexity, while for $d \ge 3$, computing ML is NP-hard [17, 6]. As such, a more practical phase transition for tensors is to characterize the algorithmic critical value $\beta_a(d, n)$ (which might depend on the tensor dimension $n$) above which the recovery of the spike is possible in polynomial time. Richard and Montanari [17] first introduced the symmetric spiked tensor model (of the form $\mathbf{Y} = \mu\mathbf{x}^{\otimes d} + \mathbf{W} \in \mathbb{T}_n^d$ with symmetric $\mathbf{W}$) and also considered the related algorithmic aspects. In particular, they used heuristics to highlight that spike recovery is possible, with Approximate Message Passing (AMP) or the tensor power iteration method, in polynomial time<fnref n="2" /> provided $\mu \gtrsim n^{\frac{d-1}{2}}$. This phase transition was later proven rigorously for AMP by [14, 12] and recently for tensor power iteration by [11].

Richard and Montanari [17] further introduced a method for tensor decomposition based
on tensor unfolding, which consists in unfolding $\mathbf{Y}$ to an $n^q \times n^{d-q}$ matrix $\text{Mat}(\mathbf{Y}) =
\mu \mathbf{x} \mathbf{y}^\top + \mathbf{Z}$ for $q \in [d-1]$, to which a SVD is then performed. Setting $q=1$, they pre-
dicted that their proposed method recovers successively the spike if $\mu \gtrsim n^{\frac{d-2}{4}}$. In a very
recent paper by Ben Arous et al. [3], a study of spiked long rectangular random matrices<fnref n="3" />
has been proposed under fairly general (bounded fourth-order moment) noise distribution as-
sumptions. They particularly proved the existence of a critical SNR for which the extreme
singular value and singular vectors exhibit a BBP-type phase transition. They applied their
result for the asymmetric rank-one spiked model in Eq. (1) (with equal dimensions) using the
tensor unfolding method, and found the exact threshold obtained by [17], i.e., $\beta \gtrsim n^{\frac{d-2}{4}}$ for
tensor unfolding to succeed in signal recovery.

<footnote n="1">¹Depending on the tensor order *d*. We will sometimes omit the dependence on *d* if there is no ambiguity.</footnote>

<footnote n="2">²Using tensor power iteration or AMP with random initialization.</footnote>

<footnote n="3">³Number of rows *m* are allowed to grow polynomially in the number of columns *n*, i.e., $\frac{m}{n} = n^\alpha$.</footnote>
For the asymmetric spiked tensor model in Eq. (1), few results are available in the literature (to the best of our knowledge only [3] considered this setting by applying the tensor unfolding method proposed by [17]). This is precisely the model we consider in the present work and our more general result is derived as follows. Given the asymmetric model from Eq. (1), the maximum likelihood (ML) estimator of the best rank-one approximation of $\mathbf{T}$ is given by

$$ (2) \qquad (\lambda_*, \mathbf{u}_*^{(i)}) = \operatorname*{arg\,min}_{\lambda \in \mathbb{R}^+, (\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)}) \in \mathbb{S}^{n_1-1} \times \dots \times \mathbb{S}^{n_d-1}} \|\mathbf{T} - \lambda \mathbf{u}^{(1)} \otimes \dots \otimes \mathbf{u}^{(d)}\|_{\mathrm{F}}^2. $$

In the above equation, $\lambda_*$ and the $\mathbf{u}_*^{(i)}$ can be respectively interpreted as the generalization to the tensor case of the concepts of dominant singular value and associated singular vectors [15]. Following variational arguments therein, Eq. (2) can be reformulated using contractions of $\mathbf{T}$ as

$$ (3) \qquad \max_{\prod_{i=1}^d \|\mathbf{u}^{(i)}\|=1} |\mathbf{T}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)})|, $$

the Lagrangian of which writes as $\mathcal{L} \equiv \mathbf{T}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)}) - \lambda (\prod_{i=1}^d \|\mathbf{u}^{(i)}\| - 1)$ with $\lambda > 0$. Hence, the stationary points $(\lambda, \mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)})$, with the $\mathbf{u}^{(i)}$'s being unitary vectors, must satisfy the Karush-Kuhn-Tucker conditions, for $i \in [d]$

$$ (4) \qquad \begin{cases} \mathbf{T}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(i-1)}, \mathbf{u}^{(i+1)}, \dots, \mathbf{u}^{(d)}) = \lambda \mathbf{u}^{(i)}, \\ \lambda = \mathbf{T}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)}). \end{cases} $$

An interesting question concerns the computation of the expected number of stationary points (local optima or saddle points) satisfying the identities in Eq. (4). [4] studied the landscape of a symmetric spiked tensor model and found that for $\beta < \beta_c$ the values of the objective function of all local maxima (including the global one) tend to concentrate on a small interval, while for $\beta > \beta_c$ the value achieved by the global maximum exits that interval and increases with $\beta$. In contrast, very recently Goulart et al. [8] have studied an order 3 symmetric spiked random tensor $\mathbf{Y}$ using a RMT approach, where it was stated that there exists a threshold $0 < \beta_s < \beta_c$ such that for $\beta \in [\beta_s, \beta_c]$ there exists a local optimum of the ML problem that correlates with the spike and such local optimum coincides with the global one for $\beta > \beta_c$. We conjecture that such observations extend to asymmetric spiked tensors, namely that there exists an order one critical value $\beta_c$ above which the ML problem in Eq. (3) admits a global maximum. As for [8], our present findings do not allow to express such $\beta_c$ and its exact characterization is left for future investigation. However, for asymmetric spiked random tensors, we also exhibit a threshold $\beta_s$ such that for $\beta > \beta_s$ there exists a local optimum of the ML objective that correlates with the true spike. Figure 1 provides an illustration of the different thresholds of $\beta$ and see the last part of Subsection 3.2 for a more detailed discussion.

**Main Contributions.** Starting from the conditions in Eq. (4), we provide an exact expression of the asymptotic singular value and alignments $\langle \mathbf{x}^{(i)}, \mathbf{u}_*^{(i)} \rangle$, when the tensor dimensions $n_i \to \infty$ with $\sum_{j=1}^d n_j \to c_i \in (0, 1)$, where the tuple $(\lambda_*, \mathbf{u}_*^{(1)}, \dots, \mathbf{u}_*^{(d)})$ is associated to a local optimum of the ML problem verifying some technical conditions (detailed in Assumption 4). We conjecture that when the SNR $\beta$ is large enough, there is a unique local optimum verifying Assumption 4 and for which our results characterize the corresponding alignments. We further conjecture that $(\lambda_*, \mathbf{u}_*^{(1)}, \dots, \mathbf{u}_*^{(d)})$ coincides with the global maximum above some$^4$ $\beta_c$ – that needs to be characterized.

<footnote>$^4$Such a critical value has been characterized by [12] for symmetric tensors. See [8] for a detailed discussion about this aspect in the case of symmetric tensors.</footnote>
<figure><img src="image_2.png" /><figcaption>FIG 1. Illustration of the different thresholds for the SNR $\beta$. Our approach exhibits $\beta_s$ such that for $\beta > \beta_s$ there exists a local optimum that correlates with the spike, the threshold $\beta_c$ is unknown for asymmetric tensors and corresponds to the ML phase transition (above $\beta_c$ the global maximum correlates with the spike), while $\beta_a$ corresponds to the algorithmic phase transition (recovery of the spike in polynomial time).</figcaption></figure>

Technically, we first show that the considered random tensor $\mathbf{T}$ can be mapped to an equiv- alent symmetric random matrix $\mathbf{T} \in \mathbb{M}_N$, constructed through contractions of $\mathbf{T}$ with $d-2$ directions among $\mathbf{u}_*^{(1)}, \dots, \mathbf{u}_*^{(d)}$. Then, leveraging on random matrix theory, we first char- acterize the limiting spectral measure of $\mathbf{T}$ and then provide estimates of the asymptotic alignments $\langle \mathbf{x}^{(i)}, \mathbf{u}_*^{(i)} \rangle$. We precisely show (see Theorem 8) that under Assumption 4, for $d \ge 3$, there exists $\beta_s > 0$ such that for $\beta > \beta_s$

$$ (5) \qquad \left\{ \begin{array}{l} \lambda_* \xrightarrow{\text{a.s.}} \lambda^\infty(\beta), \\ \left| \langle \mathbf{x}^{(i)}, \mathbf{u}_*^{(i)} \rangle \right| \xrightarrow{\text{a.s.}} q_i(\lambda^\infty(\beta)), \end{array} \right. $$

where $\lambda^\infty(\beta)$ satisfies $f(\lambda^\infty(\beta), \beta) = 0$ with $f(z, \beta) = z + g(z) - \beta \prod_{i=1}^d q_i(z)$ and

$$ q_i(z) = \sqrt{1 - \frac{g_i^2(z)}{c_i}}, \quad g_i(z) = \frac{g(z) + z}{2} - \frac{\sqrt{4c_i + (g(z) + z)^2}}{2}, $$

$g(z)$ being the solution to the fixed point equation $g(z) = \sum_{i=1}^d g_i(z)$ (for $z$ large enough); see Section B.11 for the existence of $g(z)$. Besides, for $\beta \in [0, \beta_s]$ with $^5 c_i = \frac{1}{d}$ for all $i \in [d]$

$$ (6) \qquad \lambda_* \xrightarrow{\text{a.s.}} \lambda^\infty \le 2\sqrt{\frac{d-1}{d}}, \quad \left| \langle \mathbf{x}^{(i)}, \mathbf{u}_*^{(i)} \rangle \right| \xrightarrow{\text{a.s.}} 0, $$

that is $\mathbf{u}_*^{(i)}$ (asymptotically) ceases to correlate with $\mathbf{x}^{(i)}$.

REMARK 1. Note that $q_i(z)$ can be equivalently expressed as $q_i(z) = \left( \frac{\alpha_i(z)^{d-3}}{\prod_{j \neq i} \alpha_j(z)} \right)^{\frac{1}{2d-4}}$ with $\alpha_i(z) = \frac{\beta}{z+g(z)-g_i(z)}$. Note that such expression is defined for $c_i \in [0, 1]$ with $d \ge 3$. See details in Section B.13.

We highlight that in our formulation the threshold $\beta_s$ corresponds to the minimal value of the SNR $\beta$ above which the derived asymptotic formulas are algebraically defined, which

<footnote>$^5$For arbitrary values of $c_i$'s, the upper bound of $\lambda^\infty$ can be computed numerically as the minimum non-negative real number $z$ for which Algorithm 1 converges.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG 2. Asymptotic singular value and alignments as per Eq. (5) when all the dimensions $n_i$ are equal ($c_i = \frac{1}{d}$ for all $i \in [d]$) for different values of the tensor order $d$. The matrix case corresponds to setting $d=3$ and $c_3 = 0$ (see Section 4 for more details).</figcaption></figure>

may differ from the true phase transition $\beta_c$ of the ML problem: In the case of symmetric tensors, the results form [8] seem to indicate that $\beta_s$ is slightly below $\beta_c$ obtained by [12] where the ML problem was studied.

Figure 2 notably depicts the asymptotic alignments as in Eq. (5) when all the tensor dimensions $n_i$ are equal. Since our result characterizes the alignments for $d \ge 3$, it is not possible to recover the matrix case by simply setting $d=2$ since the equations are not defined for $d=2$ (see Remark 1). However, the matrix case is recovered by considering an order 3 tensor ($d=3$) and then taking the limit $c_3 \to 0$ (equivalent to the degenerate case $n_3 = 1$ which results in a spiked random matrix model), see Section 4 for more details. From Figure 2-(b), unlike the matrix case, i.e., $d=2$, the predicted asymptotic alignments are not continuous for orders $d \ge 3$, this phenomenon has already been observed in the case of symmetric random tensors [12]. In particular, the predicted theoretical threshold $\beta_s$ in the matrix case $d=2$ coincides with the classical so-called BBP (Baik, Ben Arous, and Péché) phase transition $\beta_c(2)$ [2, 5, 7, 19]. Moreover, our result for the matrix case characterizes the asymptotic alignments for the long rectangular matrices studied by [3] and we also recover the threshold $\beta \gtrsim n^{\frac{d-2}{4}}$ for the case of tensor unfolding method (See remark 9). From a methodological view point, our results are derived based solely on Stein's lemma without the use of complex contour integrals as classically performed in RMT. In essence, we follow the approach in [8] and further introduce a new object (the mapping $\Phi_d$ defined subsequently in Eq. (31)) that simplifies drastically the use of RMT tools.

The remainder of the paper is organized as follows. Section 2 recalls some fundamental random matrix theory tools. In Section 3, we study asymmetric spiked tensors of order 3 where we provide the main steps of our approach. Section 4 characterizes the behavior of spiked random matrices given our result on spiked 3-order tensor models from Section 3. The generalization of our results to arbitrary $d$-order tensors is then presented in Section 5. Section 6 discusses the application of our findings to arbitrary rank-$k$ tensors with mutually orthogonal components. Further discussions are presented in Section 7. In Appendix A, we provide some simulations to support our findings and discuss some algorithmic aspects. Finally, Appendix B provides the proofs of the main developed results.

**2. Random matrix theory tools.** Before digging into our main findings, we briefly recall some random matrix theory results that are at the heart of our analysis. Specifically, we recall the *resolvent* formalism which allows to assess the spectral behavior of large *symmetric* random matrices. In particular, given a symmetric matrix $S \in \mathbb{M}_n$ and denoting
$\lambda_1(S) \le \dots \le \lambda_n(S)$ its $n$ eigenvalues with corresponding eigenvectors $u_i(S)$ for $i \in [n]$, its spectral decomposition writes as

$$S = \sum_{i=1}^{n} \lambda_{i}(S) u_{i}(S) u_{i}(S)^{\top}.$$

In the sequel we omit the dependence on $S$ by simply writing $\lambda_i$ and $u_i$ if there is no ambiguity.

2.1. The resolvent matrix. We start by defining the resolvent of a symmetric matrix and present its main properties.

DEFINITION 1. Given a symmetric matrix $S \in \mathbb{M}_n$, the resolvent of $S$ is defined as

$$R_S(z) = (S - zI_n)^{-1}, \quad z \in \mathbb{C} \setminus S(S),$$

where $S(S) = \{\lambda_1, \dots, \lambda_n\}$ is the spectrum of $S$.

The resolvent is a fundamental object since it retrieves the spectral characteristics (specrum and eigenvectors) of $S$. It particularly verifies the following property which we will use extensively to derive our main results,

$$(7) \qquad R_S(z) = -\frac{1}{z}I_n + \frac{1}{z}SR_S(z) = -\frac{1}{z}I_n + \frac{1}{z}R_S(z)S.$$

The above identity, coupled with Stein's Lemma (Lemma 1 in Section 2.3 below), is a fundamental tool used to derive fixed point equations that allow the evaluation of functionals of interests involving $R_S(z)$. Another interesting property of the resolvent concerns its spectral norm, which we denote by $\|R_S(z)\|$. Indeed, if the spectrum of $S$ has a bounded support, then the spectral norm of $R_S(z)$ is bounded. This is a consequence of the inequality

$$(8) \qquad \|R_S(z)\| \le \frac{1}{\text{dist}(z, S(S))}$$

where $\text{dist}(z, \cdot)$ denotes the distance of $z$ to a set. The resolvent encodes rich information about the behavior of the eigenvalues of $S$ through the so-called Stieltjes transform which we describe subsequently.

2.2. The Stieltjes transform. Random matrix theory, originally, aims at describing the limiting spectral measure of random matrices when their dimensions grow large. Typically, under certain technical conditions on $S$, the empirical spectral measure of $S$ defined as

$$(9) \qquad \nu_S = \frac{1}{n} \sum_{i=1}^{n} \delta_{\lambda_i(S)},$$

where $\delta_x$ is a Dirac mass on $x$, converges to a deterministic probability measure $\nu$. To characterize such asymptotic measure, the Stieltjes transform (defined below) approach is a widely used tool.

DEFINITION 2 (Stieltjes transform). Given a probability measure $\nu$, the Stieltjes transform of $\nu$ is defined by

$$g_{\nu}(z) \equiv \int \frac{d\nu(\lambda)}{\lambda - z}, \quad z \in \mathbb{C} \setminus S(\nu).$$
Particularly, the Stieltjes transform of $\nu_S$ is closely related to its associated resolvent $R_S(z)$ through the algebraic identity

$$g_{\nu_S}(z) = \frac{1}{n} \sum_{i=1}^{n} \int \frac{\delta_{\lambda_i(S)}(d\lambda)}{\lambda - z} = \frac{1}{n} \sum_{i=1}^{n} \frac{1}{\lambda_i(S) - z} = \frac{1}{n} \operatorname{tr} R_S(z)$$

The Stieltjes transform $g_\nu$ has several interesting analytical properties, among which, we have

1. $g_\nu$ is complex analytic on its definition domain $\mathbb{C} \setminus S(\nu)$ and $\mathfrak{F}[g_\nu(z)] > 0$ if $\mathfrak{F}[z] > 0$;

2. $g_\nu(z)$ is bounded for $z \in \mathbb{C} \setminus S(\nu)$ if $S(\nu)$ is bounded, since $|g_\nu(z)| \le \operatorname{dist}(z, S(\nu))^{-1}$;

3. Since $g'_\nu(z) = \int (\lambda - z)^{-2} d\nu(\lambda) > 0$, $g_\nu$ is monotonously increasing for $z \in \mathbb{R}$.

The Stieltjes transform $g_\nu$ admits an inverse formula which provides access to the evaluation of the underlying probability measure $\nu$, as per the following theorem.

**THEOREM 1 (Inverse formula of the Stieltjes transform).** Let $a, b$ be some continuity points of the probability measure $\nu$, then the segment $[a, b]$ is measurable with $\nu$, precisely

$$\nu([a, b]) = \frac{1}{\pi} \lim_{\epsilon \to 0} \int_a^b \mathfrak{F}[g_\nu(x + i\epsilon)] dx.$$

Moreover, if $\nu$ admits a density function $f$ at some point $x$, i.e., $\nu(x)$ is differentiable in a neighborhood of $x$ with $\lim_{\epsilon \to 0} \epsilon^{-1}\nu([x - \epsilon/2, x + \epsilon/2]) = f(x)$, then we have the inverse formula

$$f(x) = \frac{1}{\pi} \lim_{\epsilon \to 0} \mathfrak{F}[g_\nu(x + i\epsilon)].$$

When it comes to large random matrices, the Stieltjes transform admits a continuity property, in the sense that if a sequence of random probability measures converges to a deterministic measure then the corresponding Stieltjes transforms converge almost surely to a deterministic one, and vice-versa. The following theorem from [24] states precisely this continuity property.

**THEOREM 2 (Stieltjes transform continuity).** A sequence of random probability measures $\nu_n$, supported on $\mathbb{R}$ with corresponding Stieltjes transforms $g_{\nu_n}(z)$, converges almost surely weakly to a deterministic measure $\nu$, with corresponding Stieltjes transform $g_\nu$ if and only if $g_{\nu_n}(z) \to g_\nu(z)$ almost surely, for all $z \in \mathbb{R} + i\mathbb{R}_+$.

In particular, Theorems 1 and 2 combined together allow the description of the spectrum of large symmetric random matrices.

**2.3. Gaussian calculations.** The following lemma by Stein (also called Stein's identity or Gaussian integration by parts) allows to replace the expectation of the product of a Gaussian variable with a differentiable function $f$ by the variance of that variable times the expectation of $f'$.

**LEMMA 1 (Stein's Lemma [22]).** Let $X \sim \mathcal{N}(0, \sigma^2)$ and $f$ a continuously differentiable function having at most polynomial growth, then

$$\mathbb{E}[X f(X)] = \sigma^2 \mathbb{E}[f'(X)],$$

when the above expectations exist.
We will further need the Poincaré inequality which allows to control the variance of a functional of i.i.d. Gaussian random variables.

LEMMA 2 (Poincaré inequality). Let $F: \mathbb{R}^n \to \mathbb{R}$ a continuously differentiable function having at most polynomial growth and $X_1, \dots, X_n$ a collection of i.i.d. standard Gaussian variables. Then,

$$ \mathrm{Var}F(X_1, \dots, X_n) \leq \sum_{i=1}^{n} \mathbb{E} \left| \frac{\partial F}{\partial X_i} \right|^2 . $$

**3. The asymmetric rank-one spiked tensor model of order 3.** To best illustrate our approach, we start by considering the asymmetric rank-one tensor model of order 3 of the form

$$ (10) \qquad \mathbf{T} = \beta \mathbf{x} \otimes \mathbf{y} \otimes \mathbf{z} + \frac{1}{\sqrt{N}} \mathbf{X}, $$

where $\mathbf{T}$ and $\mathbf{X}$ both have dimensions $m \times n \times p$ and $N = m + n + p$. We assume that $\mathbf{x}, \mathbf{y}$ and $\mathbf{z}$ are on the unit spheres $\mathbb{S}^{m-1}, \mathbb{S}^{n-1}$ and $\mathbb{S}^{p-1}$ respectively, $\mathbf{X}$ is a Gaussian noise tensor with i.i.d. entries $X_{ijk} \sim \mathcal{N}(0, 1)$ independent from $\mathbf{x}, \mathbf{y}, \mathbf{z}$.

**3.1. Tensor singular value and vectors.** According to Eq. (4), the $l_2$-singular value and vectors, corresponding to the best rank-one approximation $\lambda_* \mathbf{u}_* \otimes \mathbf{v}_* \otimes \mathbf{w}_*$ of $\mathbf{T}$, satisfy the identities

$$ (11) \qquad \mathbf{T}(\mathbf{v})\mathbf{w} = \lambda\mathbf{u} \quad \mathbf{T}(\mathbf{u})\mathbf{w} = \lambda\mathbf{v} \quad \mathbf{T}(\mathbf{v})^{\top}\mathbf{u} = \lambda\mathbf{w} \quad \text{with} \quad (\mathbf{u}, \mathbf{v}, \mathbf{w}) \in \mathbb{S}^{m-1} \times \mathbb{S}^{n-1} \times \mathbb{S}^{p-1}, $$

where we denoted $\mathbf{T}(\mathbf{u}) = \mathbf{T}(\mathbf{u}, \cdot, \cdot) \in \mathbb{M}_{m,p}$, $\mathbf{T}(\mathbf{v}) = \mathbf{T}(\cdot, \mathbf{v}, \cdot) \in \mathbb{M}_{m,p}$ and $\mathbf{T}(\mathbf{w}) = \mathbf{T}(\cdot, \cdot, \mathbf{w}) \in \mathbb{M}_{m,n}$. Furthermore, the singular value $\lambda$ can be characterized through the contraction of $\mathbf{T}$ along all its singular vectors $\mathbf{u}, \mathbf{v}, \mathbf{w}$, i.e.

$$ (12) \qquad \lambda = \mathbf{T}(\mathbf{u}, \mathbf{v}, \mathbf{w}) = \sum_{ijk} u_i v_j w_k T_{ijk}. $$

**3.2. Associated random matrix model.** We follow the approach developed in [8] which consists in studying random matrices that are obtained through contractions of a random tensor model, and extend it to the asymmetric spiked tensor model of Eq. (10). Indeed, it has been shown in [8] that the study of a rank-one symmetric spiked tensor model $\mathbf{Y} \in \mathbb{T}_q^3$ boils down to the analysis of the symmetric random matrix$^6 \mathbf{Y}(\mathbf{a}) \in \mathbb{M}_q$ where $\mathbf{a} \in \mathbb{S}^{q-1}$ stands for the eigenvector of $\mathbf{Y}$ corresponding to the best symmetric rank-one approximation of $\mathbf{Y}$, i.e.,

$$ (13) \qquad (\mu, \mathbf{a}) = \arg \min_{\mu \in \mathbb{R}, \mathbf{a} \in \mathbb{S}^{q-1}} \|\mathbf{Y} - \mu \mathbf{a}^{\otimes 3}\|_{F}^{2} \quad \Leftrightarrow \quad \mathbf{Y}(\mathbf{a})\mathbf{a} = \mu\mathbf{a}. $$

In the asymmetric case of Eq. (10), the choice of a “relevant” random matrix to study is not trivial since the corresponding contractions $\mathbf{T}(\mathbf{u}), \mathbf{T}(\mathbf{v})$ and $\mathbf{T}(\mathbf{w})$ yield asymmetric random matrices which present more technical difficulties from the random matrix theory perspective, therefore the extension of the approach in [8] to asymmetric tensors is not straightforward. We will see in the following that such a choice of the relevant matrix to study an asymmetric $\mathbf{T}$ is naturally obtained through the use of the Pastur’s Stein approach [16].

<footnote>$^6$The contraction of the tensor $\mathbf{Y}$ with its eigenvector $\mathbf{a}$.</footnote>
As described in [8] and since the singular vectors $u, v$ and $w$ depend statistically on $X$, the first technical challenge consists in expressing the derivatives of the singular vectors $u, v$ and $w$ w.r.t. the entries of the Gaussian noise tensor $X$. Indeed, one can show that there exists a differentiable mapping $\mathcal{F}: \mathbb{T}_{m,n,p} \to \mathbb{R}^{m+n+p+1}$ that maps $X$ to $\mathcal{F}(X) = (\lambda(X), u(X), v(X), w(X))$ singular-value and vectors of $T$, since the components of $u(X), v(X)$ and $w(X)$ are bounded and $\lambda(X)$ has polynomial growth. Indeed, we have the following Lemma which is analog to [8, Lemma 8] and which justifies the application of Stein's Lemma subsequently.

**LEMMA 3.** There exists an almost everywhere continuously differentiable function $\mathcal{F}: \mathbb{T}_{m,n,p} \to \mathbb{R}^{m+n+p+1}$ such that $\mathcal{F}(X) = (\lambda(X), u(X), v(X), w(X))$ is singular-value and vectors of $T$ (for almost every $X$).

**PROOF.** The proof relies on the same arguments as [8, Lemma 8]. □

Calculus (see details in B.1) show that deriving the identities in Eq. (11) w.r.t. an entry $X_{ijk}$ of $X$ with $(i, j, k) \in [m] \times [n] \times [p]$ results in

$$ (14) \quad \left[ \frac{\partial u}{\partial X_{ijk}} \quad \frac{\partial v}{\partial X_{ijk}} \quad \frac{\partial w}{\partial X_{ijk}} \right] = - \frac{1}{\sqrt{N}} \left( \left[ \begin{array}{cc} 0_{m \times m} & \mathbf{T}(w) \mathbf{T}(v) \\ \mathbf{T}(w)^{\top} & 0_{n \times n} \end{array} \right] \mathbf{T}(u) - \lambda \mathbf{I}_N \right)^{-1} \left[ \begin{array}{c} v_j w_k (e_i^m - u_i u) \\ u_i w_k (e_j^n - v_j v) \\ u_i v_j (e_k^p - w_k w) \end{array} \right] \in \mathbb{R}^N, $$

where we recall that $e_i^d$ denotes the canonical vector in $\mathbb{R}^d$ with $[e_i^d]_j = \delta_{ij}$. We further have the identity

$$ (15) \qquad \frac{\partial \lambda}{\partial X_{ijk}} = \frac{1}{\sqrt{N}} u_i v_j w_k. $$

Denoting by $M$ the symmetric block-wise random matrix which appears in the matrix inverse in Eq. (14), the derivatives of $u, v$ and $w$ are therefore expressed in terms of the resolvent $R_M(z) = (M - z\mathbf{I}_N)^{-1}$ evaluated on $\lambda$, and we will see subsequently that the assessment of the spectral properties of $T$ boils down to the estimation of the normalized trace $\frac{1}{N} \text{tr} R_M(z)$. As such, the matrix $M$ provides the associated random matrix model that encodes the spectral properties of $T$. We will henceforth focus our analysis on this random matrix, in order to assess the spectral behavior of $T$. More generally, we will be interested in studying random matrices from the 3-order block-wise tensor contraction ensemble $\mathcal{B}_3(\mathbf{X})$ for $\mathbf{X} \sim \mathbb{T}_{m,n,p}(\mathcal{N}(0, 1))$ defined as

$$ (16) \qquad \mathcal{B}_3(\mathbf{X}) = \{ \Phi_3(\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c}) \mid (\mathbf{a}, \mathbf{b}, \mathbf{c}) \in \mathbb{S}^{m-1} \times \mathbb{S}^{n-1} \times \mathbb{S}^{p-1} \}, $$

where $\Phi_3$ is the mapping

$$ (17) \qquad \Phi_3: \mathbb{T}_{m,n,p} \times \mathbb{S}^{m-1} \times \mathbb{S}^{n-1} \times \mathbb{S}^{p-1} \to \mathbb{M}_{m+n+p} \\ \qquad \qquad \qquad (\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c}) \mapsto \left[ \begin{array}{cc} \mathbf{0}_{m \times m} & \mathbf{X}(\mathbf{c}) \mathbf{X}(\mathbf{b}) \\ \mathbf{X}(\mathbf{c})^{\top} & \mathbf{0}_{n \times n} \end{array} \right] \left[ \begin{array}{c} \mathbf{X}(\mathbf{b})^{\top} \mathbf{X}(\mathbf{a}) \\ \mathbf{X}(\mathbf{a})^{\top} \mathbf{0}_{p \times p} \end{array} \right]. $$

We associate tensor $T$ to the random matrix

$$ (18) \qquad \Phi_3(\mathbf{T}, \mathbf{u}, \mathbf{v}, \mathbf{w}) = \beta \mathbf{V} \mathbf{B} \mathbf{V}^{\top} + \frac{1}{\sqrt{N}} \Phi_3(\mathbf{X}, \mathbf{u}, \mathbf{v}, \mathbf{w}) \in \mathbb{M}_N, $$
where

$$ B = \begin{bmatrix} 0 & \langle z, w \rangle & \langle y, v \rangle \\ \langle z, w \rangle & 0 & \langle x, u \rangle \\ \langle y, v \rangle & \langle x, u \rangle & 0 \end{bmatrix} \in \mathbb{M}_3, \quad V = \begin{bmatrix} x & 0_m & 0_m \\ 0_n & y & 0_n \\ 0_p & 0_p & z \end{bmatrix} \in \mathbb{M}_{N,3}. $$

and $x, y$ and $z$ are on the unit spheres $S^{m-1}$, $S^{n-1}$ and $S^{p-1}$ respectively. Note that $\Phi_3(\mathbf{T}, u, v, w)$ is symmetric and behaves as a so-called spiked random matrix model where the signal part $\beta V B V^\top$ correlates with the true signals $x, y$ and $z$ for a sufficiently large $\beta$. However, the noise part (i.e., the term $\frac{1}{\sqrt{N}}\Phi_3(\mathbf{X}, u, v, w)$ in Eq. (18)) has a non-trivial structure due to the statistical dependencies between the singular vectors $u, v, w$ and the tensor noise $\mathbf{X}$. Despite this statistical dependency, we will show in the next subsection that the spectral measure of $\Phi_3(\mathbf{T}, u, v, w)$ (and that of $\frac{1}{\sqrt{N}}\Phi_3(\mathbf{X}, u, v, w)$) converges to a deterministic measure (see Theorem 4) that coincides with the limiting spectral measure of $\Phi_3(\mathbf{T}, a, b, c)$ where $a, b, c$ are unit vectors and independent of $\mathbf{X}$. Furthermore, the matrix $\Phi_3(\mathbf{T}, u, v, w)$ admits $2\lambda$ as an eigenvalue (regardless of the value of $\beta$), which is a simple consequence of the identities in Eq. (11), since

$$ (19) \qquad \Phi_3(\mathbf{T}, u, v, w) \begin{bmatrix} u \\ v \\ w \end{bmatrix} = \begin{bmatrix} \mathbf{T}(w)v + \mathbf{T}(v)w \\ \mathbf{T}(w)^\top u + \mathbf{T}(u)w \\ \mathbf{T}(v)^\top u + \mathbf{T}(u)^\top v \end{bmatrix} = 2\lambda \begin{bmatrix} u \\ v \\ w \end{bmatrix}. $$

Note however that the expression in Eq. (14) exists only if $\lambda$ is not an eigenvalue of $\Phi_3(\mathbf{T}, u, v, w)$. As discussed in [8] in the symmetric case, such condition is related to the locality of the maximum ML estimator. Indeed, we first have the following remark that concerns the Hessian of the underlying objective function.

**REMARK 2.** Recall the Lagrangian of the ML problem as $\mathcal{L}(u, v, w) = \mathbf{T}(u, v, w) - \lambda (\|u\|\|v\|\|w\| - 1)$ and denote $h = \frac{[u^\top, v^\top, w^\top]^\top}{\sqrt{3}}$. If $\lambda$ is a singular value of $\mathbf{T}$ with associated singular vectors $u, v, w$, then $(\lambda, h)$ is an eigenpair of the Hessian $\nabla^2 \mathcal{L} \equiv \frac{\partial^2 \mathcal{L}}{\partial h^2}$ evaluated at $h$. Indeed, $\nabla^2 \mathcal{L}(h) = \Phi_3(\mathbf{T}, u, v, w) - \lambda \mathbf{I}_N$, thus $\nabla^2 \mathcal{L}(h)h = \Phi_3(\mathbf{T}, u, v, w)h - \lambda h = \lambda h$ since $\Phi_3(\mathbf{T}, u, v, w)h = 2\lambda h$ (from Eq. (19)).

From [18, Theorem 12.5], a necessary condition for $u, v, w$ to be a local maximum of the ML problem is that

$$ \langle \nabla^2 \mathcal{L}(h) k, k \rangle \le 0, \quad \forall k \in h^\perp = \{k \in \mathbb{R}^N \mid \langle h, k \rangle = 0\}. $$

which yields (by remark 2) the condition

$$ (20) \qquad \max_{k \in S^{N-1} \cap h^\perp} \langle \Phi_3(\mathbf{T}, u, v, w)k, k \rangle \le \lambda. $$

As such, for $\lambda > 0$, $2\lambda$ must be the largest eigenvalue of $\Phi_3(\mathbf{T}, u, v, w)$ as per Eq. (19), while its second largest eigenvalue cannot exceed $\lambda$ as shown by Eq. (20). Thus our analysis applies only to some local optimum of the ML problem for which the corresponding singular value $\lambda$ lies outside the bulk of $\Phi_3(\mathbf{T}, u, v, w)$, namely we suppose$^7$ that there exists a tuple $(\lambda, u, v, w)$ verifying the identities in Eq. (11) such that $\lambda$ is not an eigenvalue of $\Phi_3(\mathbf{T}, u, v, w)$. Moreover, this allows the existence of the matrix inverse in Eq. (14) and it is not restrictive in the sense that it is satisfied for a sufficiently large value of the SNR $\beta$. In the sequel, for any $(\lambda, u, v, w)$ verifying the identities in Eq. (11), we find that the largest

<footnote>$^7$We will find that such assumption is satisfied provided $\beta > \beta_s$ for some $\beta_s > 0$ that we will determine</footnote>
<figure><img src="image_4.png" /><figcaption>FIG 3. Spectrum of $\frac{1}{\sqrt{N}}\Phi_3(\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c})$ with $\mathbf{X} \sim T_{m,n,p}(\mathcal{N}(0,1))$ and $\mathbf{a}, \mathbf{b}, \mathbf{c}$ independently sampled from the unit spheres $\mathbb{S}^{m-1}, \mathbb{S}^{n-1}, \mathbb{S}^{p-1}$ respectively. In black the semi-circle law as per Theorem 3.</figcaption></figure>

eigenvalue $2\lambda$ is always isolated from the bulk of $\Phi_3(\mathbf{T}, \mathbf{u}, \mathbf{v}, \mathbf{w})$ independently of the SNR $\beta$. In addition, there exists $\beta_s > 0$ such that for $\beta \le \beta_s$, the observed spike is not informative (in the sense that the corresponding alignments will be null) and is visible (an isolated eigenvalue appears in the spectrum of $\Phi_3(\mathbf{T}, \mathbf{u}, \mathbf{v}, \mathbf{w})$, see Figure 4) because of the statistical dependencies between $\mathbf{u}, \mathbf{v}, \mathbf{w}$ and $\mathbf{X}$. The same phenomenon has been observed in the case of symmetric random tensors by [8].

**3.3. Limiting spectral measure of block-wise 3-order tensor contractions.** We start by presenting our first result which characterizes the limiting spectral measure of the ensemble $\mathcal{B}_3(\mathbf{X})$ with $\mathbf{X} \sim T_{m,n,p}(\mathcal{N}(0,1))$. We characterize this measure in the limit when the tensor dimensions grow large as in the following assumption.

**ASSUMPTION 1 (Growth rate).** As $m, n, p \to \infty$, the dimension ratios $\frac{m}{m+n+p} \to c_1 \in (0,1)$, $\frac{n}{m+n+p} \to c_2 \in (0,1)$ and $\frac{p}{m+n+p} \to c_3 = 1 - (c_1+c_2) \in (0,1)$.

We have the following theorem which characterizes the spectrum of $\frac{1}{\sqrt{N}}\Phi_3(\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c})$ with arbitrary deterministic unit vectors $\mathbf{a}, \mathbf{b}, \mathbf{c}$.

**THEOREM 3.** Let $\mathbf{X} \sim T_{m,n,p}(\mathcal{N}(0,1))$ be a sequence of random asymmetric Gaussian tensors and $(\mathbf{a}, \mathbf{b}, \mathbf{c}) \in \mathbb{S}^{m-1} \times \mathbb{S}^{n-1} \times \mathbb{S}^{p-1}$ a sequence of deterministic vectors of increasing dimensions, following Assumption 1. Then the empirical spectral measure of $\frac{1}{\sqrt{N}}\Phi_3(\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c})$ converges weakly almost surely to a deterministic measure $\nu$ whose Stieltjes transform $g(z)$ is defined as the solution to the equation $g(z) = \sum_{i=1}^3 g_i(z)$ such that $\Im[g(z)] > 0$ for $\Im[z] > 0$ where, for $i \in [3]$, $g_i(z)$ satisfies $g_i^2(z) - (g(z)+z)g_i(z) - c_i = 0$ for $z \in \mathbb{C} \setminus \mathcal{S}(\nu)$.

**PROOF.** See Appendix B.2. $\square$

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

**REMARK 3.** The equation $g(z) = \sum_{i=1}^3 g_i(z)$ yields a polynomial equation of degree 5 in $g$ which can be solved numerically or via an iteration procedure which converges to a fixed point for $z$ large enough.

In particular, in the cubic case when $c_1 = c_2 = c_3 = \frac{1}{3}$, the empirical spectral measure of $\frac{1}{\sqrt{N}}\Phi_3(\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c})$ converges to a semi-circle law as precisely stated by the following Corollary of Theorem 3.
**COROLLARY 1.** Given the setting of Theorem 3 with $c_1 = c_2 = c_3 = \frac{1}{3}$, the empirical spectral measure of $\frac{1}{\sqrt{N}}\Phi_3(\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c})$ converges weakly almost surely to the semi-circle distribution supported on $S(\nu) = [-2\sqrt{\frac{2}{3}}, 2\sqrt{\frac{2}{3}}]$, whose density and Stieltjes transform write respectively as

$$ \nu(dx) = \frac{3}{4\pi} \sqrt{\left(\frac{8}{3} - x^2\right)^+} dx, \quad g(z) = \frac{-3z + 3\sqrt{z^2 - \frac{8}{3}}}{4}, \quad \text{where } z \in \mathbb{C} \setminus S(\nu). $$

**PROOF.** See Appendix B.3. □

Figure 3 depicts the spectrum of $\frac{1}{\sqrt{N}}\Phi_3(\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c})$ with $\mathbf{X} \sim T_{m,n,p}(\mathcal{N}(0, 1))$ and independent unit vectors $\mathbf{a}, \mathbf{b}, \mathbf{c}$, it particularly illustrates the convergence in law of this spectrum when the dimensions $m, n, p$ grow large.

**REMARK 4.** More generally, the spectral measure of $\frac{1}{\sqrt{N}}\Phi_3(\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c})$ converges to a deterministic measure with connected support if $\frac{1}{\sqrt{N}}\Phi_3(\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c})$ is almost surely full rank, i.e., if $\max(m, n) - \min(m, n) \le p \le m + n$ since the rank of $\frac{1}{\sqrt{N}}\Phi_3(\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c})$ is $\min(m, n + p) + \min(n, m + p) + \min(p, m + n)$. In contrast if it is not full rank, its spectral measure converges to a deterministic measure with unconnected support (see the case of matrices in Corollary 4 subsequently).

The analysis of the tensor $\mathbf{T}$ relies on describing the spectrum of $\Phi_3(\mathbf{T}, \mathbf{u}, \mathbf{v}, \mathbf{w})$ where the singular vectors $\mathbf{u}, \mathbf{v}, \mathbf{w}$ depend statistically on $\mathbf{X}$ (the noise part of $\mathbf{T}$). Despite these dependencies, it turns out that the spectrum of $\Phi_3(\mathbf{T}, \mathbf{u}, \mathbf{v}, \mathbf{w})$ converges in law to the same deterministic measure described by Theorem 3. Besides, we need a further technical assumption on the singular value and vectors of $\mathbf{T}$.

**ASSUMPTION 2.** We assume that there exists a sequence of critical points $(\lambda_*, \mathbf{u}_*, \mathbf{v}_*, \mathbf{w}_*)$ satisfying Eq. (11) such that $\lambda_* \xrightarrow{a.s.} \lambda^\infty(\beta)$, $|\langle \mathbf{u}_*, \mathbf{x} \rangle| \xrightarrow{a.s.} a_x^\infty(\beta)$, $|\langle \mathbf{v}_*, \mathbf{y} \rangle| \xrightarrow{a.s.} a_y^\infty(\beta)$, $|\langle \mathbf{w}_*, \mathbf{z} \rangle| \xrightarrow{a.s.} a_z^\infty(\beta)$ with $\lambda^\infty(\beta) \notin S(\nu)$ and $a_x^\infty(\beta), a_y^\infty(\beta), a_z^\infty(\beta) > 0$.

We precisely have the following result.

**THEOREM 4.** Let $\mathbf{T}$ be a sequence of random tensors defined as in Eq. (10). Under Assumptions 1 and 2, the empirical spectral measure of $\Phi_3(\mathbf{T}, \mathbf{u}_*, \mathbf{v}_*, \mathbf{w}_*)$ converges weakly almost surely to a deterministic measure $\nu$ whose Stieltjes transform $g(z)$ is defined as the solution to the equation $g(z) = \sum_{i=1}^3 g_i(z)$ such that $\Im[g(z)] > 0$ for $\Im[z] > 0$ where, for $i \in [3]$ $g_i(z)$ satisfies $g_i^2(z) - (g(z) + z)g_i(z) - c_i = 0$ for $z \in \mathbb{C} \setminus S(\nu)$.

**PROOF.** See Appendix B.4. □

However, the statistical dependency between $\mathbf{u}, \mathbf{v}, \mathbf{w}$ and $\mathbf{X}$ exhibits an isolated eigenvalue in the spectrum of $\Phi_3(\mathbf{T}, \mathbf{u}, \mathbf{v}, \mathbf{w})$ at the value $2\lambda$ independently of the value of the SNR $\beta$, which is a consequence of Eq. (19). Figure 4 depicts iterations of the power iteration method where the leftmost histogram corresponds to the fixed point solution, where one sees that the spike converges to the value $2\lambda$.
<figure><img src="image_6.png" /><figcaption>FIG 4. Spectrum of $\Phi_3(\mathbf{T}, \mathbf{u}, \mathbf{v}, \mathbf{w})$ at iterations 0, 5 and $\infty$ of the power method (see Algorithm 3) applied on $\mathbf{T}$, for $m=n=p=100$ and $\beta=0$.</figcaption></figure>

### REMARK 5.
Note that, from Assumption 2, the almost sure limit $\lambda^\infty$ of $\lambda$ must lie outside the support $S(\nu)$ of the deterministic measure $\nu$ described by Theorem 4. The same phenomenon was noticed in [8] where it is assumed that the almost sure limit $\mu^\infty$ of $\mu$ must satisfy $\mu^\infty > 2\sqrt{\frac{2}{3}}$ in the case of symmetric tensors.

Let us denote the blocks of the resolvent of $\Phi_3(\mathbf{T}, \mathbf{u}, \mathbf{v}, \mathbf{w})$ as

$$ (21) \qquad \mathbf{R}_{\Phi_3}(\mathbf{T}, \mathbf{u}, \mathbf{v}, \mathbf{w})(z) = \begin{bmatrix} \mathbf{R}^{11}(z) & \mathbf{R}^{12}(z) & \mathbf{R}^{13}(z) \\ \mathbf{R}^{12}(z)^\top & \mathbf{R}^{22}(z) & \mathbf{R}^{23}(z) \\ \mathbf{R}^{13}(z)^\top & \mathbf{R}^{23}(z)^\top & \mathbf{R}^{33}(z) \end{bmatrix}, $$

where $\mathbf{R}^{11}(z) \in \mathbb{M}_m$, $\mathbf{R}^{22}(z) \in \mathbb{M}_n$ and $\mathbf{R}^{33}(z) \in \mathbb{M}_p$. The following corollary of Theorem 4 will be useful subsequently.

### COROLLARY 2.
Recall the setting and notations of Theorem 4. For all $z \in \mathbb{C} \setminus S(\nu)$, we have

$$ \frac{1}{N} \operatorname{tr} \mathbf{R}^{11}(z) \xrightarrow{\text{a.s.}} g_1(z), \quad \frac{1}{N} \operatorname{tr} \mathbf{R}^{22}(z) \xrightarrow{\text{a.s.}} g_2(z), \quad \frac{1}{N} \operatorname{tr} \mathbf{R}^{33}(z) \xrightarrow{\text{a.s.}} g_3(z). $$

### 3.4. Concentration of the singular value and the alignments.
When the dimensions of $\mathbf{T}$ grow large under Assumption 1, the singular value $\lambda$ and the alignments $\langle \mathbf{u}, \mathbf{x} \rangle$, $\langle \mathbf{v}, \mathbf{y} \rangle$ and $\langle \mathbf{w}, \mathbf{z} \rangle$ converge almost surely to some deterministic limits. This can be shown by controlling the variances of these quantities using the Poincaré's inequality (Lemma 2). Precisely, for $\lambda$, invoking Eq. (15) we have

$$ \operatorname{Var} \lambda \le \sum_{ijk} \mathbb{E} \left| \frac{\partial \lambda}{\partial X_{ijk}} \right|^2 = \frac{1}{N} \sum_{ijk} u_i^2 v_j^2 w_k^2 = \frac{1}{N}. $$

Bounding higher order moments of $\lambda$ similarly allows to obtain the concentration of $\lambda$, e.g. by Chebyshev's inequality, we have for all $t > 0$

$$ \mathbb{P}(|\lambda - \mathbb{E}\lambda| \ge t) \le \frac{1}{Nt^2}. $$

Similarly, with Eq. (14), there exists $C > 0$ such that for all $t > 0$

$$ \mathbb{P}(|\langle \mathbf{u}, \mathbf{x} \rangle - \mathbb{E}\langle \mathbf{u}, \mathbf{x} \rangle| \ge t) \le \frac{C}{Nt^2}, \quad \mathbb{P}(|\langle \mathbf{v}, \mathbf{y} \rangle - \mathbb{E}\langle \mathbf{v}, \mathbf{y} \rangle| \ge t) \le \frac{C}{Nt^2}, \quad \mathbb{P}(|\langle \mathbf{w}, \mathbf{z} \rangle - \mathbb{E}\langle \mathbf{w}, \mathbf{z} \rangle| \ge t) \le \frac{C}{Nt^2}. $$

For the remainder of the manuscript, we denote the almost sure limits of $\lambda$ and the of alignments $\langle \mathbf{u}, \mathbf{x} \rangle$, $\langle \mathbf{v}, \mathbf{y} \rangle$ and $\langle \mathbf{w}, \mathbf{z} \rangle$ respectively as

$$ (22) \qquad \lambda^{\infty}(\beta) \equiv \lim_{N \to \infty} \lambda, \quad a_{\mathbf{x}}^{\infty}(\beta) \equiv \lim_{N \to \infty} \langle \mathbf{u}, \mathbf{x} \rangle, \quad a_{\mathbf{y}}^{\infty}(\beta) \equiv \lim_{N \to \infty} \langle \mathbf{v}, \mathbf{y} \rangle, \quad a_{\mathbf{z}}^{\infty}(\beta) \equiv \lim_{N \to \infty} \langle \mathbf{w}, \mathbf{z} \rangle. $$
3.5. Asymptotic singular value and alignments. Having the concentration result from the previous subsection, it remains to estimate the expectations $\mathbb{E}\lambda, \mathbb{E}\langle\mathbf{u}, \mathbf{x}\rangle, \mathbb{E}\langle\mathbf{v}, \mathbf{y}\rangle$ and $\mathbb{E}\langle\mathbf{w}, \mathbf{z}\rangle$. Usually, the evaluation of these quantities using tools from random matrix theory relies on computing Cauchy integrals involving $R_{\Phi_3(\mathbf{T},\mathbf{u},\mathbf{v},\mathbf{w})}(z)$ (the resolvent of $\Phi_3(\mathbf{T}, \mathbf{u}, \mathbf{v}, \mathbf{w})$). Here, we take a different (yet analytically simpler) approach by taking directly the expectation of the identities in Eq. (11) and Eq. (12), then applying Stein's Lemma (Lemma 1) and Lemma 3. For instance, for $\lambda$, we have

$$ \mathbb{E}\lambda = \frac{1}{\sqrt{N}} \sum_{ijk} \mathbb{E} \left[ v_j w_k \frac{\partial u_i}{\partial X_{ijk}} \right] + \mathbb{E} \left[ u_i w_k \frac{\partial v_j}{\partial X_{ijk}} \right] + \mathbb{E} \left[ u_i v_j \frac{\partial w_k}{\partial X_{ijk}} \right] + \beta \mathbb{E} \left[ \langle \mathbf{u}, \mathbf{x} \rangle \langle \mathbf{v}, \mathbf{y} \rangle \langle \mathbf{w}, \mathbf{z} \rangle \right]. $$

From Eq. (14), when $N \to \infty$, it turns out that the only contributing terms⁸ of the derivatives $\frac{\partial u_i}{\partial X_{ijk}}$, $\frac{\partial v_j}{\partial X_{ijk}}$ and $\frac{\partial w_k}{\partial X_{ijk}}$ in the above sum are respectively

$$ \frac{\partial u_i}{\partial X_{ijk}} \approx -\frac{1}{\sqrt{N}} v_j w_k R_{ii}^{11}(\lambda), \quad \frac{\partial v_j}{\partial X_{ijk}} \approx -\frac{1}{\sqrt{N}} u_i w_k R_{jj}^{22}(\lambda), \quad \frac{\partial w_k}{\partial X_{ijk}} \approx -\frac{1}{\sqrt{N}} u_i v_j R_{kk}^{33}(\lambda). $$

This yields

$$ \mathbb{E}\lambda = -\frac{1}{N} \left( \operatorname{tr} \mathbf{R}^{11}(\lambda) + \operatorname{tr} \mathbf{R}^{22}(\lambda) + \operatorname{tr} \mathbf{R}^{33}(\lambda) \right) + \beta \mathbb{E} \left[ \langle \mathbf{u}, \mathbf{x} \rangle \langle \mathbf{v}, \mathbf{y} \rangle \langle \mathbf{w}, \mathbf{z} \rangle \right] + \mathcal{O}(N^{-1}). $$

Therefore, the almost sure limit $\lambda^\infty(\beta)$ as $N \to \infty$ of $\lambda$ satisfies

$$ \lambda^\infty(\beta) + g(\lambda^\infty(\beta)) = \beta a_x^\infty(\beta) a_y^\infty(\beta) a_z^\infty(\beta), $$

where $a_x^\infty(\beta)$, $a_y^\infty(\beta)$ and $a_z^\infty(\beta)$ are defined in Eq. (22). From Eq. (11), proceeding similarly as above with the identities

$$ \mathbf{x}^{\top} \mathbf{T}(\mathbf{v}) \mathbf{w} = \lambda \langle \mathbf{u}, \mathbf{x} \rangle, \quad \mathbf{y}^{\top} \mathbf{T}(\mathbf{u}) \mathbf{w} = \lambda \langle \mathbf{v}, \mathbf{y} \rangle, \quad \mathbf{z}^{\top} \mathbf{T}(\mathbf{v})^{\top} \mathbf{u} = \lambda \langle \mathbf{w}, \mathbf{z} \rangle, $$

we obtain the following result.

**THEOREM 5.** *Recall the notations in Theorem 4. Under Assumptions 1 and 2, there exists $\beta_s > 0$ such that for $\beta > \beta_s$,*

$$ \left\{ \begin{array}{l} \lambda \xrightarrow{a.s.} \lambda^\infty(\beta), \\ |\langle \mathbf{u}, \mathbf{x} \rangle| \xrightarrow{a.s.} \frac{1}{\sqrt{\alpha_2(\lambda^\infty(\beta))\alpha_3(\lambda^\infty(\beta))}}, \\ |\langle \mathbf{v}, \mathbf{y} \rangle| \xrightarrow{a.s.} \frac{1}{\sqrt{\alpha_1(\lambda^\infty(\beta))\alpha_3(\lambda^\infty(\beta))}}, \\ |\langle \mathbf{w}, \mathbf{z} \rangle| \xrightarrow{a.s.} \frac{1}{\sqrt{\alpha_1(\lambda^\infty(\beta))\alpha_2(\lambda^\infty(\beta))}} \end{array} \right. $$

where $\alpha_i(z) \equiv \frac{\beta}{z+g(z)-g_i(z)}$ and $\lambda^\infty(\beta)$ satisfies $f(\lambda^\infty(\beta), \beta) = 0$ with $f(z, \beta) = z + g(z) - \frac{\beta}{\alpha_1(z)\alpha_2(z)\alpha_3(z)}$. Besides, for $\beta \in [0, \beta_s]$, $\lambda^\infty$ is bounded⁹ by an order one constant and $|\langle \mathbf{u}, \mathbf{x} \rangle|, |\langle \mathbf{v}, \mathbf{y} \rangle|, |\langle \mathbf{w}, \mathbf{z} \rangle| \xrightarrow{a.s.} 0$.

**PROOF.** See Appendix B.5. □

**REMARK 6.** *A more compact expression for the alignments is provided in Theorem 8.*
<figure><img src="image_7.png" /><figcaption>FIG 5. Asymptotic singular value and alignments of $\mathbf{T}$ as predicted by Theorem 5 for $c_1 = \frac{1}{6}$, $c_2 = \frac{1}{3}$ and $c_3 = \frac{1}{2}$.</figcaption></figure>

Figure 5 depicts the predicted asymptotic dominant singular value $\lambda^\infty(\beta)$ of $\mathbf{T}$ and the corresponding alignments from Theorem 5. As we can see, the result of Theorem 5 predicts that a non-zero correlation between the population signals $\mathbf{x}, \mathbf{y}, \mathbf{z}$ and their estimated counterparts is possible only when $\beta > \beta_s \approx 1.1134$, which corresponds to the value after which $\lambda^\infty(\beta)$ starts to increase with $\beta$.

**REMARK 7.** Note that, given $g(z)$, the inverse formula expressing $\beta$ in terms of $\lambda^\infty$ is explicit. Specifically, we have $\beta(\lambda^\infty) = \sqrt{\frac{\prod_{i=1}^3 (\lambda^\infty + g(\lambda^\infty) - g_i(\lambda^\infty))}{\lambda^\infty + g(\lambda^\infty)}}$. In particular, this inverse formula provides an estimator for the SNR $\beta$ given the largest singular value $\lambda$ of $\mathbf{T}$.

**3.6. Cubic 3-order tensors: case $c_1 = c_2 = c_3 = \frac{1}{3}$**. In this section, we study the particular case where all the tensor dimensions are equal. As such, the three alignments $\langle \mathbf{u}, \mathbf{x} \rangle, \langle \mathbf{v}, \mathbf{y} \rangle, \langle \mathbf{w}, \mathbf{z} \rangle$ converge almost surely to the same quantity. In this case, the almost sure limits of $\lambda$ and $\langle \mathbf{u}, \mathbf{x} \rangle, \langle \mathbf{v}, \mathbf{y} \rangle, \langle \mathbf{w}, \mathbf{z} \rangle$ can be obtained explicitly in terms of the signal strength $\beta$ as per the following corollary of Theorem 5.

**COROLLARY 3.** *Under Assumptions 2 and 1 with $c_1 = c_2 = c_3 = \frac{1}{3}$, for $\beta > \beta_s = \frac{2\sqrt{3}}{3}$*
$$ \lambda \xrightarrow{\text{a.s.}} \lambda^\infty(\beta) = \sqrt{\frac{\beta^2}{2} + 2 + \frac{\sqrt{3}\sqrt{(3\beta^2-4)^3}}{18\beta}}, $$
$$ |\langle \mathbf{u}, \mathbf{x} \rangle|, |\langle \mathbf{v}, \mathbf{y} \rangle|, |\langle \mathbf{w}, \mathbf{z} \rangle| \xrightarrow{\text{a.s.}} \frac{\sqrt{9\beta^2-12+\frac{\sqrt{3}\sqrt{(3\beta^2-4)^3}}{\beta}} + \sqrt{9\beta^2+36+\frac{\sqrt{3}\sqrt{(3\beta^2-4)^3}}{\beta}}}{6\sqrt{2\beta}}. $$
Besides, for $\beta \in [0, \frac{2\sqrt{3}}{3}]$, $\lambda \xrightarrow{\text{a.s.}} \lambda^\infty \le 2\sqrt{\frac{2}{3}}$ and $|\langle \mathbf{u}, \mathbf{x} \rangle|, |\langle \mathbf{v}, \mathbf{y} \rangle|, |\langle \mathbf{w}, \mathbf{z} \rangle| \xrightarrow{\text{a.s.}} 0$.

**PROOF.** See Appendix B.6. □

Figure 6 provides plots of the almost sure limits of the singular value and alignments when $\mathbf{T}$ is cubic as per Corollary 3 (see Subsection A.2 for simulations supporting the above result). In particular, this result predicts a possible correlation between the singular vectors and the underlying signal components above the value $\beta_s = \frac{2\sqrt{3}}{3} \approx 1.154$ with corresponding

<footnote n="8">⁸Yielding non-vanishing terms in the expression of $\lambda^\infty(\beta)$.</footnote>

<footnote n="9">⁹Such bound might be computed numerically and corresponds to the right edge of the limiting spectral measure in Theorem 4.</footnote>
<figure><img src="image_8.png" /><figcaption>FIG 6. Asymptotic singular value and alignments of T for $c_1 = c_2 = c_3 = \frac{1}{3}$ as per Corollary 3.</figcaption></figure>

singular value $\lambda_s = 2\sqrt{\frac{2}{3}} \approx 1.633$ with an alignment $a_s = \frac{\sqrt{2}}{2} \approx 0.707$. In addition, we can easily check from the formulas above that $\frac{\lambda^\infty(\beta)}{\beta} \to 1$ and $|\langle u, x \rangle| \to 1$ for large values of $\beta$. Besides, for values of $\beta$ around the value $\beta_s$, the expression of $\lambda^\infty(\beta)$ admits the following expansion

$$ \lambda^{\infty}(\beta) = 2\sqrt{\frac{2}{3}} + \frac{\sqrt{2}}{4}(\beta - \beta_s) + \frac{\sqrt{2}3^{\frac{1}{4}}}{4}(\beta - \beta_s)^{\frac{3}{2}} + \frac{3\sqrt{6}}{64}(\beta - \beta_s)^2 + o((\beta - \beta_s)^2), $$

whereas the corresponding alignment expands as

$$ \frac{\sqrt{2}}{2} + \frac{\sqrt{2}3^{\frac{1}{4}}}{4}\sqrt{\beta - \beta_s} - \frac{\sqrt{6}}{16}(\beta - \beta_s) - \frac{\sqrt{2}3^{\frac{3}{4}}}{16}(\beta - \beta_s)^{\frac{3}{2}} + \frac{21\sqrt{2}}{256}(\beta - \beta_s)^2 + o((\beta - \beta_s)^2). $$

**4. Random tensors meet random matrices.** In this section, we investigate the application of Theorem 5 to the particular case of spiked random matrices. Indeed, for instance when $p = 1$ the spiked tensor model in Eq. (10) becomes a spiked matrix model which we will now denote as

$$ (23) \qquad M = \beta x y^\top + \frac{1}{\sqrt{N}} X, $$

where again $x$ and $y$ are on the unit spheres $S^{m-1}$ and $S^{n-1}$ respectively, $N = m + n$ and $X$ is a Gaussian noise matrix with i.i.d. entries $X_{ij} \sim \mathcal{N}(0, 1)$.

Our approach does not apply directly to the matrix $M$ above. Indeed, the singular vectors $u, v$ of $M$ corresponding to its largest singular value $\lambda$ satisfy the identities

$$ (24) \qquad Mv = \lambda u, \quad M^\top u = \lambda v, \quad \lambda = u^\top Mv, $$

and deriving $u, v$ w.r.t. an entry $X_{ij}$ of $X$ will result in

$$ (25) \qquad \left( \begin{bmatrix} 0_{m \times m} & M \\ M^\top & 0_{n \times n} \end{bmatrix} - \lambda I_N \right) \begin{bmatrix} \frac{\partial u}{\partial X_{ij}} \\ \frac{\partial v}{\partial X_{ij}} \end{bmatrix} = - \frac{1}{\sqrt{N}} \begin{bmatrix} v_j (e_i^m - u_i u) \\ u_i (e_j^n - v_j v) \end{bmatrix}. $$

Now since $\lambda$ is an eigenvalue of $\begin{bmatrix} 0_{m \times m} & M \\ M^\top & 0_{n \times n} \end{bmatrix}$, the matrix $\left( \begin{bmatrix} 0_{m \times m} & M \\ M^\top & 0_{n \times n} \end{bmatrix} - \lambda I_N \right)$ is not invertible, which makes our approach inapplicable for $d = 2$. However, we can retrieve the behavior of the spiked matrix model in Eq. (23) by considering an order 3 tensor and setting for instance $c_3 \to 0$ as we will discuss subsequently.
<figure><img src="image_9.png" /><figcaption>FIG 7. Limiting spectral measure of $\frac{1}{\sqrt{N}}\Phi_3(\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c})$ for $c_1 = c, c_2 = 1 - c$ as per Corollary 4.</figcaption></figure>

4.1. Limiting spectral measure. Given the result of Theorem 5, the asymptotic singular value and alignments for the spiked matrix model in Eq. (23) correspond to the particular case $c_3 \to 0$. We start by characterizing the corresponding limiting spectral measure, we have the following corollary of Theorem 4 when $c_3 \to 0$.

COROLLARY 4. Given the setting and notations of Theorem 3, under Assumptions 2 and 1 with $c_1 = c, c_2 = 1 - c$ for $c \in (0, 1)$, the empirical spectral measure of $\frac{1}{\sqrt{N}}\Phi_3(\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c})$ converges weakly almost surely to a deterministic measure $\nu$ defined on the support $S(\nu) = [-\sqrt{1+2\sqrt{\eta}}, -\sqrt{1-2\sqrt{\eta}}] \cup [\sqrt{1-2\sqrt{\eta}}, \sqrt{1+2\sqrt{\eta}}]$ with $\eta = c(1-c)$, whose density function writes as

$$ \nu(dx) = \frac{1}{\pi x} \sin\left(\frac{\arctan_2(0, q_c(x))}{2}\right) \sqrt{|q_c(x)|} \operatorname{sign}\left(\frac{\sin\left(\frac{\arctan_2(0, q_c(x))}{2}\right)}{x}\right) dx + (1 - 2\min(c, 1-c)) \delta(x), $$
where $q_c(x) = (x^2 - 1)^2 + 4c(c-1)$. And the corresponding Stieltjes transform writes as

$$ g(z) = -z + \frac{\sqrt{q_c(z)}}{z}, \quad \text{for } z \in \mathbb{C} \setminus S(\nu). $$

PROOF. See Appendix B.7.

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

REMARK 8. Corollary 4 describes the limiting spectral measure of $\frac{1}{\sqrt{m+n}}\Phi_3(\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c})$ for $\mathbf{X} \sim T_{m,n,1}(\mathcal{N}(0, 1))$, $(\mathbf{a}, \mathbf{b}) \in \mathbb{S}^{m-1} \times \mathbb{S}^{n-1}$ and $c = 1$ (a scalar), which is equivalent to random matrices of the form $\frac{1}{\sqrt{m+n}} \begin{bmatrix} \mathbf{0}_{m \times m} & \mathbf{X} \\ \mathbf{X}^{\top} & \mathbf{0}_{n \times n} \end{bmatrix}$ with $\mathbf{X} \sim \mathbb{M}_{m,n}(\mathcal{N}(0, 1))$.

Figure 7 depicts the limiting spectral measures as per Corollary 4 for various values of $c$ (see also simulations in Appendix A.1). In particular, the limiting measure is a semi-circle law for square matrices (i.e., $c = \frac{1}{2}$), while it decomposes into a two-mode distribution<fnref n="10" /> for $c \in (0, 1) \setminus \{\frac{1}{2}\}$ due to the fact that the underlying random matrix model is not full rank ($\Phi_3(\mathbf{X}, \mathbf{a}, \mathbf{b}, 1)$ if of rank $2\min(m, n)$).

<footnote n="10">¹⁰With a Dirac $\delta(x)$ at 0 weighted by $1 - 2\min(c, 1-c)$, not shown in Figure 7.</footnote>
<figure><img src="image_11.png" /><figcaption>FIG 8. Asymptotic singular value and alignments of $\mathbf{T}$ (and of $\mathbf{M}$ in Eq. (23)) for $c_1 = c$ and $c_2 = 1 - c$ as per Corollary 5, for different values of $c \in \{\frac{1}{2}, \frac{1}{10}, \frac{1}{50}\}$.</figcaption></figure>

4.2. Limiting singular value and alignments. Given the limiting Stieltjes transform from the previous subsection (Corollary 4), the limiting largest singular value of $\mathbf{M}$ (in Eq. (23)) and the alignments of the corresponding singular vectors are obtained thanks to Theorem 5, yielding the following corollary.

**COROLLARY 5.** Under Assumption 1 with $c_1 = c, c_2 = 1 - c$, we have for $\beta > \beta_s = \sqrt[4]{c(1-c)}$

$$ \lambda \xrightarrow{a.s.} \lambda^{\infty}(\beta) = \sqrt{\beta^2 + 1 + \frac{c(1-c)}{\beta^2}}, \\ |\langle \mathbf{u}, \mathbf{x} \rangle| \xrightarrow{a.s.} \frac{1}{\kappa(\beta, c)}, \quad |\langle \mathbf{v}, \mathbf{y} \rangle| \xrightarrow{a.s.} \frac{1}{\kappa(\beta, 1-c)}, $$

where $\mathbf{u} \in \mathbb{S}^{m-1}$, $\mathbf{v} \in \mathbb{S}^{n-1}$ are the singular vectors of $\mathbf{M}$ corresponding to its largest singular value $\lambda$ and $\kappa(\beta, c)$ is given by

$$ \kappa(\beta, c) = \beta \sqrt{\frac{\beta^2 (\beta^2 + 1) - c(c-1)}{(\beta^4 + c(c-1))(\beta^2 + 1 - c)}} $$

Besides, for $\beta \in [0, \beta_s]$, $\lambda \xrightarrow{a.s.} \sqrt{1 + 2\sqrt{c(1-c)}}$ and $|\langle \mathbf{u}, \mathbf{x} \rangle|, |\langle \mathbf{v}, \mathbf{y} \rangle| \xrightarrow{a.s.} 0$.

**PROOF.** See Appendix B.8. $\square$

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

Figure 8 provides the curves of the asymptotic singular value and alignments of Corollary 5 (see Subsection A.1 for simulations supporting the above result). Unlike the tensor case from the previous section, we see that the asymptotic alignments $a_x^\infty(\beta)$, $a_y^\infty(\beta)$ are continuous and a positive alignment is observed for $\beta > \beta_s$ which corresponds to the classical BBP phase transition of spiked random matrix models [2, 5, 19].

**REMARK 9 (Application to tensor unfolding).** The tensor unfolding method consists in estimating the spike components of a given tensor $\mathbf{T} = \gamma \mathbf{x}^{(1)} \otimes \cdots \otimes \mathbf{x}^{(d)} + \frac{1}{\sqrt{n}} \mathbf{X} \in \mathbb{R}^{n \times n}$.
$\mathbb{T}_n^d$ with $\mathbf{X} \sim \mathbb{T}_n^d(\mathcal{N}(0,1))$ by applying an SVD to its unfolded matrices (for $i \in [d]$) $\mathrm{Mat}_i(\mathbf{T}) = \gamma \mathbf{x}^{(i)}(\mathbf{y}^{(i)})^\top + \frac{1}{\sqrt{n}}\mathbf{X}_i$ of size $n \times n^{d-1}$. Applying Corollary 5 to this case predicts a phase transition at $\beta_s = \sqrt[4]{c(1-c)}$ with$^{11}$ $c = \frac{n}{n+n^{d-1}} = \frac{1}{1+n^{d-2}} \to 0$, which yields $\beta_s = \sqrt[4]{\frac{n^{d-2}}{(1+n^{d-2})^2}}$. After re-scaling the noise component (by multiplying Eq. (23) by $\sqrt{1+\frac{m}{n}}$ with $m = n^{d-1}$) this yields $\gamma > n^{\frac{d-2}{4}}$, i.e., the phase transition for tensor unfolding obtained by [3] (see Theorem 3.3 therein). More generally, for any order-d tensor $\mathbf{T}$ of arbitrary dimensions $n_1 \times \cdots \times n_d$, the tensor unfolding method succeeds provided that $\beta \ge (\prod_{i=1}^d n_i)^{1/4} / \sqrt{\sum_{i=1}^d n_i}$.

**5. Generalization to arbitrary d-order tensors.** We now show that our approach can be generalized straightforwardly to the d-order spiked tensor model of Eq. (1). Indeed, from Eq. (4), the $\ell_2$-singular value $\lambda$ and vectors $\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)} \in \mathbb{S}^{n_1-1} \times \dots \times \mathbb{S}^{n_d-1}$, corresponding to the best rank-one approximation $\lambda \mathbf{u}^{(1)} \otimes \dots \otimes \mathbf{u}^{(d)}$ of the d-order tensor $\mathbf{T}$ in Eq. (1), satisfy the identities

$$ (26) \qquad \begin{cases} \mathbf{T}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(i-1)}, \dots, \mathbf{u}^{(i)}, \dots, \mathbf{u}^{(i+1)}, \dots, \mathbf{u}^{(d)}) = \lambda \mathbf{u}^{(i)}, \\ \lambda = \mathbf{T}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)}) = \sum_{i_1, \dots, i_d} u_{i_1}^{(1)} \cdots u_{i_d}^{(d)} T_{i_1, \dots, i_d}. \end{cases} $$

**5.1. Associated random matrix ensemble.** Let $\mathbf{T}^{ij} \in \mathbb{M}_{n_i, n_j}$ denote the matrix obtained by contracting the tensor $\mathbf{T}$ with the singular vectors in $\{\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)}\} \setminus \{\mathbf{u}^{(i)}, \mathbf{u}^{(j)}\}$, i.e.,

$$ (27) \qquad \mathbf{T}^{ij} \equiv \mathbf{T}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(i-1)}, \dots, \mathbf{u}^{(i+1)}, \dots, \mathbf{u}^{(j-1)}, \dots, \mathbf{u}^{(j)}, \dots, \mathbf{u}^{(d)}). $$

As in the order 3 case, from Eq. (26), the derivatives of the singular vectors $\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)}$ with respect to the entry $X_{i_1, \dots, i_d}$ of the noise tensor $\mathbf{X}$ express as

$$ (28) \qquad \begin{bmatrix} \frac{\partial \mathbf{u}^{(1)}}{\partial X_{i_1, \dots, i_d}} \\ \vdots \\ \frac{\partial \mathbf{u}^{(d)}}{\partial X_{i_1, \dots, i_d}} \end{bmatrix} = -\frac{1}{\sqrt{N}} \begin{bmatrix} \mathbf{0}_{n_1 \times n_1} & \mathbf{T}^{12} & \mathbf{T}^{13} & \dots & \mathbf{T}^{1d} \\ (\mathbf{T}^{12})^\top & \mathbf{0}_{n_2 \times n_2} & \mathbf{T}^{23} & \dots & \mathbf{T}^{2d} \\ (\mathbf{T}^{13})^\top & (\mathbf{T}^{23})^\top & \mathbf{0}_{n_3 \times n_3} & \dots & \mathbf{T}^{3d} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ (\mathbf{T}^{1d})^\top & (\mathbf{T}^{2d})^\top & (\mathbf{T}^{3d})^\top & \dots & \mathbf{0}_{n_d \times n_d} \end{bmatrix} - \lambda \mathbf{I}_N^{-1} \begin{bmatrix} \prod_{\ell \in \{2, \dots, d\}} u_{i_\ell}^{(\ell)} (e_{i_1}^{n_1} - u_{i_1}^{(i_1)} \mathbf{u}^{(i_1)}) \\ \vdots \\ \prod_{\ell \in \{1, \dots, d-1\}} u_{i_\ell}^{(\ell)} (e_{i_d}^{n_d} - u_{i_d}^{(i_d)} \mathbf{u}^{(i_d)}) \end{bmatrix} $$

where $N = \sum_{i \in [d]} n_i$, and the derivative of $\lambda$ w.r.t. $X_{i_1, \dots, i_d}$ writes as

$$ (29) \qquad \frac{\partial \lambda}{\partial X_{i_1, \dots, i_d}} = \frac{1}{\sqrt{N}} \prod_{\ell \in [d]} u_{i_\ell}^{(\ell)}. $$

As such, the associated random matrix model of $\mathbf{T}$ is the matrix appearing in the resolvent in Eq. (28). More generally, the d-order block-wise tensor contraction ensemble $\mathcal{B}_d(\mathbf{X})$ for $\mathbf{X} \sim \mathbb{T}_{n_1, \dots, n_d}(\mathcal{N}(0,1))$ is defined as

$$ (30) \qquad \mathcal{B}_d(\mathbf{X}) \equiv \left\{ \Phi_d(\mathbf{X}, \mathbf{a}^{(1)}, \dots, \mathbf{a}^{(d)}) \mid (\mathbf{a}^{(1)}, \dots, \mathbf{a}^{(d)}) \in \mathbb{S}^{n_1-1} \times \dots \times \mathbb{S}^{n_d-1} \right\} $$

<footnote>$^{11}$In our assumptions, we assume that $c$ is some constant that does not depend on the tensor dimensions. Still, we believe that it can be relaxed in the same vein as in [4].</footnote>
**Algorithm 1 Fixed point equation to compute the Stieltjes transform in Theorem 6**

**Input:** $z \in \mathbb{R} \setminus S(\nu)$ and tensor dimension ratios $c = [c_1, \dots, c_d]^\top \in (0, 1)^d$

**Output:** $g \in \mathbb{R}$, $g = [g_1, \dots, g_d]^\top \in \mathbb{R}^d$

**Repeat**

$$g \leftarrow \frac{g+z}{2} - \frac{\sqrt{4c+(g+z)^2}}{2}$$

$$g \leftarrow \sum_{i=1}^{d} g_i$$

until convergence of $g$

▷ Element-wise vector operation.

where $\Phi_d$ is the mapping

$$\Phi_d: \mathbb{T}_{n_1, \dots, n_d} \times \mathbb{S}^{n_1-1} \times \dots \times \mathbb{S}^{n_d-1} \to \mathbb{M}_{\sum_i n_i} \\ (31) \qquad (\mathbf{X}, \mathbf{a}^{(1)}, \dots, \mathbf{a}^{(d)}) \mapsto \begin{bmatrix} \mathbf{0}_{n_1 \times n_1} & \mathbf{X}^{12} & \mathbf{X}^{13} & \dots & \mathbf{X}^{1d} \\ (\mathbf{X}^{12})^\top & \mathbf{0}_{n_2 \times n_2} & \mathbf{X}^{23} & \dots & \mathbf{X}^{2d} \\ (\mathbf{X}^{13})^\top & (\mathbf{X}^{23})^\top & \mathbf{0}_{n_3 \times n_3} & \dots & \mathbf{X}^{3d} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ (\mathbf{X}^{1d})^\top & (\mathbf{X}^{2d})^\top & (\mathbf{X}^{3d})^\top & \dots & \mathbf{0}_{n_d \times n_d} \end{bmatrix},$$

$$\text{with } \mathbf{X}^{ij} \equiv \mathbf{X}(\mathbf{a}^{(1)}, \dots, \mathbf{a}^{(i-1)}, \dots, \mathbf{a}^{(i+1)}, \dots, \mathbf{a}^{(j-1)}, \dots, \mathbf{a}^{(j+1)}, \dots, \mathbf{a}^{(d)}) \in \mathbb{M}_{n_i, n_j}.$$

**REMARK 10.** As in the order 3 case, to ensure the existence of the matrix inverse in Eq. (28), we need to suppose that there exists a tuple $(\lambda_*, \mathbf{u}_*^{(1)}, \dots, \mathbf{u}_*^{(d)})$ verifying the identities in Eq. (26) such that $\lambda_*$ is not an eigenvalue of $\Phi_d(\mathbf{T}, \mathbf{u}_*^{(1)}, \dots, \mathbf{u}_*^{(d)})$.

5.2. Limiting spectral measure of block-wise d-order tensor contractions. In this section, we characterize the limiting spectral measure of the ensemble $\mathcal{B}_d(\mathbf{X})$ for $\mathbf{X} \sim \mathbb{T}_{n_1, \dots, n_d}(\mathcal{N}(0, 1))$ in the limit when all tensor dimensions grow as per the following assumption.

**ASSUMPTION 3.** For all $i \in [d]$, assume that $n_i \to \infty$ with $\frac{n_i}{\sum_j n_j} \to c_i \in (0, 1)$.

We thus have the following result which characterizes the spectrum of $\frac{1}{\sqrt{N}}\Phi_d(\mathbf{X}, \mathbf{a}^{(1)}, \dots, \mathbf{a}^{(d)})$ for any deterministic unit norm vectors $\mathbf{a}^{(1)}, \dots, \mathbf{a}^{(d)}$.

**THEOREM 6.** Let $\mathbf{X} \sim \mathbb{T}_{n_1, \dots, n_d}(\mathcal{N}(0, 1))$ be a sequence of random asymmetric Gaussian tensors and $(\mathbf{a}^{(1)}, \dots, \mathbf{a}^{(d)}) \in \mathbb{S}^{n_1-1} \times \dots \times \mathbb{S}^{n_d-1}$ a sequence of deterministic vectors of increasing dimensions, following Assumption 3. Then the empirical spectral measure of $\frac{1}{\sqrt{N}}\Phi_d(\mathbf{X}, \mathbf{a}^{(1)}, \dots, \mathbf{a}^{(d)})$ converges weakly almost surely to a deterministic measure $\nu$ whose Stieltjes transform $g(z)$ is defined as the solution to the equation $g(z) = \sum_{i=1}^d g_i(z)$ such that $\mathfrak{S}[g(z)] > 0$ for $\mathfrak{S}[z] > 0$ where, for $i \in [d]$ $g_i(z)$ satisfies $g_i^2(z) - (g(z) + z)g_i(z) - c_i = 0$ for $z \in \mathbb{C} \setminus S(\nu)$.

**PROOF.** See Appendix B.9. □

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

Algorithm 1 provides a pseudo-code to compute the Stieltjes transform $g(z)$ in Theorem 6 through an iterative solution to the fixed point equation $g(z) = \sum_{i=1}^d g_i(z)$. In particular, as for the order 3 case, when all the tensor dimensions are equal (i.e., $c_i = \frac{1}{d}$ for all $i \in [d]$), the spectral measure of $\frac{1}{\sqrt{N}}\Phi_d(\mathbf{X}, \mathbf{a}^{(1)}, \dots, \mathbf{a}^{(d)})$ converges to a semi-circle law. We have the following corollary of Theorem 6.
<figure><img src="image_14.png" /><figcaption>FIG 9. Spectrum of $\frac{1}{\sqrt{N}}\Phi_4(\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c}, \mathbf{d})$ with $\mathbf{X} \sim \mathbb{T}_{n_1, \dots, n_4}(\mathcal{N}(0, 1))$ and $\mathbf{a}, \mathbf{b}, \mathbf{c}, \mathbf{d}$ independently sampled from the unit spheres $\mathbb{S}^{n_1-1}, \dots, \mathbb{S}^{n_4-1}$ respectively. In black the semi-circle law as per Theorem 6.</figcaption></figure>

**COROLLARY 6.** With the setting of Theorem 6. Under Assumption 3 with $c_i = \frac{1}{d}$ for all $i \in [d]$, the empirical spectral measure of $\frac{1}{\sqrt{N}}\Phi_d(\mathbf{X}, \mathbf{a}^{(1)}, \dots, \mathbf{a}^{(d)})$ converges weakly almost surely to the semi-circle distribution support $S(\nu) = [-2\sqrt{\frac{d-1}{d}}, 2\sqrt{\frac{d-1}{d}}]$, whose density and Stieltjes transform write respectively as

$$ \nu(dx) = \frac{d}{2(d-1)\pi} \sqrt{\left(\frac{4(d-1)}{d} - x^2\right)^+}, \quad g(z) = \frac{-zd + d\sqrt{z^2 - \frac{4(d-1)}{d}}}{2(d-1)}, \quad \text{where } z \in \mathbb{C} \setminus S(\nu). $$

**PROOF.** See Appendix B.10. $\square$

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

Figure 9 provides an illustration of the convergence in law (when the dimension $n_i$ grow large) of the spectrum of $\frac{1}{\sqrt{N}}\Phi_4(\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c}, \mathbf{d})$ with a 4th-order tensor $\mathbf{X} \sim \mathbb{T}_{n_1, \dots, n_4}(\mathcal{N}(0, 1))$ and $\mathbf{a}, \mathbf{b}, \mathbf{c}, \mathbf{d}$ independently sampled from the unit spheres $\mathbb{S}^{n_1-1}, \dots, \mathbb{S}^{n_4-1}$ respectively.

**REMARK 11.** $\frac{1}{\sqrt{N}}\Phi_d(\mathbf{X}, \mathbf{a}^{(1)}, \dots, \mathbf{a}^{(d)})$ is almost surely of rank $\sum_i \min(n_i, \sum_{j \neq i} n_j)$. As in the order 3 case, when $\frac{1}{\sqrt{N}}\Phi_d(\mathbf{X}, \mathbf{a}^{(1)}, \dots, \mathbf{a}^{(d)})$ is almost surely full rank, its spectral measure converges to a semi-circle law with connected support, while in general the limiting spectral measure has unconnected support.

The result of Theorem 6 still holds for $\Phi_d(\mathbf{T}, \mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)})$ where $\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)}$ stand for the singular vectors of $\mathbf{T}$ defined through Eq. (26). Indeed, the statistical dependencies between these singular vectors and the noise tensor $\mathbf{X}$ do not affect the convergence of the spectrum of $\Phi_d(\mathbf{T}, \mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)})$ to the limiting measure described by Theorem 6. We further need the following assumption.

**ASSUMPTION 4.** We assume that there exists a sequence of critical points $(\lambda_*, \mathbf{u}_*^{(1)}, \dots, \mathbf{u}_*^{(d)})$ satisfying Eq. (26) such that $\lambda_* \xrightarrow{a.s.} \lambda^\infty(\beta)$, $|\langle \mathbf{u}_*^{(i)}, \mathbf{x}^{(i)} \rangle| \xrightarrow{a.s.} a_{\mathbf{x}^{(i)}}^\infty(\beta)$ with $\lambda^\infty(\beta) \notin S(\nu)$ and $a_{\mathbf{x}^{(i)}}^\infty(\beta) > 0$.

**THEOREM 7.** Let $\mathbf{T}$ be a sequence of spiked random tensors defined as in Eq. (1). Under Assumptions 3 and 4, the empirical spectral measure of $\Phi_d(\mathbf{T}, \mathbf{u}_*^{(1)}, \dots, \mathbf{u}_*^{(d)})$ converges weakly almost surely to a deterministic measure $\nu$ whose Stieltjes transform $g(z)$ is defined as the solution to the equation $g(z) = \sum_{i=1}^d g_i(z)$ such that $\mathfrak{F}[g(z)] > 0$ for $\mathfrak{F}[z] > 0$ where, for $i \in [d]$ $g_i(z)$ satisfies $g_i^2(z) - (g(z) + z)g_i(z) - c_i = 0$ for $z \in \mathbb{C} \setminus S(\nu)$.
<figure><img src="image_16.png" /><figcaption>FIG 10. Spectrum of $\Phi_4(\mathbf{T}, \mathbf{u}^{(1)}, \dots, \mathbf{u}^{(4)})$ at iterations 0, 5 and at convergence of the power method (see Algorithm 3) applied to $\mathbf{T}$, for $n_1 = n_2 = n_3 = n_4 = 50$ and $\beta = 0$.</figcaption></figure>

PROOF. See Appendix B.12.

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

Figure 10 depicts the spectrum of $\Phi_4(\mathbf{T}, \mathbf{u}^{(1)}, \dots, \mathbf{u}^{(4)})$ for an order 4 tensor $\mathbf{T}$ with $\beta = 0$. As we saw previously, an isolated eigenvalue pops out from the continuous bulk of $\Phi_4(\mathbf{T}, \mathbf{u}^{(1)}, \dots, \mathbf{u}^{(4)})$ because of the statistical dependencies between the tensor noise $\mathbf{X}$ and the singular vectors $\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(4)}$. More generally, for an order-$d$ tensor $\mathbf{T}$, the spectrum of $\Phi_d(\mathbf{T}, \mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)})$ admits an isolated spike at the value $(d-1)\lambda$ independently of the signal strength $\beta$. Indeed, $(d-1)\lambda$ is an eigenvalue of $\Phi_d(\mathbf{T}, \mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)})$ with corresponding eigenvector the concatenation of the singular vectors $\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)}$, i.e.,

$$ \Phi_d(\mathbf{T}, \mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)}) \begin{bmatrix} \mathbf{u}^{(1)} \\ \vdots \\ \mathbf{u}^{(d)} \end{bmatrix} = \begin{bmatrix} \sum_{i \neq 1} \mathbf{T}(:, \mathbf{u}^{(2)}, \dots, \mathbf{u}^{(i-1)}, \mathbf{u}^{(i+1)}, \dots, \mathbf{u}^{(d)}) \mathbf{u}^{(i)} \\ \vdots \\ \sum_{i \neq d} \mathbf{T}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(i-1)}, \mathbf{u}^{(i+1)}, \dots, \mathbf{u}^{(d-1)}, \mathbf{u}^{(d)})^{\top} \mathbf{u}^{(i)} \end{bmatrix} = (d-1)\lambda \begin{bmatrix} \mathbf{u}^{(1)} \\ \vdots \\ \mathbf{u}^{(d)} \end{bmatrix}. $$

5.3. Asymptotic singular value and alignments of hyper-rectangular tensors. Similarly to the 3-order case studied previously, when the dimensions of $\mathbf{T}$ grow large at a same rate, its singular value $\lambda$ and the corresponding alignments $\langle \mathbf{u}^{(i)}, \mathbf{x}^{(i)} \rangle$ for $i \in [d]$ concentrate almost surely around some deterministic quantities which we denote $\lambda^\infty(\beta) \equiv \lim_{N \to \infty} \lambda$ and $a_{\mathbf{x}^{(i)}}^\infty(\beta) = \lim_{N \to \infty} |\langle \mathbf{u}^{(i)}, \mathbf{x}^{(i)} \rangle|$ respectively. Applying again Stein's Lemma (Lemma 1) to the identities in Eq. (26), we obtain the following theorem which characterizes the aforementioned deterministic limits.

THEOREM 8. Recall the notations in Theorem 7. Under Assumptions 3 and 4, for $d \ge 3$, there exists $\beta_s > 0$ such that for $\beta > \beta_s$,

$$ \lambda \xrightarrow{a.s.} \lambda^\infty(\beta), \\ |\langle \mathbf{x}^{(i)}, \mathbf{u}^{(i)} \rangle| \xrightarrow{a.s.} q_i(\lambda^\infty(\beta)), $$

where $q_i(z)$ is given by $q_i(z) = \sqrt{1 - \frac{g_i^2(z)}{c_i}}$ and $\lambda^\infty(\beta)$ satisfies $f(\lambda^\infty(\beta), \beta) = 0$ with $f(z, \beta) = z + g(z) - \beta \prod_{i=1}^d q_i(z)$. Besides, for $\beta \in [0, \beta_s]$, $\lambda^\infty$ is bounded (in particular when $c_i = \frac{1}{d}$ for all $i \in [d]$, $\lambda \xrightarrow{a.s.} \lambda^\infty \le 2\sqrt{\frac{d-1}{d}}$) and $|\langle \mathbf{x}^{(i)}, \mathbf{u}^{(i)} \rangle| \xrightarrow{a.s.} 0$.

PROOF. See Appendix B.13.

<figure><img src="image_18.png" /></figure>
<caption>**Algorithm 2** Compute alignments as per Theorem 8</caption>

**Input:** SNR $\beta \in \mathbb{R}_+$ and tensor dimension ratios $c = [c_1, \dots, c_d]^\top \in (0, 1)^d$
**Output:** Asymptotic singular value $\lambda^\infty$ and corresponding alignments $\boldsymbol{a} = [a_{\boldsymbol{x}(1)}^\infty, \dots, a_{\boldsymbol{x}(d)}^\infty]^\top \in [0, 1]^d$
Set $\lambda^\infty$ as the solution of $f(z, \beta) = 0$ in $z$.
Set the alignments as $a_{\boldsymbol{x}(i)}^\infty \leftarrow \sqrt{1 - \frac{g_i^2(\lambda^\infty)}{c_i}}$ with $g_i(z)$ computed by Algorithm 1 for $z = \lambda^\infty$.

<figure><img src="image_19.png" /><figcaption>FIG 11. Asymptotic singular value and alignments of $\mathbf{T}$ as predicted by Theorem 8 for $c_1 = \frac{1}{10}, c_2 = \frac{1}{6}, c_3 = \frac{1}{4}$ and $c_4 = 1 - (c_1 + c_2 + c_3)$.</figcaption></figure>

**REMARK 12.** As in the order 3 case, note that the inverse formula expressing $\beta$ in terms of $\lambda^\infty$ is explicit. Specifically, we have $\beta(\lambda^\infty) = \frac{\lambda^\infty + g(\lambda^\infty)}{\prod_{i=1}^d q_i(\lambda^\infty)}$. In particular, this inverse formula provides an estimator for the SNR $\beta$ given the largest singular value $\lambda$ of $\mathbf{T}$. Algorithm 2 provides a pseudo-code to compute the asymptotic alignments.

Figure 11 depicts the asymptotic singular value and alignments for an order 4 tensor as per Theorem 8. As discussed previously, the predicted alignments are discontinuous and a strictly positive correlation between the singular vectors and the underlying signals is possible for the considered local optimum of the ML estimator only above the minimal signal strength $\beta_s \approx 1.234$ in the shown example. In the case of hyper-cubic tensors, i.e., all the dimensions $n_i$ are equal, Figure 12 depicts the minimal SNR values in terms of the tensor order $d$ and the corresponding asymptotic singular value and alignments. As such, when the tensor order increases the minimal SNR and singular value converge respectively to $\approx 1.6$ and 2, while the corresponding alignment<fnref n="12" /> gets closer to 1. The expressions of $\beta_s$ and $a^\infty(\beta_s)$ for hyper-cubic tensors of order $d$ are explicitly given as

$$ (33) \qquad \beta_s = \sqrt{\frac{d-1}{d}} \left(\frac{d-2}{d-1}\right)^{1-\frac{d}{2}}, \quad \lim_{\beta \to \beta_s} a^\infty(\beta) = \sqrt{\frac{d-2}{d-1}}. $$

**6. Generalization to rank $r$ tensor with orthogonal components.** In this section we discuss the generalization of our previous findings to rank $r$ spiked tensor model with $r > 1$ of the form

$$ (34) \qquad \mathbf{T} = \sum_{\ell=1}^{r} \beta_{\ell} \mathbf{x}_{\ell}^{(1)} \otimes \dots \otimes \mathbf{x}_{\ell}^{(d)} + \frac{1}{\sqrt{N}} \mathbf{X}, $$

<footnote n="12">¹²corresponding to the minimal theoretical SNR $\beta_s$.</footnote>
<figure><img src="image_20.png" /><figcaption>FIG 12. Minimal signal strength $\beta_s$ and the corresponding singular value and alignments in terms of the tensor order $d$, when all the tensor dimensions are equal.</figcaption></figure>

where $\beta_1 > \dots > \beta_r$ are the signal strengths, $(\mathbf{x}_{\ell}^{(1)}, \dots, \mathbf{x}_{\ell}^{(d)}) \in \mathbb{S}^{n_1-1} \times \dots \times \mathbb{S}^{n_d-1}$ for $\ell \in [r]$, $\mathbf{X} \sim T_{n_1, \dots, n_d}(\mathcal{N}(0, 1))$ and $N = \sum_{i=1}^d n_i$. Supposing that the components $\mathbf{x}_1^{(i)}, \dots, \mathbf{x}_r^{(i)}$ are mutually orthogonal, i.e., $\langle \mathbf{x}_{\ell}^{(i)}, \mathbf{x}_{\ell'}^{(i)} \rangle = 0$ for $\ell \neq \ell'$, the above $r$-rank spiked tensor model can be treated through equivalent rank-one tensors defined for each component $\beta_{\ell} \mathbf{x}_{\ell}^{(1)} \otimes \dots \otimes \mathbf{x}_{\ell}^{(d)}$ independently, by applying the result of the rank-one case established in the previous section. Precisely, the best rank-$r$ approximation of $\mathbf{T}$ corresponds to

$$ (35) \qquad \arg \min_{\lambda_{\ell} \in \mathbb{R}, (\mathbf{u}_{\ell}^{(1)}, \dots, \mathbf{u}_{\ell}^{(d)}) \in \mathbb{S}^{n_1-1} \times \dots \times \mathbb{S}^{n_d-1}} \|\mathbf{T} - \sum_{\ell=1}^{r} \lambda_{\ell} \mathbf{u}_{\ell}^{(1)} \otimes \dots \otimes \mathbf{u}_{\ell}^{(d)}\|_{F}^{2}, $$

and the $\mathbf{u}_{\ell}^{(i)}$'s correlate with $\mathbf{x}_{\ell}^{(i)}$'s as a result of the uniqueness of orthogonal tensor decomposition (see Theorem 4.1 in [1]). Therefore, the study of $\mathbf{T}$ boils down to the study of the random matrices $\Phi_d(\mathbf{T}, \mathbf{u}_{\ell}^{(1)}, \dots, \mathbf{u}_{\ell}^{(d)})$ for $\ell \in [r]$ which behave as the rank-one case treated previously with signal strength $\beta_{\ell}$ respectively. Indeed, since $\mathbf{u}_{\ell}^{(i)} \in \mathbb{S}^{n_i-1}$ by definition and given the orthogonality condition on the $\mathbf{x}_{\ell}^{(i)}$'s, for $\ell' \neq \ell$ the inner product $\langle \mathbf{u}_{\ell}^{(i)}, \mathbf{x}_{\ell'}^{(i)} \rangle_{\text{a.s.}} \to 0$ in high dimension, as such

$$ (36) \qquad \| \Phi_d(\mathbf{T}, \mathbf{u}_{\ell}^{(1)}, \dots, \mathbf{u}_{\ell}^{(d)}) - \Phi_d(\mathbf{T}_{\ell}, \mathbf{u}_{\ell}^{(1)}, \dots, \mathbf{u}_{\ell}^{(d)}) \|_{\text{a.s.}} \to 0, $$

where $\mathbf{T}_{\ell} = \beta_{\ell} \mathbf{x}_{\ell}^{(1)} \otimes \dots \otimes \mathbf{x}_{\ell}^{(d)} + \frac{1}{\sqrt{N}} \mathbf{X}$. Meaning, that the study of $\Phi_d(\mathbf{T}, \mathbf{u}_{\ell}^{(1)}, \dots, \mathbf{u}_{\ell}^{(d)})$ is equivalent to consider the rank-one spiked tensor model $\mathbf{T}_{\ell}$.

**7. Discussion.** In this work, we characterized the asymptotic behavior of spiked asymmetric tensors by mapping them to equivalent (in spectral sense) random matrices. Our starting point is mainly the identities in Eq. (26) which are verified by all critical points of the ML problem. Quite surprisingly and as also discussed in [8] for symmetric tensors, we found that our asymptotic equations describe precisely the maximum of the ML problem which correlates with the true spike. Extrapolating the findings from [12, 8] in the symmetric case, we conjure the existence of an order one threshold $\beta_c$ above which our equations describe the behavior of the global maximum of the ML problem. Unfortunately, it is still unclear how we can characterize such $\beta_c$ with our present approach, which remains an open question. The same question concerns the characterization of the algorithmic threshold $\beta_a$ which is more interesting from a practical standpoint as computing the ML solution is NP-hard for $\beta < \beta_a$.

In the present work, our results were derived under a Gaussian assumption on the tensor noise components. We believe that the derived formulas are universal in the sense that they
extend to other distributions provided that the fourth order moment is finite (as assumed by [3] for long random matrices). Other extensions concern the generalization to higher-ranks with arbitrary components and possibly correlated noise components since the present RMT-tools are more flexible than the use of tools form statistical physics.

**Acknowledgments.** This work was supported by the MIAILargeDATA Chair at University Grenoble Alpes led by R. Couillet and the UGA-HUAWEI LarDist project led by M. Guillaud. We would like to thank Henrique Goulart, Pierre Common and Gérard Ben-Aroud for valuable discussions on the topic of random tensors.

## APPENDIX A: SIMULATIONS

In this section we provide simulations to support our findings.

**A.1. Matrix case.** We start by considering the spiked random matrix model of the form

$$ (37) \qquad M = \beta x y^\top + \frac{1}{\sqrt{m+n}} X, \quad \text{with } X \sim \mathbb{M}_{m,n}(\mathcal{N}(0,1)), $$

and $x, y$ are unitary vectors of dimensions $m$ and $n$ respectively. Figure 13 depicts the spectrum of $[0_{m \times m} M \quad M^\top \quad 0_{n \times n}]$ for $\beta = 0$ for different values of the dimensions $m, n$, and the predicted limiting spectral measure as per Corollary 4. Figure 14 shows a comparison between the asymptotic singular value and alignments obtained in Corollary 5 and their simulated counterparts through SVD applied on $M$. Figure 15 further provides comparison between theory and simulations in the case of long random matrices ($m = 200, n = m^{3/2}$), which corresponds to tensor unfolding as per Remark 9, where a perfect matching is also observed between theory and simulations.

**A.2. Order 3 tensors.** Now we consider a 3-order random tensor model of the form

$$ (38) \qquad T = \beta x \otimes y \otimes z + \frac{1}{\sqrt{m+n+p}} X, \quad \text{with } X \sim \mathbb{T}_{m,n,p}(\mathcal{N}(0,1)), $$

<figure><img src="image_21.png" /><figcaption>FIG 13. Spectrum of $[0_{m \times m} M \quad M^\top \quad 0_{n \times n}]$ for $\beta = 0$ with different values of the dimensions $m, n$ and the limiting spectral measure in Corollary 4.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>FIG 14. Asymptotic singular value and alignments of the spiked random matrix in Eq. (37) as per Corollary 5, and their simulated counterparts for different dimensions *m*, *n*. Simulations are performed through SVD applied to *M*.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>FIG 15. Asymptotic singular value and alignments of the spiked random matrix in Eq. (37) as per Corollary 5, and their simulated counterparts in the long matrix regime ($m = 200, n = \frac{m^2}{3}$) which simulated tensor unfolding as in Remark 9. Simulations are performed through SVD applied to $M$.</figcaption></figure>

and $x, y, z$ are unitary vectors of dimensions $m, n$ and $p$ respectively. In our simulations the estimation of the singular vectors $u, v, w$ of $T$ is performed using the power iteration method described by Algorithm 3. We consider three initialization strategies for Algorithm 3:

(i) Random initialization by randomly sampling $u_0, v_0, w_0$ from the unitary spheres $S^{m-1}, S^{n-1}$ and $S^{p-1}$ respectively. We refer to this strategy in the figures legends by "Random init.".

(ii) Initialization with the true components $x, y, z$. We refer to this strategy in the figures legends by "Init. with $x, y, z$".

(iii) We start by running Algorithm 3 with strategy (i) for a high value of SNR $\beta \gg 1$ then progressively diminishing $\beta$ and initializing Algorithm 3 with the components obtained for the precedent value of $\beta$. We refer to this strategy in the figures legends by "Init. with $u, v, w$".
**Algorithm 3 Tensor power method [1]**

**Input:** An order $d$ tensor $\mathbf{T} \in \mathbb{T}_{n_1, \dots, n_d}$ and initial components $(\mathbf{u}_0^{(1)}, \dots, \mathbf{u}_0^{(d)}) \in \mathbb{S}^{n_1-1} \times \dots \times \mathbb{S}^{n_d-1}$

**Output:** Estimate of arg min$_{\lambda \in \mathbb{R}}$, $(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)}) \in \mathbb{S}^{n_1-1} \times \dots \times \mathbb{S}^{n_d-1} \|\mathbf{T} - \lambda \mathbf{u}^{(1)} \otimes \dots \otimes \mathbf{u}^{(d)}\|_{\mathrm{F}}^2$

**Repeat**

for $i \in [d]$ do

$$ \mathbf{u}^{(i)} \leftarrow \frac{\mathbf{T}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(i-1)}; \mathbf{u}^{(i+1)}, \dots, \mathbf{u}^{(d)})}{\|\mathbf{T}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(i-1)}; \mathbf{u}^{(i+1)}, \dots, \mathbf{u}^{(d)})\|} \quad \triangleright \text{Contract } \mathbf{T} \text{ on all the } \mathbf{u}^{(j)}\text{'s for } j \neq i $$

end for

until convergence of the $\mathbf{u}^{(i)}$'s

Figure 16 provides a comparison between the asymptotic singular value and alignments of $\mathbf{T}$ (obtained by Corollary 3) with their simulation counterparts where the singular vectors are estimated by Algorithm 3 with the initialization strategies (i) in yellow dots and (ii) in green dots. As we can see, as the tensor dimensions grow, the numerical estimates approach their asymptotic counterparts for (ii). Besides, the random initialization (i) yields poor convergence for $\beta$ around its minimal value $\beta_s$ when the tensor dimension is large enough (see $m=n=p=150$). This phenomenon is related to the algorithmic time complexity of the power method which is known to succeed in recovering the underlying signal in polynomial time provided that $\beta \gtrsim n^{\frac{d-2}{4}}$ [17, 6], which is also noticeable in our simulations (the algorithmic phase transition seems to grow with the tensor dimensions). Figure 17 further depicts comparison between theory and simulations when Algorithm 3 is initialized following strategy (iii) which allows to follow the trajectory of the global maximum of the ML problem, where we can notice a good matching between the asymptotic curves and their simulated counterparts.

## APPENDIX B: PROOFS

**B.1. Derivative of tensor singular value and vectors.** Deriving the identities in Eq. (11) and Eq. (12) w.r.t. the entry $X_{ijk}$ of the tensor noise $\mathbf{X}$, we obtain the following set of equations

$$ \left\{ \begin{aligned} & \mathbf{T}(\mathbf{w}) \frac{\partial \mathbf{v}}{\partial X_{ijk}} + \mathbf{T}(\mathbf{v}) \frac{\partial \mathbf{w}}{\partial X_{ijk}} + \frac{1}{\sqrt{N}} v_j w_k e_i^m = \frac{\partial \lambda}{\partial X_{ijk}} \mathbf{u} + \lambda \frac{\partial \mathbf{u}}{\partial X_{ijk}}, \\ & \mathbf{T}(\mathbf{w})^{\top} \frac{\partial \mathbf{u}}{\partial X_{ijk}} + \mathbf{T}(\mathbf{u}) \frac{\partial \mathbf{w}}{\partial X_{ijk}} + \frac{1}{\sqrt{N}} u_i w_k e_j^n = \frac{\partial \lambda}{\partial X_{ijk}} \mathbf{v} + \lambda \frac{\partial \mathbf{v}}{\partial X_{ijk}}, \\ & \mathbf{T}(\mathbf{v})^{\top} \frac{\partial \mathbf{u}}{\partial X_{ijk}} + \mathbf{T}(\mathbf{u})^{\top} \frac{\partial \mathbf{v}}{\partial X_{ijk}} + \frac{1}{\sqrt{N}} u_i v_j e_k^p = \frac{\partial \lambda}{\partial X_{ijk}} \mathbf{w} + \lambda \frac{\partial \mathbf{w}}{\partial X_{ijk}}, \\ & \frac{\partial \lambda}{\partial X_{ijk}} = \mathbf{T} \left( \frac{\partial \mathbf{u}}{\partial X_{ijk}}, \mathbf{v}, \mathbf{w} \right) + \mathbf{T} \left( \mathbf{u}, \frac{\partial \mathbf{v}}{\partial X_{ijk}}, \mathbf{w} \right) + \mathbf{T} \left( \mathbf{u}, \mathbf{v}, \frac{\partial \mathbf{w}}{\partial X_{ijk}} \right) + \frac{1}{\sqrt{N}} u_i v_j w_k. \end{aligned} \right. $$

Writing $\mathbf{T} \left( \frac{\partial \mathbf{u}}{\partial X_{ijk}}, \mathbf{v}, \mathbf{w} \right)$ as $\mathbf{T} \left( \frac{\partial \mathbf{u}}{\partial X_{ijk}}, \mathbf{v}, \mathbf{w} \right) = \left( \frac{\partial \mathbf{u}}{\partial X_{ijk}} \right)^{\top} \mathbf{T}(\mathbf{v})\mathbf{w}$, we can apply again the identities in Eq. (11) which results in $\mathbf{T} \left( \frac{\partial \mathbf{u}}{\partial X_{ijk}}, \mathbf{v}, \mathbf{w} \right) = \lambda \left( \frac{\partial \mathbf{u}}{\partial X_{ijk}} \right)^{\top} \mathbf{u}$. Doing similarly with $\mathbf{T} \left( \mathbf{u}, \frac{\partial \mathbf{v}}{\partial X_{ijk}}, \mathbf{w} \right)$ and $\mathbf{T} \left( \mathbf{u}, \mathbf{v}, \frac{\partial \mathbf{w}}{\partial X_{ijk}} \right)$, we have

$$ \frac{\partial \lambda}{\partial X_{ijk}} = \lambda \left( \left( \frac{\partial \mathbf{u}}{\partial X_{ijk}} \right)^{\top} \mathbf{u} + \left( \frac{\partial \mathbf{v}}{\partial X_{ijk}} \right)^{\top} \mathbf{v} + \left( \frac{\partial \mathbf{w}}{\partial X_{ijk}} \right)^{\top} \mathbf{w} \right) + \frac{1}{\sqrt{N}} u_i v_j w_k. $$

Furthermore, since $\mathbf{u}^{\top} \mathbf{u} = \mathbf{v}^{\top} \mathbf{v} = \mathbf{w}^{\top} \mathbf{w} = 1$, we have

$$ \left( \frac{\partial \mathbf{u}}{\partial X_{ijk}} \right)^{\top} \mathbf{u} = \left( \frac{\partial \mathbf{v}}{\partial X_{ijk}} \right)^{\top} \mathbf{v} = \left( \frac{\partial \mathbf{w}}{\partial X_{ijk}} \right)^{\top} \mathbf{w} = 0. $$
<figure><img src="image_24.png" /><figcaption>FIG 16. Asymptotic singular value and alignments of the order 3 tensor in Eq. (38) having equal dimensions as per Corollary 3, and their simulated counterparts for different dimensions {10, 30, 50, 100, 150}. The yellow dots correspond to a random initialization of the power iteration method strategy (i), while the green dots correspond to an initialization with the true spike components x, y and z strategy (ii).</figcaption></figure>

Thus the derivative of $\lambda$ writes simply as

$$ \frac{\partial \lambda}{\partial X_{ijk}} = \frac{1}{\sqrt{N}} u_i v_j w_k. $$

Hence, we find that

$$ \lambda \begin{bmatrix} \frac{\partial \mathbf{u}}{\partial X_{ijk}} \\ \frac{\partial \mathbf{v}}{\partial X_{ijk}} \\ \frac{\partial \mathbf{w}}{\partial X_{ijk}} \end{bmatrix} = \frac{1}{\sqrt{N}} \begin{bmatrix} v_j w_k (\mathbf{e}_i^m - u_i \mathbf{u}) \\ u_i w_k (\mathbf{e}_j^n - v_j \mathbf{v}) \\ u_i v_j (\mathbf{e}_k^p - w_k \mathbf{w}) \end{bmatrix} + \Phi_3(\mathbf{T}, \mathbf{u}, \mathbf{v}, \mathbf{w}) \begin{bmatrix} \frac{\partial \mathbf{u}}{\partial X_{ijk}} \\ \frac{\partial \mathbf{v}}{\partial X_{ijk}} \\ \frac{\partial \mathbf{w}}{\partial X_{ijk}} \end{bmatrix}. $$
<figure><img src="image_25.png" /><figcaption>FIG 17. Asymptotic singular value and alignments of the order 3 tensor in Eq. (38) having equal dimensions as per Corollary 3, and their simulated counterparts for different dimensions {10, 30, 50, 100, 150}. The green dots correspond the initialization with strategy (iii).</figcaption></figure>

Yielding the expression in Eq. (14). The same calculations apply to the more general d-order tensor case yielding the identity in Eq. (28).

**B.2. Proof of Theorem 3.** Denote the matrix model as

$$N \equiv \frac{1}{\sqrt{N}} \Phi_3(\mathbf{X}, \mathbf{a}, \mathbf{b}, \mathbf{c}),$$
where we recall $\mathbf{X} \sim \mathbb{T}_{m,n,p}(\mathcal{N}(0,1))$ and $(\mathbf{a}, \mathbf{b}, \mathbf{c}) \in \mathbb{S}^{m-1} \times \mathbb{S}^{m-1} \times \mathbb{S}^{p-1}$ are independent of $\mathbf{X}$. We further denote the resolvent matrix of $\mathcal{N}$ as

$$ \mathbf{Q}(z) \equiv (\mathbf{N} - z\mathbf{I}_{\mathcal{N}})^{-1} = \begin{bmatrix} \mathbf{Q}^{11}(z) & \mathbf{Q}^{12}(z) & \mathbf{Q}^{13}(z) \\ \mathbf{Q}^{12}(z)^{\top} & \mathbf{Q}^{22}(z) & \mathbf{Q}^{23}(z) \\ \mathbf{Q}^{13}(z)^{\top} & \mathbf{Q}^{23}(z)^{\top} & \mathbf{Q}^{33}(z) \end{bmatrix}. $$

In order to characterize the limiting Stieltjes transform $g(z)$ of $\mathcal{N}$, we need to estimate the quantity $\frac{1}{N} \mathrm{tr} \mathbf{Q}(z) \xrightarrow{\mathrm{a.s.}} g(z)$ (as a consequence of Theorem 2). We further introduce the following limits

$$ \frac{1}{N} \mathrm{tr} \mathbf{Q}^{11}(z) \xrightarrow{\mathrm{a.s.}} g_1(z), \quad \frac{1}{N} \mathrm{tr} \mathbf{Q}^{22}(z) \xrightarrow{\mathrm{a.s.}} g_2(z), \quad \frac{1}{N} \mathrm{tr} \mathbf{Q}^{33}(z) \xrightarrow{\mathrm{a.s.}} g_3(z). $$

From the identity in Eq. (7), we have

$$ \mathcal{N} \mathbf{Q}(z) - z \mathbf{Q}(z) = \mathbf{I}_{\mathcal{N}}, $$

from which we particularly have

$$ \frac{1}{\sqrt{N}} \left[ \mathbf{X}(\mathbf{c}) \mathbf{Q}^{12}(z)^{\top} \right]_{ii} + \frac{1}{\sqrt{N}} \left[ \mathbf{X}(\mathbf{b}) \mathbf{Q}^{13}(z)^{\top} \right]_{ii} - z \mathbf{Q}_{ii}^{11}(z) = 1, $$

or

$$ (39) \quad \frac{1}{N\sqrt{N}} \sum_{i=1}^{m} \left[ \mathbf{X}(\mathbf{c}) \mathbf{Q}^{12}(z)^{\top} \right]_{ii} + \frac{1}{N\sqrt{N}} \sum_{i=1}^{m} \left[ \mathbf{X}(\mathbf{b}) \mathbf{Q}^{13}(z)^{\top} \right]_{ii} - \frac{z}{N} \mathrm{tr} \mathbf{Q}^{11}(z) = \frac{m}{N}. $$

We thus need to compute the expectations of $\frac{1}{N\sqrt{N}} \sum_{i=1}^{m} \left[ \mathbf{X}(\mathbf{c}) \mathbf{Q}^{12}(z)^{\top} \right]_{ii}$ and $\frac{1}{N\sqrt{N}} \sum_{i=1}^{m} \left[ \mathbf{X}(\mathbf{b}) \mathbf{Q}^{13}(z)^{\top} \right]_{ii}$, which develop as

$$ \frac{1}{N\sqrt{N}} \sum_{i=1}^{m} \mathbb{E} \left[ \mathbf{X}(\mathbf{c}) \mathbf{Q}^{12}(z)^{\top} \right]_{ii} = \frac{1}{N\sqrt{N}} \sum_{i=1}^{m} \sum_{j=1}^{n} \sum_{k=1}^{p} c_k \mathbb{E} \left[ X_{ijk} Q_{ij}^{12} \right] = \frac{1}{N\sqrt{N}} \sum_{i=1}^{m} \sum_{j=1}^{n} \sum_{k=1}^{p} c_k \mathbb{E} \left[ \frac{\partial Q_{ij}^{12}}{\partial X_{ijk}} \right], $$

where the last equality is obtained by applying Stein's lemma (Lemma 1). For continuing the derivations, we need to express the derivative of the resolvent $\mathbf{Q}(z)$ with respect to an entry $X_{ijk}$ of the tensor noise $\mathbf{X}$. Indeed, since $\mathcal{N}\mathbf{Q}(z) - z\mathbf{Q}(z) = \mathbf{I}_{\mathcal{N}}$, we have

$$ \frac{\partial \mathcal{N}}{\partial X_{ijk}} \mathbf{Q}(z) + \mathcal{N} \frac{\partial \mathbf{Q}(z)}{\partial X_{ijk}} - z \frac{\partial \mathbf{Q}(z)}{\partial X_{ijk}} = \mathbf{0}_{N \times N} \quad \Rightarrow \quad (\mathcal{N} - z \mathbf{I}_{\mathcal{N}}) \frac{\partial \mathbf{Q}(z)}{\partial X_{ijk}} = - \frac{\partial \mathcal{N}}{\partial X_{ijk}} \mathbf{Q}(z), $$

from which we get

$$ \frac{\partial \mathbf{Q}(z)}{\partial X_{ijk}} = -\mathbf{Q}(z) \frac{\partial \mathcal{N}}{\partial X_{ijk}} \mathbf{Q}(z), $$

where

$$ \frac{\partial \mathcal{N}}{\partial X_{ijk}} = \frac{1}{\sqrt{N}} \begin{bmatrix} \mathbf{0}_{m \times m} & c_k \mathbf{e}_i^m (\mathbf{e}_j^n)^{\top} & b_j \mathbf{e}_i^m (\mathbf{e}_k^p)^{\top} \\ c_k \mathbf{e}_j^m (\mathbf{e}_i^m)^{\top} & \mathbf{0}_{n \times n} & a_i \mathbf{e}_j^n (\mathbf{e}_k^p)^{\top} \\ b_j \mathbf{e}_k^p (\mathbf{e}_i^m)^{\top} & a_i \mathbf{e}_k^p (\mathbf{e}_j^n)^{\top} & \mathbf{0}_{p \times p} \end{bmatrix}, $$

and we finally obtain the following derivatives

$$ \frac{\partial Q_{ab}^{11}}{\partial X_{ijk}} = -\frac{1}{\sqrt{N}} [a_i(Q_{aj}^{12}Q_{bk}^{13} + Q_{ak}^{13}Q_{bj}^{12}) + b_j(Q_{ai}^{11}Q_{bk}^{13} + Q_{ak}^{13}Q_{ib}^{11}) + c_k(Q_{ai}^{11}Q_{bj}^{12} + Q_{aj}^{12}Q_{ib}^{11})], $$

$$ \frac{\partial Q_{ab}^{12}}{\partial X_{ijk}} = -\frac{1}{\sqrt{N}} [a_i(Q_{aj}^{12}Q_{bk}^{23} + Q_{ak}^{13}Q_{jb}^{22}) + b_j(Q_{ai}^{11}Q_{bk}^{23} + Q_{ak}^{13}Q_{ib}^{12}) + c_k(Q_{ai}^{11}Q_{jb}^{22} + Q_{aj}^{12}Q_{ib}^{11})], $$
$$ \frac{\partial Q_{ab}^{13}}{\partial X_{ijk}} = -\frac{1}{\sqrt{N}} [a_i(Q_{aj}^{12}Q_{kb}^{33} + Q_{ak}^{13}Q_{jb}^{23}) + b_j(Q_{ai}^{11}Q_{kb}^{33} + Q_{ak}^{13}Q_{ib}^{13}) + c_k(Q_{ai}^{11}Q_{jb}^{23} + Q_{aj}^{12}Q_{ib}^{13})], $$

$$ \frac{\partial Q_{ab}^{22}}{\partial X_{ijk}} = -\frac{1}{\sqrt{N}} [a_i(Q_{aj}^{22}Q_{bk}^{23} + Q_{ak}^{23}Q_{jb}^{22}) + b_j(Q_{ia}^{12}Q_{bk}^{23} + Q_{ak}^{23}Q_{ib}^{12}) + c_k(Q_{ia}^{12}Q_{jb}^{22} + Q_{aj}^{22}Q_{ib}^{12})], $$

$$ \frac{\partial Q_{ab}^{23}}{\partial X_{ijk}} = -\frac{1}{\sqrt{N}} [a_i(Q_{aj}^{22}Q_{kb}^{33} + Q_{ak}^{23}Q_{jb}^{23}) + b_j(Q_{ia}^{12}Q_{kb}^{33} + Q_{ak}^{23}Q_{ib}^{13}) + c_k(Q_{ia}^{12}Q_{jb}^{23} + Q_{aj}^{22}Q_{ib}^{13})], $$

$$ \frac{\partial Q_{ab}^{33}}{\partial X_{ijk}} = -\frac{1}{\sqrt{N}} [a_i(Q_{ja}^{23}Q_{kb}^{33} + Q_{ak}^{33}Q_{jb}^{23}) + b_j(Q_{ia}^{13}Q_{kb}^{33} + Q_{ak}^{33}Q_{ib}^{13}) + c_k(Q_{ia}^{13}Q_{jb}^{23} + Q_{ja}^{23}Q_{ib}^{13})]. $$

In particular,

$$ \frac{\partial Q_{ij}^{12}}{\partial X_{ijk}} = -\frac{1}{\sqrt{N}} [a_i(Q_{ij}^{12}Q_{jk}^{23} + Q_{ik}^{13}Q_{jj}^{22}) + b_j(Q_{ii}^{11}Q_{jk}^{23} + Q_{ik}^{13}Q_{ij}^{12}) + c_k(Q_{ii}^{11}Q_{jj}^{22} + Q_{ij}^{12}Q_{ij}^{12})]. $$

Going back to the computation of $A \equiv \frac{1}{N\sqrt{N}} \sum_{i=1}^{m} \sum_{j=1}^{n} \sum_{k=1}^{p} c_k \mathbb{E} \left[ \frac{\partial Q_{ij}^{12}}{\partial X_{ijk}} \right]$, we will see that the only contributing term in the derivative $\frac{\partial Q_{ij}^{12}}{\partial X_{ijk}}$ is $-\frac{1}{\sqrt{N}} c_k Q_{ii}^{11} Q_{jj}^{22}$. Indeed,

$$ \begin{aligned} A &= -\frac{1}{N^2} \sum_{ijk} c_k \mathbb{E} [a_i(Q_{ij}^{12}Q_{jk}^{23} + Q_{ik}^{13}Q_{jj}^{22}) + b_j(Q_{ii}^{11}Q_{jk}^{23} + Q_{ik}^{13}Q_{ij}^{12}) + c_k(Q_{ii}^{11}Q_{jj}^{22} + Q_{ij}^{12}Q_{ij}^{12})], \\ &= -\frac{1}{N^2} \mathbb{E} [a^\top Q^{12} Q^{23} c + a^\top Q^{13} c \operatorname{tr} Q^{22} + \operatorname{tr} Q^{11} b^\top Q^{23} c + b^\top (Q^{12})^\top Q^{13} c + \operatorname{tr} Q^{11} \operatorname{tr} Q^{22} + \operatorname{tr}(Q^{12}(Q^{12})^\top)] \end{aligned} $$

Now, since the vectors $a, b, c$ are of bounded norms and assuming $Q(z)$ is of bounded spectral norm (see condition in Eq. (8)), under Assumption 1 (as $N \to \infty$), the terms $\frac{1}{N^2} a^\top Q^{12} Q^{23} c$, $\frac{1}{N^2} a^\top Q^{13} c \operatorname{tr} Q^{22}$, $\frac{1}{N^2} \operatorname{tr} Q^{11} b^\top Q^{23} c$, $\frac{1}{N^2} b^\top (Q^{12})^\top Q^{13} c$ and $\frac{1}{N^2} \operatorname{tr}(Q^{12}(Q^{12})^\top)$ are vanishing almost surely. As such, we find that

$$ \frac{1}{N\sqrt{N}} \sum_{i=1}^{m} \left[ \mathbf{x}(c) Q^{12}(z)^{\top} \right]_{ii} = -\frac{1}{N} \operatorname{tr} Q^{11}(z) \frac{1}{N} \operatorname{tr} Q^{22}(z) + \mathcal{O}(N^{-1}) \xrightarrow{\text{a.s.}} -g_1(z)g_2(z) + \mathcal{O}(N^{-1}). $$

Similarly, we find that

$$ \frac{1}{N\sqrt{N}} \sum_{i=1}^{m} \left[ \mathbf{x}(b) Q^{13}(z)^{\top} \right]_{ii} = -\frac{1}{N} \operatorname{tr} Q^{11}(z) \frac{1}{N} \operatorname{tr} Q^{33}(z) + \mathcal{O}(N^{-1}) \xrightarrow{\text{a.s.}} -g_1(z)g_3(z) + \mathcal{O}(N^{-1}). $$

From Eq. (39), $g_1(z)$ satisfies

$$ -g_1(z)(g_2(z) + g_3(z)) - zg_1(z) = c_1, $$

where we recall $c_1 = \lim_{N} \frac{m}{N}$. Similarly, $g_2(z)$ and $g_3(z)$ satisfy

$$ \begin{cases} -g_2(z)(g_1(z) + g_3(z)) - zg_2(z) = c_2, \\ -g_3(z)(g_1(z) + g_2(z)) - zg_3(z) = c_3, \end{cases} $$

where we recall again $c_2 = \lim_{N} \frac{n}{N}$ and $c_3 = \lim_{N} \frac{p}{N}$. Moreover, by definition, $g(z) = \sum_{i=1}^{3} g_i(z)$, thus we have for each $i \in [3]$

$$ g_i(z)(g(z) - g_i(z)) + zg_i(z) + c_i = 0, $$

yielding

$$ g_i(z) = \frac{g(z) + z}{2} - \frac{\sqrt{4c_i + (g(z) + z)^2}}{2}, $$

with $g(z)$ solution of the equation $g(z) = \sum_{i=1}^{3} g_i(z)$ satisfying $\mathfrak{F}[g(z)] > 0$ for $z \in \mathbb{C}$ with $\mathfrak{F}[z] > 0$ (see Property 2 of the Stieltjes transform in Subsection 2.2).
**B.3. Proof of Corollary 1.** Given the result of Theorem 3, setting $c_1 = c_2 = c_3 = \frac{1}{3}$, we have for all $i \in [3]$

$$g_i(z) = \frac{g(z) + z}{2} - \frac{\sqrt{\frac{4}{3} + (g(z) + z)^2}}{2},$$

where $g(z)$ satisfies $g(z) = \sum_{i=1}^3 g_i(z)$, thus $g(z)$ is the solution to

$$z + \frac{g(z) + z}{2} - \frac{3\sqrt{\frac{4}{3} + (g(z) + z)^2}}{2} = 0,$$

solving in $g(z)$ yields

$$g(z) \in \left\{ -\frac{3z}{4} - \frac{\sqrt{3\sqrt{3z^2 - 8}}}{4}, -\frac{3z}{4} + \frac{\sqrt{3\sqrt{3z^2 - 8}}}{4} \right\},$$

and the limiting Stieltjes transform (with $\Im[g(z)] > 0$ for $z$ with $\Im(z) > 0$ by property 2 of the Stieltjes transform, see Subsection 2.2) is therefore

$$g(z) = -\frac{3z}{4} + \frac{\sqrt{3\sqrt{3z^2 - 8}}}{4}.$$

**B.4. Proof of Theorem 4.** Given the random tensor model in Eq. (10) and its singular vectors characterized by Eq. (11), we denote the associated random matrix model as

$$T \equiv \Phi_3(\mathbf{T}, \mathbf{u}, \mathbf{v}, \mathbf{w}) = \beta V B V^\top + N,$$

where

$$N = \frac{1}{\sqrt{N}} \Phi_3(\mathbf{X}, \mathbf{u}, \mathbf{v}, \mathbf{w}), \quad B \equiv \begin{bmatrix} 0 & \langle z, w \rangle \langle y, v \rangle \\ \langle z, w \rangle & 0 \\ \langle y, v \rangle & \langle x, u \rangle \end{bmatrix} \in \mathbb{M}_3, \quad V \equiv \begin{bmatrix} \mathbf{x} & \mathbf{0}_m & \mathbf{0}_m \\ \mathbf{0}_n & \mathbf{y} & \mathbf{0}_n \\ \mathbf{0}_p & \mathbf{0}_p & \mathbf{z} \end{bmatrix} \in \mathbb{M}_{N,3}.$$

We further denote the resolvents of $T$ and $N$ respectively as

$$R(z) = (T - zI_N)^{-1} = \begin{bmatrix} R^{11}(z) & R^{12}(z) & R^{13}(z) \\ R^{12}(z)^\top & R^{22}(z) & R^{23}(z) \\ R^{13}(z)^\top & R^{23}(z)^\top & R^{33}(z) \end{bmatrix},$$

$$Q(z) = (N - zI_N)^{-1} = \begin{bmatrix} Q^{11}(z) & Q^{12}(z) & Q^{13}(z) \\ Q^{12}(z)^\top & Q^{22}(z) & Q^{23}(z) \\ Q^{13}(z)^\top & Q^{23}(z)^\top & Q^{33}(z) \end{bmatrix}.$$

By Woodbury matrix identity (Lemma 5), we have

$$(40) \qquad R(z) = Q(z) - Q(z)V \left( \frac{1}{\beta} B^{-1} + V^\top Q(z)V \right)^{-1} V^\top Q(z),$$

In particular, taking the normalized trace operator, we get

$$\frac{1}{N} \operatorname{tr} R(z) = \frac{1}{N} \operatorname{tr} Q(z) - \frac{1}{N} \operatorname{tr} \left[ \left( \frac{1}{\beta} B^{-1} + V^\top Q(z) V \right)^{-1} V^\top Q^2(z) V \right], \\ = \frac{1}{N} \operatorname{tr} Q(z) + \mathcal{O}(N^{-1}),$$

since the matrix $(\frac{1}{\beta}B^{-1} + V^\top Q(z)V)^{-1}V^\top Q^2(z)V$ is of bounded spectral norm (assuming $\|Q(z)\|$ is bounded, see condition in Eq. (8)) and being of finite size (3 × 3 matrix).
As such, the asymptotic spectral measure of $T$ is the same as the one of $N$ which can be estimated though $\frac{1}{N} \operatorname{tr} Q(z)$. Comparing to the result from Appendix B.2, now the singular vectors $u, v, w$ depend statistically on the tensor noise $X$ which needs to be handled.

From the identity Eq. (7), we have $NQ(z) - zQ(z) = I_N$, from which

$$ (41) \quad \frac{1}{N\sqrt{N}} \sum_{i=1}^{m} \left[ \mathbf{X}(\mathbf{w})\mathbf{Q}^{12}(z)^{\top} \right]_{ii} + \frac{1}{N\sqrt{N}} \sum_{i=1}^{m} \left[ \mathbf{X}(\mathbf{v})\mathbf{Q}^{13}(z)^{\top} \right]_{ii} - \frac{z}{N} \operatorname{tr} \mathbf{Q}^{11}(z) = \frac{m}{N}. $$

We thus need to compute the expectations of $\frac{1}{N\sqrt{N}} \sum_{i=1}^{m} [\mathbf{X}(\mathbf{w})\mathbf{Q}^{12}(z)^{\top}]_{ii}$ and $\frac{1}{N\sqrt{N}} \sum_{i=1}^{m} [\mathbf{X}(\mathbf{v})\mathbf{Q}^{13}(z)^{\top}]_{ii}$. In particular,

$$ A \equiv \frac{1}{N\sqrt{N}} \sum_{i=1}^{m} \mathbb{E} \left[ \mathbf{X}(\mathbf{w}) \mathbf{Q}^{12}(z)^{\top} \right]_{ii} = \frac{1}{N\sqrt{N}} \sum_{ijk} \mathbb{E} \left[ X_{ijk} w_k Q_{ij}^{12} \right] = \frac{1}{N\sqrt{N}} \sum_{ijk} \mathbb{E} \left[ \frac{\partial(w_k Q_{ij}^{12})}{\partial X_{ijk}} \right], $$

where the last equality is obtained by Stein's lemma (Lemma 1). Due to the statistical dependency between $w$ and $X$, the above sum decomposes in two terms which are

$$ A = \frac{1}{N\sqrt{N}} \sum_{ijk} \mathbb{E} \left[ w_k \frac{\partial Q_{ij}^{12}}{\partial X_{ijk}} \right] + \frac{1}{N\sqrt{N}} \sum_{ijk} \mathbb{E} \left[ \frac{\partial w_k}{\partial X_{ijk}} Q_{ij}^{12} \right] = A_1 + A_2, $$

where the first term $A_1$ has already been handled in the previous subsection (if replacing $c$ with $w$). We now show that the second term $A_2$ is asymptotically vanishing under Assumption 1. Indeed, by Eq. (14), we have

$$ \frac{\partial w_k}{\partial X_{ijk}} = -\frac{1}{\sqrt{N}} \left( v_j w_k (R_{ik}^{13}(\lambda) - u_i \mathbf{u}^{\top} \mathbf{R}_{:,k}^{13}(\lambda)) \right) \\ - \frac{1}{\sqrt{N}} \left( u_i w_k (R_{jk}^{23}(\lambda) - v_j \mathbf{v}^{\top} \mathbf{R}_{:,k}^{23}(\lambda)) \right) - \frac{1}{\sqrt{N}} \left( u_i v_j (R_{kk}^{33}(\lambda) - w_k \mathbf{w}^{\top} \mathbf{R}_{:,k}^{33}(\lambda)) \right). $$

As such $A_2$ decomposes in three terms $A_2 = A_{21} + A_{22} + A_{23}$, where

$$ \begin{aligned} A_{21} &= -\frac{1}{N^2} \sum_{ijk} \mathbb{E} [v_j w_k R_{ik}^{13}(\lambda) Q_{ij}^{12}] + \frac{1}{N^2} \sum_{ijkl} \mathbb{E} [u_i v_j w_k u_l R_{lk}^{13}(\lambda) Q_{ij}^{12}], \\ &= -\frac{1}{N^2} \mathbb{E} [\mathbf{v}^{\top} \mathbf{Q}^{12}(z)^{\top} \mathbf{R}^{13}(\lambda) \mathbf{w}] + \frac{1}{N^2} \mathbb{E} [\mathbf{u}^{\top} \mathbf{Q}^{12}(z) \mathbf{v} \mathbf{u}^{\top} \mathbf{R}^{13}(\lambda) \mathbf{w}] \to 0, \end{aligned} $$

as $N \to \infty$, since the singular vectors $u, v, w$ are of bounded norms and assuming the resolution $Q(z)$ and $R(\lambda)$ are of bounded spectral norms ($R(\lambda)$ has bounded spectral norm by Assumption 2 and Eq. (8)). Similarly, we further have

$$ \begin{aligned} A_{22} &= -\frac{1}{N^2} \sum_{ijk} \mathbb{E} [u_i w_k R_{jk}^{23}(\lambda) Q_{ij}^{12}] + \frac{1}{N^2} \sum_{ijkl} \mathbb{E} [u_i v_j w_k v_l R_{lk}^{23}(\lambda) Q_{ij}^{12}], \\ &= -\frac{1}{N^2} \mathbb{E} [\mathbf{u}^{\top} \mathbf{Q}^{12}(z) \mathbf{R}^{23}(\lambda) \mathbf{w}] + \frac{1}{N^2} \mathbb{E} [\mathbf{u}^{\top} \mathbf{Q}^{12}(z) \mathbf{v} \mathbf{v}^{\top} \mathbf{R}^{23}(\lambda) \mathbf{w}] \to 0. \end{aligned} $$

And finally,

$$ \begin{aligned} A_{23} &= -\frac{1}{N^2} \sum_{ijk} \mathbb{E} [u_i v_j R_{kk}^{33}(\lambda) Q_{ij}^{12}] + \frac{1}{N^2} \sum_{ijkl} \mathbb{E} [u_i v_j w_k w_l R_{lk}^{33}(\lambda) Q_{ij}^{12}], \\ &= -\frac{1}{N} \mathbb{E} \left[ \mathbf{u}^{\top} \mathbf{Q}^{12}(z) \mathbf{v} \frac{1}{N} \operatorname{tr} \mathbf{R}^{33}(\lambda) \right] + \frac{1}{N^2} \mathbb{E} \left[ \mathbf{u}^{\top} \mathbf{Q}^{12}(z) \mathbf{v} \mathbf{w}^{\top} \mathbf{R}^{33}(\lambda) \mathbf{w} \right] \to 0. \end{aligned} $$
Therefore,

$$A = \frac{1}{N\sqrt{N}} \sum_{ijk} \mathbb{E} \left[ w_k \frac{\partial Q_{ij}^{12}}{\partial X_{ijk}} \right] + \mathcal{O}(N^{-1}).$$

As in the previous subsection, the derivative of $Q(z)$ w.r.t. the entry $X_{ijk}$ expresses as
$$ \frac{\partial Q(z)}{\partial X_{ijk}} = -Q(z) \frac{\partial N}{\partial X_{ijk}} Q(z) \text{ but now with} $$

$$ \frac{\partial N}{\partial X_{ijk}} = \frac{1}{\sqrt{N}} \left[ \begin{array}{ccc} \mathbf{0}_{m \times m} & w_k e_i^m (e_j^n)^\top & v_j e_i^m (e_k^p)^\top \\ w_k e_j^n (e_i^m)^\top & \mathbf{0}_{n \times n} & u_i e_j^p (e_k^p)^\top \\ v_j e_k^p (e_i^m)^\top & u_i e_k^p (e_j^n)^\top & \mathbf{0}_{p \times p} \end{array} \right] + \frac{1}{\sqrt{N}} \Phi_3 \left( \mathbf{x}, \frac{\partial \mathbf{u}}{\partial X_{ijk}}, \frac{\partial \mathbf{v}}{\partial X_{ijk}}, \frac{\partial \mathbf{w}}{\partial X_{ijk}} \right), $$

where $O \equiv \frac{1}{\sqrt{N}} \Phi_3(\mathbf{x}, \frac{\partial \mathbf{u}}{\partial X_{ijk}}, \frac{\partial \mathbf{v}}{\partial X_{ijk}}, \frac{\partial \mathbf{w}}{\partial X_{ijk}})$ is of vanishing spectral norm. Indeed, from Eq. (14), there exists $C > 0$ independent of $N$ such that $\|\frac{\partial \mathbf{u}}{\partial X_{ijk}}\|, \|\frac{\partial \mathbf{v}}{\partial X_{ijk}}\|, \|\frac{\partial \mathbf{w}}{\partial X_{ijk}}\| \le \frac{C}{\sqrt{N}}$ yielding that the spectral norm of $O$ is bounded by $\frac{C'}{N^{3/2}}$ for some constant $C' > 0$ independent of $N$. Therefore, we find that $A \to -g_1(z)g_2(z) + \mathcal{O}(N^{-1})$ (with $g_1(z), g_2(z)$ the almost sure limits of $\frac{1}{N} \text{tr} \mathbf{Q}^{11}(z)$ and $\frac{1}{N} \text{tr} \mathbf{Q}^{22}(z)$ respectively), thus yielding the same limiting Stieltjes transform as the one obtained in Appendix B.2.

**B.5. Proof of Theorem 5.** Given the identities in Eq. (11), we have

$$ (42) \qquad \lambda \langle \mathbf{u}, \mathbf{x} \rangle = \mathbf{x}^{\top} \mathbf{T}(\mathbf{v}) \mathbf{w} = \beta \langle \mathbf{v}, \mathbf{y} \rangle \langle \mathbf{w}, \mathbf{z} \rangle + \frac{1}{\sqrt{N}} \mathbf{x}^{\top} \mathbf{X}(\mathbf{v}) \mathbf{w}, $$

with $\lambda, \langle \mathbf{u}, \mathbf{x} \rangle, \langle \mathbf{v}, \mathbf{y} \rangle$ and $\langle \mathbf{w}, \mathbf{z} \rangle$ converging almost surely to their asymptotic limits $\lambda^\infty(\beta), a_x^\infty(\beta), a_y^\infty(\beta)$ and $a_z^\infty(\beta)$ respectively given the concentration properties in Subsection 3.4. To characterize such limits we need to evaluate the expectation of $\frac{1}{\sqrt{N}} \mathbf{x}^\top \mathbf{X}(\mathbf{v}) \mathbf{w}$. Indeed,

$$ \frac{1}{\sqrt{N}} \mathbb{E} \left[ \mathbf{x}^{\top} \mathbf{X}(\mathbf{v}) \mathbf{w} \right] = \frac{1}{\sqrt{N}} \sum_{ijk} x_i \mathbb{E} \left[ v_j w_k X_{ijk} \right] = \frac{1}{\sqrt{N}} \sum_{ijk} x_i \mathbb{E} \left[ \frac{\partial v_j}{\partial X_{ijk}} w_k \right] + \frac{1}{\sqrt{N}} \sum_{ijk} x_i \mathbb{E} \left[ v_j \frac{\partial w_k}{\partial X_{ijk}} \right] $$

where the last equality is obtained by Stein's lemma (Lemma 1). From Eq. (14), we have

$$ \frac{\partial v_j}{\partial X_{ijk}} = -\frac{1}{\sqrt{N}} \left( v_j w_k (R_{ij}^{12}(\lambda) - u_i \mathbf{u}^{\top} \mathbf{R}_{:,j}^{12}(\lambda)) \right) \\ -\frac{1}{\sqrt{N}} \left( u_i w_k (R_{jj}^{22}(\lambda) - v_j \mathbf{v}^{\top} \mathbf{R}_{:,j}^{22}(\lambda)) \right) - \frac{1}{\sqrt{N}} \left( u_i v_j (R_{jk}^{23}(\lambda) - w_k \mathbf{w}^{\top} \mathbf{R}_{:,j}^{32}(\lambda)) \right). $$

Hence $A \equiv \frac{1}{\sqrt{N}} \sum_{ijk} x_i \mathbb{E} \left[ \frac{\partial v_j}{\partial X_{ijk}} w_k \right] = A_1 + A_2 + A_3$ decomposes in three terms $A_1, A_2$ and $A_3$. The terms $A_1$ and $A_3$ will be vanishing asymptotically and only $A_2$ contains non-vanishing terms. Indeed,

$$ \begin{aligned} A_1 &= -\frac{1}{N} \sum_{ijk} \mathbb{E} \left[ x_i w_k v_j w_k R_{ij}^{12}(\lambda) \right] + \frac{1}{N} \sum_{ijkl} \mathbb{E} \left[ x_i w_k v_j w_k u_i u_l R_{lj}^{12}(\lambda) \right], \\ &= -\frac{1}{N} \mathbb{E} \left[ \mathbf{x}^{\top} \mathbf{R}^{12}(\lambda) \mathbf{v} \right] + \frac{1}{N} \mathbb{E} \left[ \langle \mathbf{x}, \mathbf{u} \rangle \mathbf{u}^{\top} \mathbf{R}^{12}(\lambda) \mathbf{v} \right] \to 0, \end{aligned} $$

as $N \to \infty$ since $\mathbf{x}, \mathbf{u}$ and $\mathbf{v}$ are of bounded norms and $\mathbf{R}(\lambda)$ being of bounded spectral norm for $\lambda$ outside the support of $\frac{1}{\sqrt{N}} \Phi_3(\mathbf{x}, \mathbf{u}, \mathbf{v}, \mathbf{w})$ (through the identity in Eq. (40)). Similarly
with $A_3$, we have

$$
\begin{align*}
A_3 &= -\frac{1}{N} \sum_{ijk} \mathbb{E} [x_i w_k u_i v_j R_{jk}^{23}(\lambda)] + \frac{1}{N} \sum_{ijkl} \mathbb{E} [x_i w_k u_i v_j w_k w_l R_{jl}^{23}(\lambda)], \\
&= -\frac{1}{N} \mathbb{E} [\langle \mathbf{x}, \mathbf{u} \rangle \mathbf{v}^\top \mathbf{R}^{23}(\lambda) \mathbf{w}] + \frac{1}{N} \mathbb{E} [\langle \mathbf{x}, \mathbf{u} \rangle \mathbf{v}^\top \mathbf{R}^{23}(\lambda) \mathbf{w}] \to 0.
\end{align*}
$$

Now $A_2$ is not vanishing, precisely,

$$
\begin{align*}
A_2 &= -\frac{1}{N} \sum_{ijk} \mathbb{E} [x_i w_k u_i w_k R_{jj}^{22}(\lambda)] + \frac{1}{N} \sum_{ijkl} \mathbb{E} [x_i w_k u_i w_k v_j v_l R_{jl}^{22}(\lambda)], \\
&= -\frac{1}{N} \mathbb{E} [\langle \mathbf{x}, \mathbf{u} \rangle \operatorname{tr} \mathbf{R}^{22}(\lambda)] + \frac{1}{N} \mathbb{E} [\langle \mathbf{x}, \mathbf{u} \rangle \mathbf{v}^\top \mathbf{R}^{22}(\lambda) \mathbf{v}], \\
&\to -g_2(\lambda) \mathbb{E} [\langle \mathbf{x}, \mathbf{u} \rangle] + \mathcal{O}(N^{-1}),
\end{align*}
$$

where the last line results from the fact that $\frac{1}{N} \operatorname{tr} \mathbf{R}^{22}(\lambda) \xrightarrow{\text{a.s.}} g_2(\lambda)$ as we saw in the previous subsection. Similarly, we find that

$$
B \equiv \frac{1}{\sqrt{N}} \sum_{ijk} x_i \mathbb{E} \left[ v_j \frac{\partial w_k}{\partial X_{ijk}} \right] \to -g_3(\lambda) \mathbb{E} \left[ \langle \mathbf{x}, \mathbf{u} \rangle \right] + \mathcal{O}(N^{-1}).
$$

Therefore, by Eq. (42), the almost sure limits $\lambda^\infty(\beta)$, $a_x^\infty(\beta)$, $a_y^\infty(\beta)$ and $a_z^\infty(\beta)$ satisfy the equation

$$
\lambda^{\infty}(\beta) a_{\mathbf{x}}^{\infty}(\beta) = \beta a_{\mathbf{y}}^{\infty}(\beta) a_{\mathbf{z}}^{\infty}(\beta) - [g_2(\lambda^{\infty}(\beta)) + g_3(\lambda^{\infty}(\beta))] a_{\mathbf{x}}^{\infty}(\beta).
$$

Hence,

$$
a_x^\infty(\beta) = \alpha_1(\lambda^\infty(\beta))a_y^\infty(\beta)a_z^\infty(\beta),
$$

with $\alpha_1(z) \equiv \frac{\beta}{z+g(z)-g_1(z)}$. Similarly, we find that

$$
\begin{cases}
a_{\mathbf{y}}^{\infty}(\beta) = \alpha_2(\lambda^{\infty}(\beta))a_{\mathbf{x}}^{\infty}(\beta)a_{\mathbf{z}}^{\infty}(\beta), \\
a_{\mathbf{z}}^{\infty}(\beta) = \alpha_3(\lambda^{\infty}(\beta))a_{\mathbf{x}}^{\infty}(\beta)a_{\mathbf{y}}^{\infty}(\beta),
\end{cases}
$$

with $\alpha_i(z) \equiv \frac{\beta}{z+g(z)-g_i(z)}$. Solving the above system of equations provides the final asymptotic alignments. Proceeding similarly with Eq. (12) we obtain an estimate of the asymptotic singular value $\lambda^\infty(\beta)$, thereby ending the proof.

**B.6. Proof of Corollary 3.** By Corollary 1, the limiting Stieltjes transform is given by

$$
g(z) = -\frac{3z}{4} + \frac{\sqrt{3\sqrt{3z^2 - 8}}}{4},
$$

and since $g(z) = \sum_{i=1}^{3} g_i(z)$ with all $g_i(z)$ equal, then $g_i(z) = \frac{g(z)}{3}$ for all $i \in [3]$. Hence, for each $i \in [3]$, $\alpha_i(z)$ defined in Theorem 5 writes as

$$
\alpha_i(z) = \frac{\beta}{\frac{z}{2} + \frac{\sqrt{3\sqrt{3z^2-8}}}{6}},
$$

and $\lambda^\infty(\beta)$ is solution to the equation

$$
\frac{z}{4} + \frac{\sqrt{3\sqrt{3z^2 - 8}}}{4} - \frac{\left(\frac{z}{2} + \frac{\sqrt{3\sqrt{3z^2 - 8}}}{6}\right)^3}{\beta^2} = 0.
$$
First we compute the critical value of $\beta$ by solving the above equation in $\beta$ and taking the limit when $z$ tends to the right edge of the semi-circle law (i.e., $z \to 2\sqrt{\frac{2}{3}}$). Indeed, solving the above equation in $\beta$ yields

$$ \beta(z) = \sqrt{\frac{2z^3 + \frac{2z^2\sqrt{9z^2-24}}{3} - 4z - \frac{4\sqrt{9z^2-24}}{9}}{z + \sqrt{3}\sqrt{3z^2-8}}, $$

hence

$$ \beta_s = \lim_{z \to 2\sqrt{\frac{2}{3}}} \beta(z) = \frac{2\sqrt{3}}{3}. $$

Now to express $\lambda^\infty(\beta)$ in terms of $\beta$, we solve the equation $\beta - \beta(z) = 0$ in $z$ and choose the unique non-decreasing (in $\beta$) and positive solution, which yields

$$ \lambda^\infty(\beta) = \sqrt{\frac{\beta^2}{2} + 2 + \frac{\sqrt{3}\sqrt{(3\beta^2 - 4)^3}}{18\beta}}. $$

Plugging the above expression of $\lambda^\infty(\beta)$ into the expressions of the asymptotic alignments in Theorem 5, we obtain for all $i \in [3]$

$$ \alpha_i(\lambda^\infty(\beta)) = \frac{6\sqrt{2}\beta}{\sqrt{9\beta^2 - 12 + \frac{\sqrt{3}\sqrt{(3\beta^2 - 4)^3}}{\beta}} + \sqrt{9\beta^2 + 36 + \frac{\sqrt{3}\sqrt{(3\beta^2 - 4)^3}}{\beta}}}, $$

yielding the final result.

**B.7. Proof of Corollary 4.** Setting $c_1 = c$, $c_2 = 1-c$ and $c_3 = 0$, we get

$$ \left\{ \begin{array}{l} g_1(z) = \frac{g(z)+z}{2} - \frac{\sqrt{4c+(g(z)+z)^2}}{2}, \\ g_2(z) = \frac{g(z)+z}{2} - \frac{\sqrt{4(1-c)+(g(z)+z)^2}}{2}, \\ g_3(z) = 0. \end{array} \right. $$

And since $g(z) = \sum_{i=1}^3 g_i(z)$, then $g(z)$ satisfies the equation

$$ z - \frac{\sqrt{4c + (g(z) + z)^2}}{2} - \frac{\sqrt{4(1-c) + (g(z) + z)^2}}{2} = 0, $$

the solution of which belongs to

$$ g(z) \in \left\{ -z - \frac{\sqrt{4c(c-1) + (z^2-1)^2}}{z}, -z + \frac{\sqrt{4c(c-1) + (z^2-1)^2}}{z} \right\}, $$

thus the limiting Stieltjes transform (having $\Im[g(z)] > 0$) is given by

$$ g(z) = -z + \frac{\sqrt{4c(c-1) + (z^2-1)^2}}{z}. $$

In particular, the edges of the support of the corresponding limiting distribution $\nu$ are the roots of $4c(c-1) + (z^2-1)^2$, yielding

$$ S(\nu) = \left[ -\sqrt{1 + 2\sqrt{c(1-c)}}, -\sqrt{1 - 2\sqrt{c(1-c)}} \right] \cup \left[ \sqrt{1 - 2\sqrt{c(1-c)}}, \sqrt{1 + 2\sqrt{c(1-c)}} \right]. $$
And the density function of $\nu$ is obtained by computing the limit $\lim_{\epsilon \to 0} |\mathfrak{S}[g(x + i\epsilon)]|$ yielding

$$ \nu(dx) = \frac{1}{\pi x} \sin\left(\frac{\arctan_2(0, q_c(x))}{2}\right) \sqrt{|q_c(x)|} \operatorname{sign}\left(\frac{\sin\left(\frac{\arctan_2(0, q_c(x))}{2}\right)}{x}\right) dx + (1 - 2\min(c, 1-c)) \delta(x), $$

where $q_c(x) = (x^2 - 1)^2 + 4c(c - 1)$. The Dirac $\delta(x)$ component in the above expression corresponds to the fact that the corresponding matrix model is of rank $2\min(m, n)$.

**B.8. Proof of Corollary 5.** From Subsection B.7, plugging the expression of the limiting Stieltjes transform $g(z)$ into the expressions defining $g_1(z)$ and $g_2(z)$ yields

$$ \left\{ \begin{aligned} g_1(z) &= \frac{\sqrt{4c(c-1)+(z^2-1)^2}}{2z} - \frac{\sqrt{4c(c-1+z^2)+(z^2-1)^2}}{2z}, \\ g_2(z) &= \frac{\sqrt{4c(c-1)+(z^2-1)^2}}{2z} - \frac{\sqrt{4c(c-1-z^2)+(z^2+1)^2}}{2z}. \end{aligned} \right. $$

Thus, by Theorem 5, we have

$$ \left\{ \begin{aligned} \alpha_1(z) &= \frac{\beta}{\frac{\sqrt{4c^2-4c+z^4-2z^2+1}}{2z} + \frac{\sqrt{4c^2+4cz^2-4c+z^4-2z^2+1}}{2z}}, \\ \alpha_2(z) &= \frac{\beta}{\frac{\sqrt{4c^2-4c+z^4-2z^2+1}}{2z} + \frac{\sqrt{4c^2-4cz^2-4c+z^4+2z^2+1}}{2z}}, \\ \alpha_3(z) &= \frac{\beta z}{\sqrt{4c^2-4c+z^4-2z^2+1}}. \end{aligned} \right. $$

Then solving the equation $z + g(z) - \frac{\beta}{\alpha_1(z)\alpha_2(z)\alpha_3(z)}$ in $z$ provides the almost sure limit of $\lambda$ in terms of $\beta$. Specifically,

$$ \lambda^{\infty}(\beta) = \sqrt{\beta^2 + 1 + \frac{c(1-c)}{\beta^2}}. $$

Plugging the above expression of $\lambda^\infty(\beta)$ into $\frac{1}{\sqrt{\alpha_2(z)\alpha_3(z)}}$ by replacing $z$ with the expression of $\lambda^\infty(\beta)$ provides the asymptotic alignment $a_x^\infty(\beta)$ as $\frac{1}{\kappa(\beta, c)}$, with $\kappa(\beta, c) = \sqrt{\alpha_2(\lambda^\infty(\beta))\alpha_3(\lambda^\infty(\beta))}$ given by the following expression

$$ \kappa(\beta, c) = \sqrt{ \frac{2\sqrt{\frac{\beta^2(\beta^2+1)-c^2+c}{\beta^2}(\beta^2(\beta^2+1)-c(c-1))} \left( \frac{\sqrt{4\beta^4c(c-1)+(\beta^4-c^2+c)^2}}{\beta^4} - \frac{\sqrt{\frac{\beta^2(\beta^2+1)-c(c-1)}{\beta^2}(\beta^2(\beta^2+1)-c(c-1))}}{\beta^4} \right) + \frac{\sqrt{\frac{\beta^2c(\beta^2(\beta^2-c^2)-c(c-1))+(\beta^2(\beta^2+2)-c(c-1))^2}{\beta^4}} + \frac{\sqrt{\frac{\beta^4c(c-1)+(\beta^4-c(c-1))^2}{\beta^4}}}{\beta^4} \right) } $$

From the identities

$$ \begin{aligned} & 4\beta^4 c(c-1) + (\beta^4 - c(c-1))^2 = (\beta^4 + c(c-1))^2 \\ & - 4\beta^2 c(\beta^2(\beta^2 - c + 2) - c(c-1)) + (\beta^2(\beta^2 + 2) - c(c-1))^2 = (\beta^4 + \beta^2(2-2c) - c^2 + c)^2 \\ \kappa(\beta, c) & \text{simplifies as} \end{aligned} $$

$$ \kappa(\beta, c) = \beta \sqrt{\frac{\beta^2 (\beta^2 + 1) - c(c - 1)}{(\beta^4 + c(c - 1)) (\beta^2 + 1 - c)}}. $$

The asymptotic alignment $a_y^\infty(\beta)$ is given by $\frac{1}{\kappa(\beta, 1-c)}$ since the dimension ratio of the component $y$ is $c_2 = 1 - c$. Moreover, the critical value of $\beta$ is obtained by solving the equation $\lambda^\infty(\beta_s) = \sqrt{1 + 2\sqrt{c(1-c)}}$, i.e., when the limiting singular value gets closer to
the right edge of the support of the corresponding limiting spectral measure (see Corollary 4). Finally, similarly as above, we can check that $\sqrt{\alpha_1(\lambda^\infty(\beta))\alpha_2(\lambda^\infty(\beta))}$ is equal to 1 for $\beta \ge \beta_s$ (i.e., the alignment along the third dimension is 1, corresponding to the component $z$ in Eq. (10)).

**B.9. Proof of Theorem 6.** Denote the matrix model as

$$N \equiv \frac{1}{\sqrt{N}} \Phi_d(\mathbf{X}, \mathbf{a}^{(1)}, \dots, \mathbf{a}^{(d)}),$$

with $\mathbf{X} \sim \mathbb{T}_{n_1, \dots, n_d}(\mathcal{N}(0, 1))$ and $(\mathbf{a}^{(1)}, \dots, \mathbf{a}^{(d)}) \in \mathbb{S}^{n_1-1} \times \dots \times \mathbb{S}^{n_d-1}$ are independent of $\mathbf{X}$. We further denote the resolvent of $N$ as

$$Q(z) = (N - zI_N)^{-1} = \begin{bmatrix} Q^{11}(z) & Q^{12}(z) & Q^{13}(z) & \dots & Q^{1d}(z) \\ Q^{12}(z)^\top & Q^{22}(z) & Q^{23}(z) & \dots & Q^{2d}(z) \\ Q^{13}(z)^\top & Q^{23}(z)^\top & Q^{33}(z) & \dots & Q^{3d}(z) \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ Q^{1d}(z)^\top & Q^{2d}(z)^\top & Q^{3d}(z)^\top & \dots & Q^{dd}(z) \end{bmatrix}.$$

By Borel-Cantelli lemma, we have $\frac{1}{N} \text{tr} \, Q(z) \xrightarrow{\text{a.s.}} g(z)$ and for all $i \in [d]$, $\frac{1}{N} \text{tr} \, Q^{ii}(z) \xrightarrow{\text{a.s.}} g_i(z)$. Applying the identity in Eq. (7) to the symmetric matrix $N$ we get $NQ(z) - zQ(z) = I_N$, from which we particularly get

$$\frac{1}{\sqrt{N}} \sum_{j=2}^{d} \left[ \mathbf{x}^{1j} \mathbf{Q}^{1j}(z)^{\top} \right]_{i_1 i_1} - z \mathbf{Q}^{11}_{i_1 i_1}(z) = 1,$$

or

$$(43) \qquad \frac{1}{N\sqrt{N}} \sum_{j=1}^{2} \sum_{i_1=1}^{n_1} \left[ \mathbf{x}^{1j} \mathbf{Q}^{1j}(z)^{\top} \right]_{i_1 i_1} - \frac{z}{N} \text{tr} \, \mathbf{Q}^{11}(z) = \frac{n_1}{N},$$

where we recall that $\mathbf{x}^{ij} \equiv \mathbf{X}(\mathbf{a}^{(1)}, \dots, \mathbf{a}^{(i-1)}, \dots, \mathbf{a}^{(i+1)}, \dots, \mathbf{a}^{(j-1)}, \dots, \mathbf{a}^{(j+1)}, \dots, \mathbf{a}^{(d)}) \in \mathbb{M}_{n_1, n_d}$.

We thus need to compute the expectation of $\frac{1}{N\sqrt{N}} \sum_{i_1=1}^{n_1} \left[ \mathbf{x}^{1j} \mathbf{Q}^{1j}(z)^{\top} \right]_{i_1 i_1}$ which develops as

$$A_j \equiv \frac{1}{N\sqrt{N}} \sum_{i_1=1}^{n_1} \mathbb{E} \left[ \mathbf{x}^{1j} \mathbf{Q}^{1j}(z)^{\top} \right]_{i_1 i_1} = \frac{1}{N\sqrt{N}} \sum_{i_1 i_2} \mathbb{E} \left[ \left[ \mathbf{x}^{1j} \right]_{i_1 i_2} \mathbf{Q}^{1j}_{i_1 i_2} \right] = \frac{1}{N\sqrt{N}} \sum_{i_1, \dots, i_d} \prod_{k \neq 1, k \neq j} a_{i_k}^{(k)} \mathbb{E} \left[ \frac{\partial \mathbf{Q}^{1j}_{i_1 i_j}}{\partial X_{i_1, \dots, i_d}} \right]$$

where the last equality follows from Stein's lemma (Lemma 1). In particular, as in Appendix B.2 for the 3-order case, it turns out that the only contributing term in the derivative $\frac{\partial \mathbf{Q}^{1j}_{i_1 i_j}}{\partial X_{i_1, \dots, i_d}}$ is $-\frac{1}{\sqrt{N}} \prod_{k \neq 1, k \neq j} a_{i_k}^{(k)} \mathbf{Q}^{11}_{i_1 i_1} \mathbf{Q}^{jj}_{i_j i_j}$ with the other terms yielding quantities of order $\mathcal{O}(N^{-1})$. Therefore, we find that

$$\begin{align*}
A_j &= \frac{1}{N\sqrt{N}} \sum_{i_1, \dots, i_d} \prod_{k \neq 1, k \neq j} a_{i_k}^{(k)} \mathbb{E} \left[ \frac{\partial \mathbf{Q}^{1j}_{i_1 i_j}}{\partial X_{i_1, \dots, i_d}} \right] = -\frac{1}{N^2} \sum_{i_1, \dots, i_d} \prod_{k \neq 1, k \neq j} \left( a_{i_k}^{(k)} \right)^2 \mathbb{E} \left[ \mathbf{Q}^{11}_{i_1 i_1} \mathbf{Q}^{jj}_{i_j i_j} \right] + \mathcal{O}(N^{-1}) \\
&= -\frac{1}{N^2} \sum_{i_1 i_j} \mathbb{E} \left[ \mathbf{Q}^{11}_{i_1 i_1} \mathbf{Q}^{jj}_{i_j i_j} \right] + \mathcal{O}(N^{-1}) \\
&= -\frac{1}{N} \text{tr} \, \mathbf{Q}^{11}(z) \frac{1}{N} \text{tr} \, \mathbf{Q}^{jj}(z) + \mathcal{O}(N^{-1}) \\
&\xrightarrow{\text{a.s.}} -g_1(z)g_j(z) + \mathcal{O}(N^{-1}).
\end{align*}$$
From Eq. (43), $g_1(z)$ satisfies $-g_1(z) \sum_{j \neq i} g_j(z) - zg_1(z) = c_1$ with $c_1 = \lim_{N \to \infty} \frac{n_1}{N}$. Similarly, for all $i \in [d]$, $g_i(z)$ satisfies $-g_i(z) \sum_{j \neq i} g_j(z) - zg_i(z) = c_i$ with $c_i = \lim_{N \to \infty} \frac{n_i}{N}$. And since, $g(z) = \sum_{i=1}^d g_i(z)$, we have for each $i \in [d]$, $g_i(z)(g(z) - g_i(z)) + zg_i(z) + c_i = 0$, yielding

$$g_i(z) = \frac{g(z) + z}{2} - \frac{\sqrt{4c_i + (g(z) + z)^2}}{2},$$

with $g(z)$ solution to the equation $g(z) = \sum_{i=1}^d g_i(z)$ satisfying $\mathfrak{I}[g(z)] > 0$ for $\mathfrak{I}[z] > 0$.

**B.10. Proof of Corollary 6.** Given Theorem 6 and setting $c_i = \frac{1}{d}$ for all $i \in [d]$, we have $g_i(z) = \frac{g(z)}{d}$, thus $g(z)$ satisfies the equation

$$\frac{g(z)}{d} = \frac{g(z) + z}{2} - \frac{\sqrt{\frac{4}{d} + (g(z) + z)^2}}{2}.$$

Solving in $g(z)$ yields

$$g(z) \in \left\{ -\frac{dz}{2(d-1)} - \frac{\sqrt{d(dz^2 - 4d + 4)}}{2(d-1)}, -\frac{dz}{2(d-1)} + \frac{\sqrt{d(dz^2 - 4d + 4)}}{2(d-1)} \right\},$$

and the limiting Stieltjes transform with $\mathfrak{I}[g(z)] \ge 0$ is

$$g(z) = -\frac{dz}{2(d-1)} + \frac{\sqrt{d(dz^2 - 4d + 4)}}{2(d-1)}.$$

**B.11. Existence of $g(z)$.** Define the following function for $(g, z) \in \mathbb{R} \times \mathbb{R}_+$ and $d \ge 3$

$$k(g, z) = g - \frac{g + z}{2}d + \frac{1}{2} \sum_{i=1}^{d} \sqrt{4c_i + (g + z)^2},$$

with $\sum_{i=1}^d c_i = 1$. By concavity of the function $\sqrt{\cdot}$ we have $\sum_{i=1}^d \sqrt{4c_i + x^2} \le d\sqrt{\frac{4}{d} + x^2}$ for all $x \in \mathbb{R}$. Therefore, $k(g, z)$ is bounded as

$$k(g, z) \le \bar{k}(g, z) \equiv g - \frac{g + z}{2}d + \frac{d}{2}\sqrt{\frac{4}{d} + (g + z)^2}.$$

From Section B.10, we have, for any $z > 2\sqrt{\frac{d-1}{d}}$, there exists $g_* \in \mathbb{R}$ such that $\bar{k}(g_*, z) = 0$. Besides, we further have $\lim_{g \to -\infty} k(g, z) = \lim_{g \to +\infty} k(g, z) = +\infty$ and hence by continuity, there exists $g(z)$ such that $k(g(z), z) = 0$ for $z$ large enough (e.g., $z > 2\sqrt{\frac{d-1}{d}}$).

**B.12. Proof of Theorem 7.** Given the random tensor model in Eq. (1) and its singular vectors characterized by Eq. (26), we denote the associated random matrix model as

$$T \equiv \Phi_d (\mathbf{T}, \mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)}) = \beta V B V^\top + N,$$

where $N = \frac{1}{\sqrt{N}} \Phi_d (\mathbf{X}, \mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)})$, $\mathbf{B} \in \mathbb{M}_d$ with entries $B_{ij} = (1 - \delta_{ij}) \prod_{k \neq i, j} \langle \mathbf{x}^{(k)}, \mathbf{u}^{(k)} \rangle$

$$V \equiv \begin{bmatrix} \mathbf{x}^{(1)} & \mathbf{0}_{n_1} & \cdots & \mathbf{0}_{n_1} \\ \mathbf{0}_{n_2} & \mathbf{x}^{(2)} & \cdots & \mathbf{0}_{n_2} \\ \vdots & \vdots & \ddots & \mathbf{0}_{n_3} \\ \mathbf{0}_{n_d} & \mathbf{0}_{n_d} & \mathbf{0}_{n_d} & \mathbf{x}^{(d)} \end{bmatrix} \in \mathbb{M}_{N,d}.$$
We further denote the resolvent of $T$ and $N$ respectively as

$$R(z) = (T - zI_N)^{-1} = \begin{bmatrix} R^{11}(z) & R^{12}(z) & R^{13}(z) & \cdots & R^{1d}(z) \\ R^{12}(z)^\top & R^{22}(z) & R^{23}(z) & \cdots & R^{2d}(z) \\ R^{13}(z)^\top & R^{23}(z)^\top & R^{33}(z) & \cdots & R^{3d}(z) \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ R^{1d}(z)^\top & R^{2d}(z)^\top & R^{3d}(z)^\top & \cdots & R^{dd}(z) \end{bmatrix}.$$

$$Q(z) = (N - zI_N)^{-1} = \begin{bmatrix} Q^{11}(z) & Q^{12}(z) & Q^{13}(z) & \cdots & Q^{1d}(z) \\ Q^{12}(z)^\top & Q^{22}(z) & Q^{23}(z) & \cdots & Q^{2d}(z) \\ Q^{13}(z)^\top & Q^{23}(z)^\top & Q^{33}(z) & \cdots & Q^{3d}(z) \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ Q^{1d}(z)^\top & Q^{2d}(z)^\top & Q^{3d}(z)^\top & \cdots & Q^{dd}(z) \end{bmatrix}.$$

Similarly as in 3-order case, by Woodbury matrix identity (Lemma 5), we have

$$\frac{1}{N} \operatorname{tr} \mathbf{R}(z) = \frac{1}{N} \operatorname{tr} \mathbf{Q}(z) - \frac{1}{N} \operatorname{tr} \left[ \left( \frac{1}{\beta} \mathbf{B}^{-1} + \mathbf{V}^{\top} \mathbf{Q}(z) \mathbf{V} \right)^{-1} \mathbf{V}^{\top} \mathbf{Q}^2(z) \mathbf{V} \right], \\ = \frac{1}{N} \operatorname{tr} \mathbf{Q}(z) + \mathcal{O}(N^{-1}),$$

since the perturbation matrix $(\frac{1}{\beta}\mathbf{B}^{-1} + \mathbf{V}^{\top}\mathbf{Q}(z)\mathbf{V})^{-1}\mathbf{V}^{\top}\mathbf{Q}^2(z)\mathbf{V}$ is of bounded spectral norm (if $\|\mathbf{Q}(z)\|$ is bounded, see condition in Eq. (8)) and has finite size ($d \times d$ matrix). Therefore, the characterization of the spectrum of $T$ boils down to the estimation of $\frac{1}{N} \operatorname{tr} \mathbf{Q}(z)$. Now we are left to handle the statistical dependency between the tensor noise $\mathbf{X}$ and the singular vectors of $T$. Recalling the proof of Appendix B.9, we have again

$$\frac{1}{N\sqrt{N}} \sum_{j=1}^{2} \sum_{i_1=1}^{n_1} \left[ \mathbf{X}^{1j} \mathbf{Q}^{1j}(z)^{\top} \right]_{i_1 i_1} - \frac{z}{N} \operatorname{tr} \mathbf{Q}^{11}(z) = \frac{n_1}{N},$$

with $\mathbf{X}^{ij} = \mathbf{X}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(i-1)}, \mathbf{u}^{(i+1)}, \dots, \mathbf{u}^{(j-1)}, \mathbf{u}^{(j+1)}, \dots, \mathbf{u}^{(d)}) \in \mathbb{M}_{n_i, n_j}$. Taking the expectation of $\frac{1}{N\sqrt{N}} \sum_{i_1=1}^{n_1} [\mathbf{X}^{1j} \mathbf{Q}^{1j}(z)^{\top}]_{i_1 i_1}$ yields

$$A_j = \frac{1}{N\sqrt{N}} \sum_{i_1, \dots, i_d} \mathbb{E} \left[ \frac{\partial Q_{i_1 i_j}^{1j}}{\partial X_{i_1, \dots, i_d}} \prod_{k \neq 1, k \neq j} u_{i_k}^{(k)} \right] + \frac{1}{N\sqrt{N}} \sum_{i_1, \dots, i_d} \mathbb{E} \left[ Q_{i_1 i_j}^{1j} \prod_{k \neq 1, k \neq j} \frac{\partial u_{i_k}^{(k)}}{\partial X_{i_1, \dots, i_d}} \right] = A_{j1} + A_{j2},$$

where we already computed the first term ($A_{j1}$) in Appendix B.9. Now we will show that $A_{j2}$ is asymptotically vanishing under Assumption 3. Indeed, by Eq. (28), the higher order terms arise from the term $-\frac{1}{\sqrt{N}} \prod_{\ell \neq k} u_{i_\ell}^{(\ell)} R_{i_k i_k}^{kk}(\lambda)$, we thus only show that the contribution of this term is also vanishing. Precisely,

$$\begin{align*}
A_{j2} &= -\frac{1}{N^2} \sum_{i_1, \dots, i_d} \mathbb{E} \left[ Q_{i_1 i_j}^{1j} \prod_{k \neq 1, k \neq j} \prod_{\ell \neq k} u_{i_\ell}^{(\ell)} R_{i_k i_k}^{kk}(\lambda) \right] + \mathcal{O}(N^{-1}) \\
&= -\frac{1}{N^2} \sum_{i_1, \dots, i_d} \mathbb{E} \left[ \left( u_{i_1}^{(1)} \right)^{d-2} Q_{i_1 i_j}^{1j} \left( u_{i_j}^{(j)} \right)^{d-2} \prod_{k \neq 1, k \neq j} u_{i_k}^{(k)} R_{i_k i_k}^{kk}(\lambda) \right] + \mathcal{O}(N^{-1}) \\
&= -\frac{1}{N^2} \mathbb{E} \left[ \left( (\mathbf{u}^{(1)})^{\odot d-2} \right)^{\top} \mathbf{Q}^{1j}(z) \left( \mathbf{u}^{(j)} \right)^{\odot d-2} \sum_{i_2, \dots, i_{j-1}, i_{j+1}, \dots, i_d} \prod_{k \neq 1, k \neq j} u_{i_k}^{(k)} R_{i_k i_k}^{kk}(\lambda) \right] + \mathcal{O}(N^{-1})
\end{align*}$$
where $\mathbf{a}^{\odot q}$ denotes the vector with entries $a_i^q$. As such, $A_{j2}$ is vanishing asymptotically since the singular vectors $\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)}$ are unitary and since their entries are bounded by 1.

We finally need to check that the derivative of $Q(z)$ w.r.t. the entry $X_{i_1, \dots, i_d}$ has the same expression asymptotically as the one in the independent case of Appendix B.9. Indeed, we have $\frac{\partial Q(z)}{\partial X_{i_1, \dots, i_d}} = -Q(z) \frac{\partial N}{\partial X_{i_1, \dots, i_d}} Q(z)$ with

$$ \frac{\partial N}{\partial X_{i_1, \dots, i_d}} = \begin{bmatrix} \mathbf{0}_{n_1 \times n_1} & \mathbf{C}_{12} & \cdots & \mathbf{C}_{1d} \\ \mathbf{C}_{12}^{\top} & \mathbf{0}_{n_2 \times n_2} & \cdots & \mathbf{C}_{2d} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{C}_{1d}^{\top} & \mathbf{C}_{2d}^{\top} & \cdots & \mathbf{0}_{n_d \times n_d} \end{bmatrix} + \frac{1}{\sqrt{N}} \Phi_d \left( \mathbf{X}, \frac{\partial \mathbf{u}^{(1)}}{\partial X_{i_1, \dots, i_d}}, \dots, \frac{\partial \mathbf{u}^{(d)}}{\partial X_{i_1, \dots, i_d}} \right), $$

where $C_{ij} = \prod_{k \neq i, k \neq j} u_{ik}^{(k)} e_{i_i}^{n_i} (e_{ij}^{n_j})^{\top}$ and $O = \frac{1}{\sqrt{N}} \Phi_d (\mathbf{X}, \frac{\partial \mathbf{u}^{(1)}}{\partial X_{i_1, \dots, i_d}}, \dots, \frac{\partial \mathbf{u}^{(d)}}{\partial X_{i_1, \dots, i_d}})$ is of vanishing norm. Indeed, by Eq. (28), there exists $C > 0$ independent of $N$ such that $\|\frac{\partial \mathbf{u}^{(1)}}{\partial X_{i_1, \dots, i_d}}\|, \dots, \|\frac{\partial \mathbf{u}^{(d)}}{\partial X_{i_1, \dots, i_d}}\| \le \frac{C}{\sqrt{N}}$, therefore, the spectral norm of $O$ is bounded by $\frac{C'}{N^{d/2}}$ for some constant $C' > 0$ independent of $N$. Finally, $A_j \to -g_1(z)g_j(z) + O(N^{-1})$ (with $g_1(z), g_j(z)$ the almost sure limits of $\frac{1}{N} \operatorname{tr} Q^{11}(z)$ and $\frac{1}{N} \operatorname{tr} Q^{jj}(z)$ respectively), hence yielding the same limiting Stieltjes transform as the one obtained in Appendix B.9.

**B.13. Proof of Theorem 8.** Given the identities in Eq. (26), we have for all $i \in [d]$

$$ \frac{1}{\sqrt{N}} \sum_{j_1, \dots, j_d} x_{j_i}^{(i)} \prod_{k \neq i} u_{j_k}^{(k)} X_{j_1, \dots, j_d} + \beta \prod_{k \neq i} \langle \mathbf{u}^{(k)}, \mathbf{x}^{(k)} \rangle = \lambda \langle \mathbf{u}^{(i)}, \mathbf{x}^{(i)} \rangle, $$

with $\lambda$ and $\langle \mathbf{u}^{(i)}, \mathbf{x}^{(i)} \rangle$ concentrate almost surely around their asymptotic denoted $\lambda^{\infty}(\beta)$ and $a_{\mathbf{x}^{(i)}}^{\infty}(\beta)$ respectively. Taking the expectation of the first term and applying Stein's lemma (Lemma 1), we get

$$ A_i = \frac{1}{\sqrt{N}} \sum_{j_1, \dots, j_d} x_{j_i}^{(i)} \mathbb{E} \left[ \frac{\partial \left( \prod_{k \neq i} u_{j_k}^{(k)} \right)}{\partial X_{j_1, \dots, j_d}} \right] = \frac{1}{\sqrt{N}} \sum_{j_1, \dots, j_d} x_{j_i}^{(i)} \sum_{k \neq i} \mathbb{E} \left[ \frac{\partial u_{j_k}^{(k)}}{\partial X_{j_1, \dots, j_d}} \prod_{\ell \neq k, \ell \neq i} u_{j_\ell}^{(\ell)} \right], $$

where the only contributing term in the expression of $\frac{\partial u_{j_k}^{(k)}}{\partial X_{j_1, \dots, j_d}}$ from Eq. (28) is $-\frac{1}{\sqrt{N}} \prod_{\ell \neq k} u_{j_\ell}^{(\ell)} R_{j_k j_k}^{kk}(\lambda)$ which yields

$$ \begin{aligned} A_i &= -\frac{1}{N} \sum_{j_1, \dots, j_d} x_{j_i}^{(i)} \sum_{k \neq i} \mathbb{E} \left[ R_{j_k j_k}^{kk}(\lambda) \prod_{\ell \neq k} u_{j_\ell}^{(\ell)} \prod_{\ell \neq k, \ell \neq i} u_{j_\ell}^{(\ell)} \right] + \mathcal{O}(N^{-1}) \\ &= -\mathbb{E} \left[ \langle \mathbf{u}^{(i)}, \mathbf{x}^{(i)} \rangle \sum_{k \neq i} \frac{1}{N} \operatorname{tr} \mathbf{R}^{kk}(\lambda) \right] + \mathcal{O}(N^{-1}) \\ &\to -\mathbb{E} \left[ \langle \mathbf{u}^{(i)}, \mathbf{x}^{(i)} \rangle \right] \sum_{k \neq i} g_k(\lambda). \end{aligned} $$

Therefore, the almost sure limits $\lambda^{\infty}$ and $a_{\mathbf{x}^{(i)}}^{\infty}$ for each $i \in [d]$ satisfy

$$ \lambda^{\infty} a_{\mathbf{x}^{(i)}}^{\infty} = \beta \prod_{k \neq i} a_{\mathbf{x}^{(k)}}^{\infty} - a_{\mathbf{x}^{(i)}}^{\infty} \sum_{k \neq i} g_k(\lambda^{\infty}), $$

therefore

$$ a_{\mathbf{x}^{(i)}}^{\infty} = \alpha_i(\lambda^{\infty}) \prod_{k \neq i} a_{\mathbf{x}^{(k)}}^{\infty}, \quad \text{with} \quad \alpha_i(z) = \frac{\beta}{z + g(z) - g_i(z)}, $$
since $g(z) = \sum_{k=1}^{d} g_i(z)$. To solve the above equation, we simply write $x_i = a_{x^{(i)}}^\infty$ and $\alpha_i = \alpha_i(z)$ by omitting the dependence on $z$. We therefore have

$$x_i = \alpha_i \prod_{k \neq i} x_k \quad \Rightarrow \quad x_i = \alpha_i x_j \prod_{k \neq i, k \neq j} x_k = \alpha_i x_j \prod_{k \neq i, k \neq j} \left( \alpha_k \prod_{\ell \neq k} x_\ell \right),$$

from which we have

$$x_i = x_j \left( \prod_{k \neq j} \alpha_k \right) \left( \prod_{k \neq i, k \neq j} \prod_{\ell \neq k} x_\ell \right) = x_j \left( \prod_{k \neq j} \alpha_k \right) \left( \prod_{k \neq i, k \neq j} \prod_{\ell \neq k, \ell \neq i} x_\ell \right) x_i^{d-2},$$

thus

$$x_j \left( \prod_{k \neq j} \alpha_k \right) \left( \prod_{k \neq i, k \neq j} \prod_{\ell \neq k, \ell \neq i} x_\ell \right) x_i^{d-3} = 1,$$

and we remark that $(\prod_{k \neq i, k \neq j} \prod_{\ell \neq k, \ell \neq i} x_\ell) x_i^{d-3} = (\frac{x_j}{\alpha_j})^{d-3} x_j^{d-2}$, hence $x_j$ is given by

$$x_j = \left( \frac{\alpha_j^{d-3}}{\prod_{k \neq j} \alpha_k} \right)^{\frac{1}{2d-4}},$$

which ends the proof.

Alternative expression of $q_i(z)$. From the above, we have $\lambda^\infty + g(\lambda^\infty) = \beta \prod_{i=1}^d x_i$ and

$$(\lambda^\infty + g(\lambda^\infty) - g_i(\lambda^\infty))x_i = \beta \prod_{j \neq i} x_j \quad \Rightarrow \quad (\lambda^\infty + g(\lambda^\infty) - g_i(\lambda^\infty))x_i^2 = \beta \prod_{i=1}^d x_i.$$

Therefore,

$$x_i = \sqrt{\frac{\lambda^\infty + g(\lambda^\infty)}{\lambda^\infty + g(\lambda^\infty) - g_i(\lambda^\infty)}} = \sqrt{1 + \frac{g_i(\lambda^\infty)}{\lambda^\infty + g(\lambda^\infty) - g_i(\lambda^\infty)}},$$

with $z + g(z) - g_i(z) = \frac{-c_i}{g_i(z)}$ since $g_i(z)$ satisfies $g_i^2(z) - (g(z) + z)g_i(z) - c_i = 0$. Hence, we find

$$x_i = \sqrt{1 - \frac{g_i^2(\lambda^\infty)}{c_i}}.$$

### B.14. Additional lemmas.

LEMMA 4 (Spectral norm of random Gaussian tensors). Let $\mathbf{X} \sim T_{n_1, \dots, n_d}(\mathcal{N}(0, 1))$, then the spectral norm of $\|\mathbf{X}\|$ can be bounded, with probability at least $1 - \delta$ for $\delta > 0$, as

$$\|\mathbf{X}\| \le \sqrt{2 \left[ \left( \sum_{i=1}^{d} n_i \right) \log \left( \frac{2d}{\log(3/2)} \right) + \log \left( \frac{2}{\delta} \right) \right]}.$$
PROOF. By definition, the spectral norm of $\mathbf{X}$ is given as

$$ (44) \qquad \|\mathbf{X}\| = \sup_{\mathbf{u}^{(i)} \in \mathbb{S}^{n_i-1}, i \in [d]} \mathbf{X}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)}) $$

For some $\varepsilon > 0$, let $\mathcal{E}_1, \dots, \mathcal{E}_d$ be $\varepsilon$-nets of $\mathbb{S}^{n_1-1}, \dots, \mathbb{S}^{n_d-1}$ respectively. Since $\mathbb{S}^{n_1-1} \times \dots \times \mathbb{S}^{n_d-1}$ is compact, there exists a maximizer of Eq. (44). And with the $\varepsilon$-net argument, there exists $\mathbf{e}^{(i)} \in \mathcal{E}_i$ for each $i \in [d]$ such that

$$ \|\mathbf{X}\| = \mathbf{X}(\mathbf{e}^{(1)} + \mathbf{\delta}^{(1)}, \dots, \mathbf{e}^{(d)} + \mathbf{\delta}^{(d)}), $$

such that $\|\mathbf{\delta}^{(i)}\| \le \varepsilon$ for $i \in [d]$. Therefore, one has

$$ \|\mathbf{X}\| \le \mathbf{X}(\mathbf{e}^{(1)}, \dots, \mathbf{e}^{(d)}) + \left( \sum_{i=1}^{d} \varepsilon^i \binom{d}{i} \right) \|\mathbf{X}\|. $$

For $\varepsilon = \frac{\log(3/2)}{d}$, one has $\sum_{i=1}^{d} \varepsilon^i \binom{d}{i} \le \sum_{i=1}^{d} \frac{(\varepsilon d)^i}{i!} \le e^{\varepsilon d} - 1 = \frac{1}{2}$. As such, the spectral norm of $\mathbf{X}$ can be bounded as

$$ (45) \qquad \|\mathbf{X}\| \le 2 \max_{\mathbf{e}^{(i)} \in \mathcal{E}_i, i \in [d]} \mathbf{X}(\mathbf{e}^{(1)}, \dots, \mathbf{e}^{(d)}). $$

Since the entries of $\mathbf{X}$ are i.i.d. standard Gaussian random variables, we have $\mathbb{E}[e^{t\mathbf{X}_{i_1 \dots i_d}}] \le e^{t^2/2}$, hence using Hoeffding's inequality we have for any $\mathbf{u}^{(i)} \in \mathbb{S}^{n_i-1}$ with $i \in [d]$

$$ \begin{aligned} \mathbb{P}\{\mathbf{X}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)}) \ge t\} &\le \mathbb{P}\{e^{s\mathbf{X}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)})} \ge e^{st}\} \le e^{-st} \mathbb{E}[e^{s\mathbf{X}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)})}] \\ &\le \exp\left\{-st + \frac{s^2}{2} \sum_{i_1, \dots, i_d} (u_{i_1}^{(1)})^2 \dots (u_{i_d}^{(d)})^2\right\} = \exp\left\{-st + \frac{s^2}{2}\right\}. \end{aligned} $$

Minimizing the right-hand side w.r.t. $s$ yields $\mathbb{P}\{\mathbf{X}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)}) \ge t\} \le e^{-t^2/2}$. With the same arguments, we further have $\mathbb{P}\{\mathbf{X}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)}) \le -t\} \le e^{-t^2/2}$. And taking the union of the two cases yields

$$ \mathbb{P}\{|\mathbf{X}(\mathbf{u}^{(1)}, \dots, \mathbf{u}^{(d)})| \ge t\} \le 2e^{-t^2/2}. $$

Back to Eq. (45), since $|\mathcal{E}_i| \le (\frac{2}{\varepsilon})^{n_i}$, involving the union bound gives us

$$ \begin{aligned} \mathbb{P}\{\|\mathbf{X}\| \ge t\} &\le \sum_{\mathbf{e}^{(1)} \in \mathcal{E}_1, \dots, \mathbf{e}^{(d)} \in \mathcal{E}_d} \mathbb{P}\{\mathbf{X}(\mathbf{e}^{(1)}, \dots, \mathbf{e}^{(d)}) \ge \frac{t}{2}\} \\ &\le \left(\frac{2d}{\log(3/2)}\right)^{\sum_{i=1}^{d} n_i} 2 \exp\left(-\frac{t^2}{8}\right), \end{aligned} $$

which yields the final bound for an appropriate choice of $t$. $\square$

LEMMA 5 (Woodbury matrix identity). Let $\mathbf{A} \in \mathbb{M}_n$, $\mathbf{B} \in \mathbb{M}_k$, $\mathbf{U} \in \mathbb{M}_{n,k}$ and $\mathbf{V} \in \mathbb{M}_{k,n}$, we have

$$ (\mathbf{A} + \mathbf{U}\mathbf{B}\mathbf{V})^{-1} = \mathbf{A}^{-1} - \mathbf{A}^{-1}\mathbf{U}(\mathbf{B}^{-1} + \mathbf{V}\mathbf{A}^{-1}\mathbf{U})^{-1}\mathbf{V}\mathbf{A}^{-1}. $$
REFERENCES

[1] ANANDKUMAR, A., GE, R., HSU, D., KAKADE, S. M. and TELGARSKY, M. (2014). Tensor decompositions for learning latent variable models. *Journal of machine learning research* **15** 2773–2832.

[2] BAIK, J., AROUS, G. B. and PÉCHÉ, S. (2005). Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices. *The Annals of Probability* **33** 1643–1697.

[3] BEN AROUS, G., HUANG, D. Z. and HUANG, J. (2021). Long Random Matrices and Tensor Unfolding. *arXiv preprint arXiv:2110.10210*.

[4] BEN AROUS, G., MEI, S., MONTANARI, A. and NICA, M. (2019). The landscape of the spiked tensor model. *Communications on Pure and Applied Mathematics* **72** 2282–2330.

[5] BENAYCH-GEORGES, F. and NADAKUDITI, R. R. (2011). The eigenvalues and eigenvectors of finite, low rank perturbations of large random matrices. *Advances in Mathematics* **227** 494–521.

[6] BIROLI, G., CAMMAROTA, C. and RICCI-TERSENGHI, F. (2020). How to iron out rough landscapes and get optimal performances: averaged gradient descent and its application to tensor PCA. *Journal of Physics A: Mathematical and Theoretical* **53** 174003.

[7] CAPITAINE, M., DONATI-MARTIN, C. and FÉRAL, D. (2009). The largest eigenvalues of finite rank deformation of large Wigner matrices: convergence and nonuniversality of the fluctuations. *The Annals of Probability* **37** 1–47.

[8] GOULART, J. H. D. M., COUILLET, R. and COMON, P. (2021). A Random Matrix Perspective on Random Tensors. *arXiv preprint arXiv:2108.00774*.

[9] HANDSCHY, M. C. (2019). Phase Transition in Random Tensors with Multiple Spikes, PhD thesis, University of Minnesota.

[10] HITCHCOCK, F. L. (1927). The expression of a tensor or a polyadic as a sum of products. *Journal of Mathematics and Physics* **6** 164–189.

[11] HUANG, J., HUANG, D. Z., YANG, Q. and CHENG, G. (2020). Power iteration for tensor pca. *arXiv preprint arXiv:2012.13669*.

[12] JAGANNATH, A., LOPATTO, P. and MIOLANE, L. (2020). Statistical thresholds for tensor PCA. *The Annals of Applied Probability* **30** 1910–1933.

[13] LANDSBERG, J. M. (2012). Tensors: geometry and applications. *Representation theory* **381** 3.

[14] LESIEUR, T., MIOLANE, L., LELARGE, M., KRZAKALA, F. and ZDEBOROVÁ, L. (2017). Statistical and computational phase transitions in spiked tensor estimation. In *Proc. IEEE International Symposium on Information Theory (ISIT)* 511–515.

[15] LIM, L.-H. (2005). Singular values and eigenvalues of tensors: a variational approach. In *Proc. IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing* 129–132.

[16] MARČENKO, V. A. and PASTUR, L. A. (1967). Distribution of eigenvalues for some sets of random matrices. *Mathematics of the USSR-Sbornik* **1** 457.

[17] MONTANARI, A. and RICHARD, E. (2014). A statistical model for tensor PCA. *arXiv preprint arXiv:1411.1076*.

[18] NOCEDAL, J. and WRIGHT, S. (2006). *Numerical optimization*. Springer Science & Business Media.

[19] PÉCHÉ, S. (2006). The largest eigenvalue of small rank perturbations of Hermitian random matrices. *Probability Theory and Related Fields* **134** 127–173.

[20] PERRY, A., WEIN, A. S. and BANDEIRA, A. S. (2020). Statistical limits of spiked tensor models. In *Annales de l'Institut Henri Poincaré, Probabilités et Statistiques* **56** 230–264. Institut Henri Poincaré.

[21] RABANSER, S., SHCHUR, O. and GÜNNEMANN, S. (2017). Introduction to tensor decompositions and their applications in machine learning. *arXiv preprint arXiv:1711.10781*.

[22] STEIN, C. M. (1981). Estimation of the mean of a multivariate normal distribution. *The Annals of Statistics* **11** 135–151.

[23] SUN, W. W., HAO, B. and LI, L. (2014). Tensors in Modern Statistical Learning. *Wiley StatsRef: Statistics Reference Online* 1–25.

[24] TAO, T. (2012). *Topics in random matrix theory* **132**. American Mathematical Soc.

[25] ZARE, A., OZDEMIR, A., IWEN, M. A. and AVIYENTE, S. (2018). Extension of PCA to higher order data structures: An introduction to tensors, tensor decompositions, and tensor PCA. *Proceedings of the IEEE* **106** 1341–1358.