A Random Matrix Analysis of Symmetric Tensor Unfolding

Introduction

We are interested in the following tensor model:

where $\lambda> 0$ is a signal-to-noise ratio (SNR), $v \in S^{n-1}$ is a planted signal and

is a random tensor of order $p$ and dimension $n$, obtained by symmetrizing a tensor $\mathbf{G}$ with independent standard Gaussian entries, ensuring rotational invariance. The goal is to quantify the performance as $n \to\infty$ of a certain estimator of $v$ called the unfolding algorithm—which amounts to extracting the leading singular vector of a matrix reshaping of the tensor—as a function of the SNR $\lambda$. As usual in this setting, we are also interested in determining the SNR regime under which one can partially recover the signal. An illustration of the problem is given in Figure 1.

Motivation

Model (1) is a spiked tensor model ([46]) and the problem of estimating its parameters $\lambda$ and $\mathbf{v}$ is commonly referred to as tensor principal component analysis (Tensor PCA). This generalization of PCA to tensors comes with many applications in machine learning (see reviews from Kolda and Bader [38] and Sidiropoulos et al. [54]). In particular, it is central in the estimation of latent variable models via the method of moments ([51]), with applications to independent component analysis, mixture models, topic models and hidden Markov models ([2]). It also arises in community detection, via hypergraph adjacency tensors ([67, 37, 43, 44]), in collaborative filtering, via tensor completion ([42, 47]) and many image processing tasks going beyond two-dimensional images ([61]). Finally, it serves as a typical model for complex high-dimensional nonconvex landscapes, such as those encountered in deep learning ([35, 66, 45, 10, 4] (Ge et al. 2018)).

Unfolding (sometimes called matricization, reshaping or flattening) is a classical operation in the analysis of tensor-based methods. It is a natural way of transforming a tensor into a matrix, enabling the use of well-established matrix analysis tools, such as eigenvalue decompositions, concentration inequalities [59] and more general random matrix results [58, 19]. Consequently, unfolding is widely used as a spectral relaxation of low-rank tensor decomposition [53, 46] (Mu et al. 2014). It is also used as a correlated initialization in many iterative algorithms—such as power iteration [46, 36], recursive unfolding [46] and higher-order orthogonal iteration (HOOI) [20]—that achieve exact recovery beyond a critical SNR threshold [41] (Feldman and Donoho 2023).¹ Hence, studying unfolding under spiked tensor models such as (1) is essential for understanding the statistical performance of low-rank tensor decomposition methods, coming with plenty of applications in machine learning.

State of the Art

Many works have analyzed spiked tensor models, identifying statistical and algorithmic thresholds on the SNR for signal recovery, under diverse estimation procedures. A statistical threshold refers here to the SNR above which a hypothesis test can distinguish a spiked model $\lambda v^{\otimes p}+Z$ from pure noise $Z$. An algorithmic threshold refers to the SNR above which a specified computationally-tractable class of estimators can achieve partial recovery. These notions need not coincide and are model- and algorithm-class dependent. The following recovery results are voluntarily simplified for expository purposes.²

Model (1) was originally introduced in Montanari and Richard [46], who proposed several estimators for the planted signal $v \in S^{n-1}$, as $n \to\infty$. They showed that maximum likelihood estimation (MLE) achieves partial recovery at $\lambda= \Omega(1)$—that is the statistical threshold [52, 34, 26]. However, solving the MLE problem itself is NP-hard [29], leaving the question of the existence of tractable estimators achieving partial recovery at $\lambda= \Omega(1)$ open.

Several polynomial-time alternatives for Tensor PCA have then been analyzed. Local methods such as power iteration, approximate message passing, Langevin dynamics and gradient descent have been shown to achieve partial recovery at $\lambda= \Omega(n^{(p-2)/2})$ under random initialization [46, 11]. Unfolding, partial tracing and local methods with “smart” initializations have been shown to achieve partial recovery at $\lambda= \Omega(n^{(p-2)/4})$ [46, 31, 32, 3, 12] (Feldman and Donoho 2023). Several analyses based on sum-of-squares relaxations and low-degree methods suggest that the latter threshold cannot be significantly improved with polynomial-time algorithms [31, 39]. This leads to a conjectured statistical–computational gap between intractable MLE and tractable algorithms like unfolding that remains open today; see Figure 2 for an illustration.

Looking more closely into the assumptions of Montanari and Richard [46] regarding tensor unfolding, the main results are actually derived under an asymmetric noise model, where the entries of the noise tensor are independent standard Gaussians. This importantly differs from the original noise model (2)—called Gaussian Orthogonal Tensor Ensemble (GOTE, Gurau [27])—whose symmetry yields dependent entries, preventing the application of classical results from random matrix theory. Montanari and Richard [46] conjecture that the same order of phase transition should hold for the symmetric model.

Some later works partially address this conjecture. Notably, Hopkins et al. [31] and Wein et al. [63] study symmetric models³ and obtain $\Omega(n^{(p-2)/4})$ recovery guarantees under variants of the tensor unfolding algorithm, while not providing either sharp transitions or precise overlaps. Follow-up works obtain sharp unfolding transitions and overlaps for even and odd orders [12, 41] (Feldman and Donoho 2023), but all under an asymmetric noise model that is not adapted to symmetric signals.

Hence, a random matrix analysis of the unfolding algorithm on the genuinely symmetric spiked tensor model of Montanari and Richard [46] is currently lacking in the Tensor PCA literature.

Contributions

We perform a random matrix analysis of the unfolding algorithm under the spiked tensor model (1), with even order ($p := 2q$) and square unfolding ($n^q \times n^q$). Specifically, we derive a deterministic equivalent (Couillet and Liao [19]) of the resolvent and deduce the weak convergence of the empirical spectral distribution towards a generalization of the Wigner semicircle distribution [64], that is a mixture of a scaled semicircle distribution and a Dirac distribution, whose weights depend on the tensor order $p$. Despite the non-classical structure of the resolvent and the limiting spectral distribution, we show that the unfolding algorithm undergoes a classical phase transition regarding the recovery of the lifted signal $\boldsymbol{v}^{\boxtimes q} \in\mathbb{S}^{n^q-1}$, with explicit critical value and overlap—$\boxtimes$ denotes the Kronecker product. See § 2 for notation and definitions.

More precisely, let $\mathbf{Z}$ be a Gaussian orthogonal tensor defined as in (2) and $\mathbf{Z}$ its $n^q \times n^q$ unfolding. Denoting by $f_{\mathrm{sc}}(x) := \sqrt{(4-x^2)_+}/(2\pi)$ the standard Wigner semicircle density and $\alpha_{p,q} := \sqrt{q!/(p-1)!}$, the empirical spectral measure of $\alpha_{p,q}^{-1}n^{-(p-2)/4}\mathbf{Z}$ converges weakly almost surely to $q!^{-1}f_{\mathrm{sc}}(x)\,dx+(1-q!^{-1})\delta_0(dx)$. Let $\mathbf{X}=\lambda\boldsymbol{v}^{\otimes p}+\mathbf{Z}$ be a spiked tensor as in (1), $\mathbf{X}$ its $n^q\times n^q$ unfolding and $(\hat{\lambda},\hat{\boldsymbol{v}})\in\mathbb{R}_+\times\mathbb{S}^{n^q-1}$ its leading eigenpair. Writing $\tau:=\alpha_{p,q}^{-1}n^{-(p-2)/4}\lambda$ and $\hat{\tau}:=\alpha_{p,q}^{-1}n^{-(p-2)/4}\hat{\lambda}$, one has, almost surely, as $n\to\infty$, the following sharp phase transitions: for $\tau\leq1$, $\hat{\tau}\to2$ and $|\hat{\boldsymbol{v}}^{\mathsf{T}}\boldsymbol{v}^{\boxtimes q}|^2\to0$ while for $\tau>1$, $\hat{\tau}\to\tau+\tau^{-1}$ and $|\hat{\boldsymbol{v}}^{\mathsf{T}}\boldsymbol{v}^{\boxtimes q}|^2\to1-\tau^{-2}$.

Our results confirm the open conjectures from Montanari and Richard [46] and Hopkins et al. [31] on symmetric spiked tensor unfolding: there is indeed a phase transition for partial recovery at $\lambda\asymp n^{(p-2)/4}$ and we now know its precise constant: $\alpha_{p,q}:=\sqrt{(p/2)!/(p-1)!}$. Hence, the high-order symmetries do not affect significantly the order of magnitude of the phase transition for partial recovery with respect to asymmetric models. They however significantly affect the spectral distribution, which becomes a mixture of a semicircle and a Dirac distribution. Moreover, although the phase transition is at the conjectured order of magnitude, $\lambda\asymp n^{(p-2)/4}$, it appears at a different scale from the asymmetric model and the effect increases with the tensor order according to $\alpha_{p,q}$. ⁴ Therefore, a simplified analysis focused on a “symmetric signal + asymmetric noise” model (Montanari and Richard, 2014; Kim et al., 2017; Feldman and Donoho, 2023) does not lead to the same phase transition thresholds and overlaps obtained under the genuine spiked tensor model (1) and misses important spectral phenomena caused by the loss of independence induced by the symmetries.

The proof relies on random matrix theory (RMT), exposing important phenomena that cannot be revealed via classical machine learning tools—notably, a sharp phase transition in spectral behavior (the so-called BBP transition of Baik et al. 2005), with precise eigenvector overlaps. The high-order symmetries in the GOTE break the fundamental independence assumption required in standard random matrix theory. To alleviate these dependencies, we project the unfolded GOTE onto symmetric tensor subspaces, which boils down to removing duplicated rows and columns and applying a specific scaling on the entries. The whole procedure, that we name symmetric unfolding, preserves the eigenstructure (up to a simple isometry), while removing the most impacting dependencies. The resulting matrix spectrally behaves like a standard Wigner matrix, from which we can precisely characterize the limiting spectral distribution and overlaps.

The paper is organized as follows. In § 2, we introduce our notation, the spiked tensor model, the unfolding map and the main random matrix tools needed to establish our results. In § 3, we state the main results of the paper: the analysis of unfolded symmetric tensors under the framework of hypersymmetric matrices, the deterministic equivalent, the limiting spectral distribution and the eigenvector overlap. We run simulations and show a strong match with theory, even for relatively small dimensions, $n \in\{20, 50, 100\}$. In § 4, we discuss and put the previous results into perspective. In Appendix § 6, we give the proofs of the main results. In Appendix § 7, we extend our results and simulations to the i.i.d. symmetric standard Gaussian tensor model. In Appendix § 8, we extend our results to multiple spikes, under the canonical polyadic (CP) model. In Appendix § 9, we analyze a recursive unfolding procedure to estimate the planted signal $\boldsymbol{v} \in\mathbb{S}^{n-1}$ from the lifted signal $\boldsymbol{v}^{\boxtimes q} \in\mathbb{S}^{n^q-1}$; we conjecture a sharp phase transition at $\lambda= \alpha_{p,q} n^{(p-2)/4}$, from null to exact recovery. In Appendix § 10, we study odd order $p$ and non-square unfolding ($n^k \times n^{p-k}$, with $1 \leq k < p/2$); we conjecture a sharp phase transition at $\lambda= \gamma_{p,k} n^{(p-2)/4}$, with $\gamma_{p,k} := \left(k!(p-k)!/(p-1)!^2\right)^{1/4}$, resulting in interesting statistical–computational tradeoffs.

Preliminaries

Notation

Throughout this paper, we adopt the following notation. The dimension of the tensor is denoted by $n$ and assumed to be large, while the order is denoted by $p$ and assumed to be fixed and even: $p := 2q$ for some (fixed) positive integer $q$. For convenience, integer intervals are shortened as $[n] := \{1, \ldots, n\}$ and multi-indices as $\boldsymbol{i} := (i_1, \ldots, i_p) \in[n]^p$. The set of permutations of a multi-index $\{i_1, \ldots, i_p\}$ is written $\mathfrak{S}_{i_1,\ldots,i_p}$; for short, we write $\mathfrak{S}_p := \mathfrak{S}_{1,\ldots,p}$. The cardinal of a set is denoted by $\#$. We also denote by $\boldsymbol{i} \sqcup\boldsymbol{j}$ the multiset union of two multi-indices $\boldsymbol{i}$ and $\boldsymbol{j}$. Scalars, vectors, matrices and tensors are respectively written with the following fonts: $a$, $\boldsymbol{a}$, $A$ and $\mathcal{A}$. We follow the common convention of omitting the indices in sequences of (random or deterministic) matrices: $a \in\mathbb{R}^n$ and $A \in\mathbb{R}^{n \times n}$ should be respectively understood as $a_n \in\mathbb{R}^n$ and $A_n \in\mathbb{R}^{n \times n}$.

Outer and Kronecker products of vectors $\boldsymbol{a}_1, \ldots, \boldsymbol{a}_p \in\mathbb{R}^n$ are respectively denoted by $\boldsymbol{a}_1 \otimes\cdots\otimes\boldsymbol{a}_p \in(\mathbb{R}^n)^{\otimes p}$ and $\boldsymbol{a}_1 \boxtimes\cdots\boxtimes\boldsymbol{a}_p \in\mathbb{R}^{n^p}$; in addition, we write $\boldsymbol{a}^{\otimes p} = \boldsymbol{a} \otimes\cdots\otimes\boldsymbol{a}$ and $\boldsymbol{a}^{\boxtimes p} = \boldsymbol{a} \boxtimes\cdots\boxtimes\boldsymbol{a}$, where the product is repeated $p$ times. The spectrum of a symmetric matrix $M \in\mathbb{R}^{n \times n}$ is denoted by $\operatorname{Sp}(M)$; its eigenvalues, $\lambda_1(M) \geq\cdots\geq\lambda_n(M)$, are repeated according to their multiplicity and sorted in decreasing order.

As we work in the asymptotic regime where $n \to\infty$, we use the classical Bachmann–Landau “big $O$” notation. We write $f(n) = O(g(n))$ if there exists a positive constant $C$ such that $|f(n)| \leq Cg(n)$ as $n \to\infty$. Similarly, we write $f(n) = \Omega(g(n))$ if $g(n) = O(f(n))$, $f(n) \asymp g(n)$ if $f(n) = O(g(n))$ and $g(n) = O(f(n))$. Finally, we write $f(n) = o(g(n))$ if $f(n)/g(n) \to0$, $f(n) = \omega(g(n))$ if $f(n)/g(n) \to\infty$ and $f(n) \sim g(n)$ if $f(n)/g(n) \to1$.

Symmetric Tensors

The spiked tensor model (1) we consider in this work has the following entries:

with $\lambda> 0$, $v \in S^{n-1}$ and $G_{i_1\cdots i_p} \overset{\mathrm{iid}}{\sim} \mathcal{N}(0,1)$. The random tensor $\mathbf{Z} := \frac{1}{p!}\sqrt{\frac{p}{n}}\sum_{\sigma\in\mathfrak{S}_p}\mathbf{G}^{\sigma}$ belongs to the Gaussian Orthogonal Tensor Ensemble; we write $\mathbf{Z} \sim\mathrm{GOTE}_{p,n}$.⁵

One can show that the entries of $\mathbf{Z}$ are independent normal—up to tensor symmetries—and their variance depends on the multiplicities of the multi-index [15]. More precisely, for $i_1 \le\cdots\le i_p$ and $\sigma\in\mathfrak{S}_p$, one has

Hence, the GOTE is a generalization of the (matrix) Gaussian Orthogonal Ensemble (GOE), which exhibits variance $1/n$ off-diagonal and $2/n$ on-diagonal, to $p \ge2$.

From a combinatorial viewpoint, writing $\mathbf{m}$ the multiset of multiplicities of the indices in $i \in[n]^p$, one has

The latter expression is called a multinomial coefficient—a multi-way generalization of a binomial coefficient. The set of $p$-indices with multiplicities $\mathbf{m}$ therefore has cardinal $\binom{n}{\#\mathbf{m}}\binom{p}{\mathbf{m}}$. This implies that multi-indices with distinct indices predominate in high dimensions; indeed, whenever $\mathbf{m} \ne\{1,\ldots,1\}$, one has $\#\mathbf{m} < p$ therefore $\binom{n}{\#\mathbf{m}}/\binom{n}{p} \to0$.

Let $\operatorname{Sym}^{p}(\mathbb{R}^{n})$ denote the space of symmetric tensors of order $p$ and dimension $n$—we remove the superscript when $p=2$ for symmetric matrices. Then $\operatorname{Sym}^{p}(\mathbb{R}^{n})$ is a linear subspace of $(\mathbb{R}^{n})^{\otimes p}$ of dimension $n_p := \dim(\operatorname{Sym}^{p}(\mathbb{R}^{n})) = \binom{n+p-1}{p}$; this can be established via standard “stars-and-bars” arguments. A canonical orthonormal basis of $\operatorname{Sym}^{p}(\mathbb{R}^{n})$ is

where $\mathbf{E}_j \in(\mathbb{R}^{n})^{\otimes p}$ is the unit tensor with zeros everywhere except in position $(j_1,\ldots,j_p)$.

Tensor Unfolding

Let us define the following bijective function between multi-indices and indices:

The map $\varphi_k$ induces a natural bijection between order-$p$ tensors and matrices that is commonly called unfolding:

We identify $\operatorname{Mat}_p$ with the classical vectorization $\operatorname{vec}: (\mathbb{R}^{n})^{\otimes p} \to\mathbb{R}^{n^p}$. In this paper, we choose $k = \lfloor p/2 \rfloor= q$. This corresponds to a square unfolding, yielding an $n^q \times n^q$ symmetric matrix—like a classical Wigner matrix but with extra symmetries.

Remark 1 (Choice of Unfolding) Other choices of $k$ for unfolding would have led to so-called long [12] or disproportionate [22] random matrices. In that case, the dimensions $(n^k \times n^{p-k})$ are no longer comparable; one cannot apply classical RMT results [14] and needs to use more sophisticated tools [1]. In addition, the choice of an unfolding that is “as square as possible” is recommended in Montanari and Richard [46] and Mu et al. (2014), as it minimizes the operator norm of the unfolded noise. It is worth noting that this heuristic is later nuanced by Ben Arous et al. [12], which states that the SNR threshold does not depend on the choice of $k$ in the asymmetric model and therefore recommend $k = 1$ to avoid an extra recursive unfolding step to recover the planted signal. We shed light on these opposite points of view in § 4 and Appendix § 10.

With a slight abuse of notation, we will often replace a multi-index $\boldsymbol{i} := (i_1, \ldots, i_k) \in[n]^k$ by its associated index $i := \varphi_k(\boldsymbol{i}) \in[n^k]$ and vice-versa. Likewise, we will often assimilate a tensor $\mathbf{X} \in(\mathbb{R}^{n})^{\otimes p}$ to its associated matrix $\mathbf{X} := \operatorname{Mat}_k(\mathbf{X}) \in\mathbb{R}^{n^k \times n^{p-k}}$.

Useful Tools from Random Matrix Theory

Random matrix theory is a rich mathematical field at the intersection of probability theory, linear algebra and complex analysis. It provides tools to study the eigenstructure of large matrices with random entries—notably their limiting spectral distribution, isolated eigenvalues and associated eigenspaces—enabling the theoretical performance analysis of many machine learning algorithms. The RMT framework offers several advantages over more classical theoretical tools in machine learning, often referred to as the blessing of dimensionality [21]. In high dimensions, a priori random spectral quantities often become asymptotically deterministic through concentration of measure, while universality makes their limiting behavior largely insensitive to the precise noise distribution. Moreover, RMT accommodates structured and dependent noise, as well as certain nonlinear spectral statistics that are central in machine learning (kernel, activation function). These features make it particularly well suited to characterizing phase transitions and partial recovery in high-dimensional planted-signal models. This subsection recalls the main results from RMT that will be used throughout this paper. Its content is mostly based on Couillet and Liao [19], where the reader can find more details.

The first object we review is the resolvent, which is central in RMT, as it provides information on the eigenstructure.

Hence, the resolvent gives access to the empirical spectral measure and the eigenvector overlap of deterministic matrices.⁶ When these are random, the eigenstructure can be characterized via deterministic equivalents [25, 28], i.e. deterministic matrices that “behave” similarly asymptotically with respect to certain linear and bilinear forms, as defined next.

Such deterministic equivalents can be obtained by showing the concentration of the expressions in (8) around their expectation. To that end, we state two fundamental probabilistic identities, where $\mathbb{E}[\cdot]$ and $\mathbb{V}[\cdot]$ denote respectively the expectation and the variance.

Now that the main tools are given, let us outline a typical approach to characterize eigenstructure via RMT. The first step is to establish a deterministic equivalent for the resolvent matrix $\mathbf{Q}_{\mathrm{M}} \in\mathbb{R}^{n \times n}$, a natural candidate being its expectation, $\mathbb{E}[\mathbf{Q}_{\mathrm{M}}]$. To achieve this, we derive a self-consistent equation for the resolvent and simplify it using Stein’s Lemma (Lemma 4). This yields an explicit asymptotic expression for the expected resolvent. Next, we demonstrate the concentration of linear and bilinear forms of the resolvent around their expectation via the Nash–Poincaré inequality (Lemma 5). If the resulting variance bound drops fast enough, then this allows establishing almost sure deterministic equivalence (Definition 3). The latter yields the limiting spectral distribution via inverse Stieltjes transform and the eigenvector overlap $|\widehat{\mathbf{u}}^{\mathsf{T}}\mathbf{u}|$ via Cauchy’s integral formula (see Lemma 2). Finally, the integral expression can be further simplified using the Sherman–Morrison formula, which relates the resolvent of the spiked model to that of the noise model. We specialize this scheme to the symmetric unfolding in Appendix § 6.

Main Results

This section states the main results, namely the impact of tensor symmetries on the spectrum, the deterministic equivalent of the resolvent—consequently, the limiting spectral distribution—and the spiked eigenvalues and overlaps. The proofs are given in Appendix § 6.

Hypersymmetric Matrices and Symmetric Unfolding

Let $\mathbf{Z} \sim\operatorname{GOTE}_{p,n}$ and $\mathbf{Z} := \operatorname{Mat}_{q}(\mathbf{Z}) \in\mathbb{R}^{n^{q} \times n^{q}}$. As $\mathbf{Z}$ is a symmetric matrix with (mostly) independent Gaussian entries, one could think that its eigenvalues are almost-surely simple and nonzero—by applying Sard’s theorem to the discriminant polynomial [18]. However, the symmetries in $\mathbf{Z}$ lead its unfolding $\mathbf{Z}$ to have hypersymmetries—i.e., symmetries beyond the classical matrix ones—including notably duplicated columns. The next definition formalizes this notion of hypersymmetry in symmetric tensors.

An important consequence of hypersymmetry is duplicated rows and columns: for $(i',j') \in\mathfrak{S}_{i} \times\mathfrak{S}_{j}$, one has $M_{i\cdot} = M_{i'\cdot}$ and $M_{\cdot j} = M_{\cdot j'}$. Intuitively, removing these duplicates does not change the eigenstructure. We formalize this observation by projecting the row and column spaces onto symmetric tensor subspaces.

Let $\mathcal{I}_{q} := \{(i_{1},\ldots,i_{q}) \in[n]^{q}: i_{1} \leq\cdots\leq i_{q}\}$ and $n_{q} := |\mathcal{I}_{q}| = \binom{q+n-1}{q}$. For $i \in\mathcal{I}_{q}$, define $\mathbf{u}_{i} := \operatorname{vec}(\mathbf{U}_{i}) \in\mathbb{S}^{n^{q}-1}$ and $\mathbf{U} := (\mathbf{u}_{i})_{i \in\mathcal{I}_{q}} \in\mathbb{R}^{n^{q} \times n_{q}}$. The vectors $(\mathbf{u}_{i})_{i \in\mathcal{I}_{q}}$ form an orthonormal basis of $\operatorname{Sym}^{q}(\mathbb{R}^{n})$; hence $\mathbf{U}^{\mathsf{T}}\mathbf{U} = \mathbf{I}_{n_{q}}$ and $\mathbf{\Pi} := \mathbf{U}\mathbf{U}^{\mathsf{T}} = \left(\frac{1}{\sqrt{\#\mathfrak{S}_{i}\#\mathfrak{S}_{j}}}\mathbf{1}_{\{j\in\mathfrak{S}_{i}\}}\right)_{i,j\in[n^{q}]}$ is the orthogonal projector onto $\operatorname{Sym}^{q}(\mathbb{R}^{n})$.

We call $\Phi_{q} : \mathbf{T} \mapsto\mathbf{U}^{\mathsf{T}}\operatorname{Mat}_{q}(\mathbf{T})\mathbf{U}$ the symmetric unfolding. The previous result implies that the nonzero eigenpairs of $\operatorname{Mat}_{q}(\mathbf{T})$ and $\Phi_{q}(\mathbf{T})$ are isometric; we may therefore abusively switch between the two representations in the following.

Although elementary, Proposition 7 is actually extremely informative: the eigenstructure of $\operatorname{Mat}_{q}(\mathbf{T})$—notably the limiting spectral distribution, the spiked eigenvalues and overlaps—can be simply obtained via $\Phi_{q}(\mathbf{T})$, a compact version obtained by removing duplicated columns and scaling by the number of permutations of the indices.⁷ From a practical standpoint, this yields a computational speedup by a factor $(n^q/n_q)^3 \sim q!^3$ for eigendecomposition. From a theoretical standpoint, the projection onto the symmetric subspace reduces the number of dependencies between the entries. Let

represent the equivalence class of entries that are determined by the entry $(i,j) \in\mathcal{I}_q^2$. Then $\#\mathcal{E}(i,j) \le\binom{p}{q} = p!/q!^2$ and we later show that these extra dependencies do not impact the asymptotic eigenstructure; see § 4 for a broader discussion.

Deterministic Equivalent

The tools presented in previous sections (2.4 and 3.1) to analyze random matrices and hypersymmetric ones can now be combined to characterize the spectrum of square unfoldings of symmetric tensors under the spiked model (1). The following theorem gives a deterministic equivalent of the noise resolvent.

The theorem says that, on the symmetric tensor subspace, the normalized GOTE unfolding has the same deterministic resolvent as a classical GOE matrix. As a consequence of Proposition 7, we get a deterministic equivalent of the GOTE unfolding.

The proof follows the scheme exposed in § 2.4. The symmetric unfolding removes the deterministic null space created by repeated rows and columns. After the normalization by $\alpha_{p,q}n^{(p-2)/4}$, the entries have the same leading variance as those of a standard Wigner matrix of dimension $n^q$; the remaining correlations are indexed by the finite equivalence classes (13). The proof shows, using Stein’s lemma, that all dependencies other than the usual matrix transpose symmetry are asymptotically negligible in the resolvent equation.

Remark 11 (Normalization Factor). One could expect a normalization by $\sqrt{n^q}$ given the resolvent matrix size ($n^q \times n^q$). The $n^{(p-2)/4}=\sqrt{n^q}/\sqrt{n}$ factor we get instead comes from the $1/\sqrt{n}$ factor in the definition of the noise model (2). Although slightly more cumbersome, we decide to keep this normalization in order to stick with the original work of Montanari and Richard [46]. Several follow-up works indeed stick with this normalization, making the conjectured expression “$n^{(p-2)/4}$” central in the Tensor PCA community [10, 52, 26, 12, 41], while some other remove the noise tensor normalization by $\sqrt{n}$ for conciseness (Hopkins et al., 2015; Kim et al., 2017; Feldman and Donoho, 2023), getting “$n^{p/4}$”-like thresholds.

Remark 12 (Spectral Norm). The weak convergence of the empirical spectral distribution towards the limiting spectral distribution is not sufficient to guarantee that no empirical eigenvalue “escapes” from the semicircle bulk support, which is crucial for estimation (see Couillet and Liao 2022, § 2.3). In Appendix § 6.§ 3, we show that the spectral edge of the (normalized) symmetrically unfolded GOTE is the standard Wigner spectral edge [5], by a high-moment argument adapted to the high-order dependencies (13).

Phase Transition

The previous results characterize the deterministic equivalent and the limiting spectral distribution of the unfolded noise tensor $\mathbf{Z} \sim\mathrm{GOTE}_{p,n}$. Coming back to the spiked tensor model $\mathbf{X}=\lambda v^{\otimes p}+\mathbf{Z}$ (1), we now study the recovery of the lifted signal $v^{\otimes q}\in S^{n^q-1}$ by eigendecomposition of the unfolded tensor $\mathbf{X}:=\operatorname{Mat}_{q}(\mathbf{X})$. The next theorem establishes a sharp phase transition with precise overlaps, as a function of the SNR $\lambda$.

Writing $\tau:= \lambda/(\alpha_{p,q} n^{(p-2)/4})$ and $\hat{\tau} := \hat{\lambda}/(\alpha_{p,q} n^{(p-2)/4})$, with $\alpha_{p,q} := \sqrt{q!/(p-1)!}$, one has

with $v^{\odot q}:=U^{\top}v^{\otimes q}\in\mathbb{S}^{n_q-1}$. These results apply verbatim to $\mathbf{X}:=\operatorname{Mat}_q(\mathbf{X})\in\operatorname{Sym}(n^q)$.

The previous results generalize the GOE to higher orders. Indeed, when $p=2$, one has $q!=1$ and $\alpha_{p,q}=1$, therefore the empirical spectral distribution is precisely the Wigner semicircle distribution and the associated overlaps follow [13]. We note that the support of the semicircle shrinks—as $\sqrt{q!/(p-1)!}$—and the proportion of zero eigenvalues grows—as $1-q!^{-1}$—with the tensor order $p$.

Simulations

To confront the asymptotic theoretical results with finite-sample observations, we consider a spiked tensor model $\mathbf{X}=\lambda v^{\otimes p}+\mathbf{Z}$ (1) with $p=4$. We report the theoretical and simulated spectral distribution of the symmetrically unfolded GOTE, $\mathbf{M}:=\alpha_{p,q}^{-1}n^{-(p-2)/4}\Phi_q(\mathbf{Z})$, in Figure 3, with $n=100$. Moreover, we report the eigenvalues and overlaps in Figure 4, with $n\in\{20,50,100\}$. For that experiment, we fix the signal $v\in\mathbb{S}^{n-1}$ via uniform sampling on the sphere, then, ten times independently, we sample $\mathbf{Z}$ according to (2) and increase $\tau:=\lambda/(\alpha_{p,q}n^{(p-2)/4})$ linearly from 0 to 2.5. We see that the empirical eigenvalues are well estimated by the limiting spectral distribution. Both leading eigenvalue and overlap are increasingly well estimated by the theoretical ones beyond phase transition, as $\lambda$ grows. The phase transition is notably already well captured for $n=50$.

Discussion

We have studied the spectral behavior of the unfolding algorithm under the spiked symmetric tensor model (1). The strategy was to project the usual unfolding onto symmetric tensor subspaces, after which the matrix spectrally behaves like a Wigner matrix; this was rigorously proven using Stein’s lemma and Cauchy–Schwarz’s inequality on tensor slices with shared indices. We showed that, after normalization, the symmetric unfolding has a semicircular limiting spectrum and that its low-rank perturbation undergoes the usual BBP transition, with a sharp threshold $\lambda= \alpha_{p,q} n^{(p-2)/4}$ and an explicit overlap, settling the conjectured scaling of Montanari and Richard (2014). Hence, despite its simplicity, the unfolding algorithm ranks among the best known polynomial-time estimators for symmetric Tensor PCA problems. To the best of our knowledge, this is the first random matrix analysis of the unfolding algorithm under the genuinely symmetric model of Montanari and Richard (2014) at that level of precision.

Several questions remain open. First, our analysis suggests a broader universality principle for tensor unfoldings. For a general order $p$ and matricization $\operatorname{Mat}_k(\mathbf{T})$, one can similarly project the lines and columns onto $\operatorname{Sym}^k(\mathbb{R}^n)$ and $\operatorname{Sym}^{p-k}(\mathbb{R}^n)$, obtaining a rectangular matrix $\mathbf{\Phi}_k(\mathbf{T}) \in\mathbb{R}^{n^k \times n^{p-k}}$. For $k < p-k$, a direct computation of the entrywise variance shows that, after normalization, the resulting model can be written in the form $n_{p-k}^{-1/2}\mathbf{X} = \theta\mathbf{v}^{\odot k}\mathbf{v}^{\odot p-k\mathsf{T}} + n_{p-k}^{-1/2}\mathbf{Z}$, with $\theta:= \lambda\sqrt{(p-1)!/(k!n^{p-k-1})}$, where $\mathbf{Z}$ has “mostly” independent unit-variance entries. This is precisely the scaling of the disproportionate spiked matrix model of Feldman (2023), with matrix dimensions $n^k \times n^{p-k}$. If the remaining dependencies induced by tensor symmetry (13) are spectrally negligible, like in Lemma 17, their Theorem 1 would therefore suggest a sharp transition at

and analogous formulas for the outlier singular value and the overlap with the left singular vector, corresponding to the shorter side. Interestingly, the polynomial order $n^{(p-2)/4}$ is independent of the unfolding, as in the asymmetric setting [12] (Feldman and Donoho 2023), but the sharp constant is not: it is minimized by an “as square as possible” unfolding.⁸ Moreover, since the random matrix results hold under substantially more general assumptions than Gaussianity—see Benaych-Georges and Nadakuditi [13] for square unfolding and Feldman [22] for non-square unfolding—it is natural to expect that the same phenomenon extends to non-Gaussian symmetric noise under suitable moment and tail conditions. We conjecture that the symmetrically projected unfoldings indeed satisfy the corresponding disproportionate random matrix asymptotics; numerical simulations and additional comments are provided in Appendix § 10.

It would also be interesting to extend the present framework beyond the additive spiked model. Natural examples include moment and cumulant tensors [40, 57] arising in latent variable models [2]. In each case, the main question is whether an appropriate symmetric unfolding reduces the problem to a tractable random matrix model with a sharp spectral transition.

Symmetric unfolding is related to other spectral approaches to Tensor PCA, notably the Kikuchi hierarchy [63]. While symmetric unfolding exploits a single structured matricization, Kikuchi matrices introduce richer families of tensor transformations and provide a hierarchy of statistical and computational tradeoffs. Understanding the relation between these constructions, possibly through the recent random matrix and free-probabilistic techniques of Bandeira et al. [8, 7], could help understanding the impact of high-order symmetries on spectral transitions. Related questions arise for trace-invariant methods [50], which exploit tensor-network contractions.

Finally, our results raise the more fundamental question of whether symmetry can significantly modify the spectral behavior. In the present setting, the dependencies created by high-order symmetries (13) disappear at leading order after projection, resulting in a similar spectral behavior as the corresponding asymmetric model. This may not hold in general: localized signal priors or more-structured noise distributions may lead high-order dependencies to survive asymptotically, significantly impacting the eigenstructure.

Acknowledgments and Disclosure of Funding

The authors would like to thank Afonso S. Bandeira for insightful discussions around symmetric unfolding, Kikuchi matrices and intrinsic freeness. This work was supported by the French National Research Agency (ANR) via project LATENT (grant ANR-25-CE48-7508-01) and by the AI Interdisciplinary Institute ANITI funded by the France 2030 program (grant ANR-23-IACL-0002).

Proofs

Symmetric Unfolding (Proposition 7)

Proof. By definition of $\mathbf{U}$, one has $\Pi_{ij} := [\mathbf{U}\mathbf{U}^{\mathsf{T}}]_{ij} = \mathbf{1}_{\{j\in\mathcal{S}_i\}}/\#\mathcal{S}_i$. Moreover, as $\mathbf{M}$ is hypersymmetric, every row and column is constant on permutation orbits, so $\mathbf{M} = \mathbf{\Pi}\mathbf{M}\mathbf{\Pi} = \mathbf{U}\mathbf{U}^{\mathsf{T}}\mathbf{M}\mathbf{U}\mathbf{U}^{\mathsf{T}}$, with

Finally, since $\mathbf{U}^{\mathsf{T}}\mathbf{U} = \mathbf{I}_{n_q}$, the nonzero eigenpairs are transported by the isometric embedding $\mathbf{U}$. The remaining $n^q-n_q$ dimensions belong to $\operatorname{ker}(\mathbf{M})$.

Deterministic Equivalent (Theorem 9)

Proof. The starting point is the resolvent identity $\mathbf{M}\mathbf{Q} = \mathbf{I}_{n_q} + z\mathbf{Q}$. Applying Stein’s equality (Lemma 4), one gets

with

Noting that $\#\mathcal{S}_{i\sqcup k}=p!\#\mathcal{S}_i\#\mathcal{S}_k/q!$, for any “generic” $k\in\mathcal{I}_q$ such that $i\sqcap k=\varnothing$, one gets $p\#\mathcal{S}_i\#\mathcal{S}_k/(\alpha_{p,q}^{2}n^{p/2}\#\mathcal{S}_{i\sqcup k})=q!/n^{p/2}\sim n_q^{-1}$. Consequently,

Equation (20) is rather simple, but the sum over permutations of multi-indices it involves prevents us from obtaining a self-consistent equation on $\mathbf{Q}$. Fortunately, we can show that the only permutation yielding non-negligible terms is the trivial one, matching $(i,k)$ with $(i',k')$; indeed, otherwise $Q_{kk'}$ and $Q_{i'j}$ have shared indices, leading the associated summation over $k$ to be negligible with respect to the right-hand side.

Proof. One has

with the tensor slice notation

Let $\mathbf{A}'\in(\mathbb{R}^{n})^{\otimes s}$ be such that $A'_{a_1\cdots a_s}=A_{a_1\cdots a_s}$ if $(a_1,\ldots,a_s)\in\mathcal{I}_s$ and zero otherwise. Writing $\mathbf{A}:=\operatorname{Mat}_{n^{\lfloor s/2\rfloor}}(\mathbf{A}')$, one has $\|\mathbf{A}\|_{\mathrm{F}}=\|\mathbf{A}\|_{\mathrm{F}}\leq\sqrt{n^{\lfloor s/2\rfloor}}\|\mathbf{A}\|_{\mathrm{op}}=O(\sqrt{n^{\lfloor s/2\rfloor}})$. Similarly, one gets $\|\mathbf{B}\|_{\mathrm{F}}=O(\sqrt{n^{\lfloor s/2\rfloor}})$. We conclude via the Cauchy–Schwarz inequality:

Hence, as soon as $k$ and $i'$ share one index, the associated sum over $k$ is $O(n^{q-s}n^{\lfloor s/2\rfloor})=O(n^{q\eta-\lceil s/2\rceil})=O(n^{q-1})$. Consequently, one has

Lemma 17 implies that high-order dependencies are negligible, from the expected resolvent’s point of view. Applying it to the equation (20) yields

Writing $s(z):=n_q^{-1}\operatorname{tr}(\mathbf{Q})$ and using $\mathbf{M}\mathbf{Q}=\mathbf{I}_{n_q}+z\mathbf{Q}$, we get

Taking the normalized trace yields, as $n\to\infty$, the self-consistent equation $s^2+zs+1=0$, from which we recognize the classical Stieltjes transform of the semicircular distribution.

Substituting this relation back into (23) yields, for every deterministic $\mathbf{A},a,b$ such that $\|\mathbf{A}\|_{\mathrm{op}}=\|a\|_2=\|b\|_2=1$, $\frac{1}{n_q}\operatorname{tr}\left(\mathbf{A}\left(\mathbb{E}[\mathbf{Q}]-s_{\mathrm{sc}}(z)\mathbf{I}\right)\right)\to0$ and $a^\top\left(\mathbb{E}[\mathbf{Q}]-s_{\mathrm{sc}}(z)\mathbf{I}\right)b\to0$.

Hence, $\overline{\mathbf{Q}}:=s_{\mathrm{sc}}\mathbf{I}_{n_q}$ is a good candidate for deterministic equivalence of $\mathbf{Q}$ (Definition 3). However, the previous results hold in expectation; to show concentration, we apply the Nash–Poincaré inequality (Lemma 5) to the functions $f(x)=\operatorname{tr}(\mathbf{A}\mathbf{Q})$ and $g(x)=a^\top\mathbf{Q}b$, where $\mathbf{x}:=\operatorname{vec}(M_i)_{i\in\mathcal{I}_p}$. Using again the derivative (19), we get

Noticing that $\mathbf{a}^{\mathsf{T}}\mathbf{Q}\mathbf{b}=\operatorname{tr}(\mathbf{b}\mathbf{a}^{\mathsf{T}}\mathbf{Q})$, we can apply the same sequence of equalities for the concentration of bilinear forms of the resolvent. Writing $A_{kl}=b_k a_l$, we get

Hence, one has $\nabla\left[\frac{1}{n^q}\operatorname{tr}(\mathbf{A}(\mathbf{Q}-\overline{\mathbf{Q}}))\right]=O(n^{-p})$ and $\nabla[\mathbf{a}^{\mathsf{T}}(\mathbf{Q}-\overline{\mathbf{Q}})\mathbf{b}]=O(n^{-q})$. The variance bounds are both summable for $q\ge2$, therefore the convergence is almost sure by virtue of Chebyshev’s inequality and Borel–Cantelli’s lemma. For $q=1$, the matrix is an ordinary GOE and the almost-sure resolvent convergence is classical (see, for instance, Theorem 2.5 of Couillet and Liao 2022).

This proves the deterministic equivalent (14). The almost-sure convergence of the Stieltjes transform to $s_{\mathrm{sc}}$ follows. The Stieltjes continuity theorem then gives weak convergence to the standard semicircle law (15).

Proof [Proof of Corollary 10] The full unfolding satisfies $\alpha_{p,q}n^{(p-2)/4}\mathbf{M}=\mathbf{U}\Phi_q(\mathbf{Z})\mathbf{U}^{\top}$ on $\operatorname{Im}(\mathbf{U})$ and is zero on its orthogonal complement. Since $n_q/n^q\to1/q!$, its empirical law is the mixture of the semicircle scaled by $\alpha_{p,q}$ and a Dirac mass of weight $1-1/q!$.

Spectral Norm (Remark 12)

We establish the spectral norm of the symmetrically unfolded GOTE by adapting the standard high-moment argument [see 58] and showing that the additional identifications induced by tensor symmetry are negligible.

Proof Expanding the trace and applying Wick’s formula gives

where $\mathcal{P}_2(2k)$ denotes the set of pair partitions of $[2k]$ and $\mathbf{i}_{2k+1}=\mathbf{i}_1$. A covariance is nonzero only if $\mathbf{i}_a\sqcup\mathbf{i}_{a+1}=\mathbf{i}_b\sqcup\mathbf{i}_{b+1}$. The standard identifications are $(\mathbf{i}_a,\mathbf{i}_{a+1})=(\mathbf{i}_b,\mathbf{i}_{b+1})$ and $(\mathbf{i}_a,\mathbf{i}_{a+1})=(\mathbf{i}_{b+1},\mathbf{i}_b)$; all other identifications are called non-standard.

We first show that, for every $\varepsilon>0$, uniformly for $k\le c_\varepsilon\ln n$,

To this end, write $W(\pi,\mathbf{i}_1,\ldots,\mathbf{i}_{2k}):=\prod_{\{a,b\}\in\pi}\mathbb{E}[\mathbf{M}_{\mathbf{i}_a\mathbf{i}_{a+1}}\mathbf{M}_{\mathbf{i}_b\mathbf{i}_{b+1}}]$. The covariance structure of the symmetric unfolding gives, for some constants $A_q,B_q$,

Indeed, each non-standard identification leads to a factor $O(n^{-1})$ with respect to the associated standard identification—see Lemma 17—and the $k^{B_q}$ factor accounts for the possible choices of locations of the non-standard covariances. The standard Wigner moment count gives, uniformly for $k\le c_\varepsilon\ln n$,

Together with (29), this proves (28) for all sufficiently large $n$. Now let $k_n=\lfloor c_\varepsilon\ln n\rfloor$, with $c_\varepsilon$ sufficiently large. Since $\|\mathbf{M}\|_{\mathrm{op}}^{2k_n}\le\operatorname{tr}(\mathbf{M}^{2k_n})$, Markov’s inequality and (28) give

Since $n_q = O(n^q)$, $c_\varepsilon$ can be chosen so that the right-hand side is summable in $n$. Hence, by Borel–Cantelli, $\limsup_{n\to\infty}\|\mathbf{M}\|_{\mathrm{op}} \le2$, almost surely. Finally, the empirical spectral distribution converges almost surely to the semicircle law. Since the latter has positive mass in $(2-\varepsilon,2)$, almost surely $\mathbf{M}$ has an eigenvalue larger than $2-\varepsilon$ for all sufficiently large $n$. Thus $\liminf_{n\to\infty}\|\mathbf{M}\|_{\mathrm{op}} \ge2-\varepsilon$. Letting $\varepsilon\downarrow0$ concludes that $\|\mathbf{M}\|_{\mathrm{op}} \to2$, almost surely. ■

Phase Transition (Theorem 13)

Proof Similarly to Theorem 9, let $\mathbf{M} := \alpha_{p,q}^{-1}n^{-(p-2)/4}\Phi_q(\mathbf{Z})$ and $\mathbf{Q}(z) := (\mathbf{M}-z\mathbf{I}_{n_q})^{-1}$. By Proposition 18, $\lambda_{\max}(\mathbf{M}) \to2$ almost surely.

Suppose first that $\tau> 1$ and let $z \in\mathbb{C}\setminus[-2,2]$. Then,

By Theorem 9, one has $\mathbf{v}^{\odot q\mathsf{T}}\mathbf{Q}(z)\mathbf{v}^{\odot q}\to s_{\mathrm{sc}}(z)$. Since $1+\tau s_{\mathrm{sc}}(z)=0$ has a unique solution outside $[-2,2]$, that is $\rho(\tau):=\tau+\tau^{-1}$, $\mathbf{X}$ has a unique eigenvalue $\hat{\tau}$ converging to $\rho(\tau)>2$. Let $\Gamma$ be a small contour around $\rho(\tau)$. By Cauchy’s integral formula (Lemma 2) and Sherman–Morrison’s formula, one has

The integrand has a unique simple pole at $\rho(\tau)$, so the residue theorem yields

Now suppose that $\tau\le1$ and fix $\varepsilon>0$. For $z\ge2+\varepsilon$, one has $1+\tau s_{\mathrm{sc}}(z)\ge1+s_{\mathrm{sc}}(2+\varepsilon)>0$. By Theorem 9, $1+\tau\mathbf{v}^{\odot q\mathsf{T}}\mathbf{Q}(z)\mathbf{v}^{\odot q}$ is therefore bounded away from zero, hence the determinant identity (30) shows that $\mathbf{X}$ has no eigenvalue larger than $2+\varepsilon$. By Weyl’s inequality, $\lambda_{\max}(\mathbf{X})\ge\lambda_{\max}(\mathbf{M})\to2$; eventually, $\lambda_{\max}(\mathbf{X})\to2$. Finally, the subcritical finite-rank deformation result for Wigner matrices, together with the isotropic deterministic equivalent proved in Theorem 9, gives $|\hat{\mathbf{v}}^{\mathsf{T}}\mathbf{v}^{\odot q}|^2\to0$.

The results apply similarly to $\operatorname{Mat}_q(\mathbf{X})$, as $(\mathbf{U}\hat{\mathbf{v}})^{\mathsf{T}}(\mathbf{U}\mathbf{v}^{\odot q})=\hat{\mathbf{v}}^{\mathsf{T}}\mathbf{U}^{\mathsf{T}}\mathbf{U}\mathbf{v}^{\odot q}=\hat{\mathbf{v}}^{\mathsf{T}}\mathbf{v}^{\odot q}$. ■

Standard Gaussian Noise

This section extends the main text results to the spiked tensor model $\mathbf{X}:=\lambda\mathbf{v}^{\otimes p}+\mathbf{Z}$, where $\mathbf{Z}\sim\mathrm{SG}_{p,n}$ has independent standard Gaussian entries up to tensor symmetry:

Proof [Theorem 19, Corollary 20 and Theorem 21] The proof is identical in structure to the one of $\mathbf{Z}\sim\mathrm{GOTE}_{p,n}$. After multiplication by the symmetric basis weights and division by $\beta n^{p/4}$, a generic entry of $\mathbf{M}$ has variance $q!^2/(\beta^2n^{p/2})=n_q^{-1}(1+O(n^{-1}))$. The dependency classes are exactly the same, whence the negligibility of non-standard symmetries. Eventually, the Stein, Nash–Poincaré and Sherman–Morrison arguments apply similarly.

We report the theoretical and simulated spectral distribution of the symmetrically unfolded standard Gaussian model, $\mathbf{M}:=\beta^{-1}n^{-p/4}\Phi_q(\mathbf{Z})$, in Figure 5, with $n=100$. Moreover, we report the eigenvalues and overlaps in Figure 6, with $n\in\{20,50,100\}$, in a similar experimental setting to § 3.4. We see that the empirical eigenvalues are well estimated by the limiting spectral distribution. Both leading eigenvalue and overlap are increasingly well estimated by the theoretical ones beyond phase transition, as $\lambda$ grows.

The match between theory and simulations is slightly better under $\mathbf{Z}\sim\mathrm{SG}_{p,n}$ than $\mathbf{Z}\sim\mathrm{GOTE}_{p,n}$; see Figure 3 and Figure 4 for comparison. This may be due to the extra asymptotic approximations that are performed in the proof under the GOTE noise, notably $\#\mathcal{S}_{i\sqcup k}=p!\#\mathcal{S}_i\#\mathcal{S}_k/q!$, which holds only for generic indices.

Finite-Rank Extensions

Theorem 13 on phase transition extends to a finite number of orthogonal planted signals.

Proof. After the symmetric unfolding and normalization, the model is $\mathbf{M}+\sum_{k=1}^{r}\tau_k\mathbf{u}_k\mathbf{u}_k^{\top}$, with $\mathbf{u}_k=\mathbf{v}_k^{\otimes q}$. The vectors $\mathbf{u}_k$ are orthonormal because $\langle\mathbf{u}_k,\mathbf{u}_l\rangle=\langle\mathbf{v}_k,\mathbf{v}_l\rangle^q$. The deformation has fixed rank. The isotropic deterministic equivalent from Theorem 9 reduces its outlier equations to $1+\tau_k s_{\mathrm{sc}}(z)=0$ after applying Woodbury’s formula—a rank-$r$ version of Sherman–Morrison’s. Distinct supercritical strengths therefore give separated outliers at $\tau_k+\tau_k^{-1}$ and the residue computation is the same as in the rank-one proof, yielding $|\widehat{\mathbf{u}}_k^{\top}\mathbf{u}_k|^2\to(1-\tau_k^{-2})^{+}$ and vanishing cross-overlaps. The following result on subspace alignment follows from the convergence of the overlap matrix.

We report the theoretical and simulated spectral distribution and subspace overlap of the symmetrically unfolded $CP$ model (37), $\mathbf{M}:=\alpha_{p,q}^{-1}n^{-(p-2)/4}\Phi_q(\mathbf{X})$, in Figure 7, with $p=4$, $r=4$ and $n\in\{20,50,100\}$. We see that the $r$ largest eigenvalues escape the support of the limiting spectral density; they are therefore above the phase transition threshold. The theoretical subspace overlap curve undergoes multiple phase transitions, one for each spike. It increasingly well estimates the finite-sample simulations as $n$ and $\lambda$ grow.

The preceding subspace recovery result suggests an extension to other tensor models, such as the Tucker decomposition [60]:

where $\mathbf{H}\in\operatorname{Sym}^{p}(\mathbb{R}^{r})$, $\mathbf{v}_1,\ldots,\mathbf{v}_r\in\mathbb{S}^{n-1}$ are orthonormal and $\mathbf{Z}\sim\operatorname{GOTE}_{p,n}$. Such a model is notably relevant in applications where the signal is described by a low-dimensional

subspace rather than rank-one components [54, 41]. However, the resulting lifted subspace estimate—an $\binom{n+q-1}{q} \times\binom{r+q-1}{q}$ matrix—would not straightforwardly lead to an estimate of the subspace span($v_1,\ldots,v_r$). We therefore leave the Tucker extension for future research.

Recursive Unfolding

The main text analysis concerns a single unfolding of the tensor, resulting in overlap formulas with the lifted vector $v^{\boxtimes q} \in\mathbb{S}^{n^q-1}$. It is natural to ask whether the corresponding recursive unfolding procedure, introduced and analyzed in the asymmetric setting by Montanari and Richard [46], also achieves exact recovery in the genuinely symmetric model. The similarity between the two Gaussian tensor models—notably enabled by Lemma 17 on the negligible high-order symmetries—suggests that it should.

Proof [Idea of proof] The supercritical case can be understood by separating the signal component of the leading eigenvector from its orthogonal remainder. Write

with $u = v^{\boxtimes q}$, $g \perp u$ and $\lVert g\rVert= 1$ and set $\mathbf{G} = \operatorname{Mat}_1(g)$. By Theorem 13, $\rho_n^2 \to1-\tau^{-2} > 0$ whenever $\tau> 1$. Since $\operatorname{Mat}_1(u) = vv^{\boxtimes(q-1)\mathsf{T}}$, we consequently have

The first term in (42) is rank one with singular value $\rho_n$, bounded away from zero. Thus, if one could establish

then the Davis–Kahan/Wedin perturbation bound would imply $\lvert\hat{v}^{\mathsf{T}}v\rvert\to1$.

In the asymmetric model, Montanari and Richard [46] establish the analogue of (43) in the proof of their Theorem 5, with the sharper bound $\lVert\mathbf{G}'\rVert_{\mathrm{op}} = O(n^{-1/2})$, almost surely. We expect the same estimate to hold in the symmetric model, with a negligibility argument similar to Lemma 17. Establishing (43) rigorously would, however, require an operator-norm estimate for the matricized eigenvector remainder that is not provided by our present deterministic-equivalent analysis. We therefore leave this step as a conjecture.

The numerical experiments in Figure 8 support Conjecture 23.

Choice of Unfolding

We saw in § 4 that, if the dependencies induced by the high-order symmetries are spectrally negligible for the $n_k \times n_{p-k}$ symmetric unfolding, then the disproportionate random matrix results of Feldman [22] suggest the following sharp phase transition.

This conjecture is conditional on a negligibility result analogous to Lemma 17, showing that the contractions associated with indices $(i',k') \in\mathcal{I}_k \times\mathcal{I}_{p-k}$ such that $i' \sqcup k' = i \sqcup k$ are spectrally negligible. By studying the function $k \mapsto k!(p-k)!$, we see that the sharp constant in the phase-transition threshold for the short-side singular vector is minimized by an “as square as possible” unfolding. This is consistent with the recommendation of Montanari and Richard [46], while highlighting an important distinction with the asymmetric-noise setting of Ben Arous et al. [12]. More precisely, the result of Montanari and Richard [46] concerns partial recovery on the long side of the unfolding, for which the relevant threshold has a larger order of magnitude in the disproportionate regime, namely $n^{(p-k-1)/2}$. Thus, their recommendation should be interpreted with some caution in light of the distinction between the two singular-vector sides. Nevertheless, our results show that, for the short-side singular vector, their “as square as possible” prescription is indeed optimal at the level of the sharp constant (although not at the level of the polynomial order, which is the same for all unfoldings, namely $n^{(p-2)/4}$).

We compare this conjecture with numerical simulations in Figure 9, for $(p,k,n) = (3,1,500)$ and $(p,k,n) = (4,\{1,2\},100)$. The simulations match the conjectured overlaps. For $p = 4$, the “as square as possible” unfolding ($k = 2$) undergoes a phase transition before the “as large as possible” unfolding ($k = 1$), resulting in better partial recovery. This superiority is reinforced with recursive unfolding; see Figure 10 for numerical simulations.

From the computational standpoint, the comparison is somewhat different. The square unfolding involves extracting the leading $r$ eigenpairs of a symmetric $n^{p/2} \times n^{p/2}$ matrix, which costs $O(rn^{p/2}) = O(rn^p/(p/2)!^2)$ using an iterative eigensolver—up to the number of iterations required for convergence. A non-square unfolding involves extracting the leading $r$ left singular vectors of a $n_k \times n_{p-k}$ matrix. An iterative SVD requires $O(rn_kn_{p-k})$ arithmetic operations per iteration, giving $O(rn^p/(k!(p-k)!))$ up to the same iteration-count considerations. Thus, in principle both types of unfolding have the same polynomial computational complexity, but the leading combinatorial factor is minimized by the “as large as possible” unfolding ($k = 1$). Moreover, when $k = 1$, the left singular vectors directly provide an estimate of $v \in S^{n-1}$, without requiring a subsequent recursive unfolding step.

Overall, this suggests a tradeoff between statistical and computational efficiency. From a recovery standpoint, an “as square as possible” unfolding ($k = \lfloor p/2 \rfloor$) is preferable because it minimizes the conjectured phase-transition threshold, while recursive unfolding can then be used to recover the planted signal $v$. From a computational standpoint, an “as large as possible” unfolding ($k = 1$) is preferable: it has the smallest leading combinatorial factor and directly produces an estimator in the original signal space. To the best of our knowledge, this is the first analysis of the choice of unfolding to exhibit such a tradeoff between computational efficiency and estimation performance. It highlights a qualitative distinction between symmetric and asymmetric tensor noise: in the asymmetric setting, the recovery threshold is independent of the choice of unfolding [12] (Feldman and Donoho 2023). It is also reminiscent of the hierarchy of statistical–computational tradeoffs provided by Kikuchi matrices [63, 7].