Landscape analysis of an improved power method for tensor decomposition
Year: 2022

Abstract:
In this work, we consider the optimization formulation for symmetric tensor decomposition recently introduced in the Subspace Power Method (SPM) of Kileel and Pereira. Unlike popular alternative functionals for tensor decomposition, the SPM objective function has the desirable properties that its maximal value is known in advance, and its global optima are exactly the rank-1 components of the tensor when the input is sufficiently low-rank. We analyze the non-convex optimization landscape associated with the SPM objective. Our analysis accounts for working with noisy tensors. We derive quantitative bounds such that any second-order critical point with SPM objective value exceeding the bound must equal a tensor component in the noiseless case, and must approximate a tensor component in the noisy case. For decomposing tensors of size D×m, we obtain a near-global guarantee up to rank õ(D⌊m/2⌋) under a random tensor model, and a global guarantee up to rank O(D) assuming deterministic frame conditions. This implies that SPM with suitable initialization is a provable, efficient, robust algorithm for low-rank symmetric tensor decomposition. We conclude with numerics that show a practical preferability for using the SPM functional over a more established counterpart.

heading: 1 Introduction
text: From applied and computational mathematics [27] to machine learning [36] and multivariate statistics [29] to many-body quantum systems [41], high-dimensional data sets arise that are naturally organized into higher-order arrays. Frequently these arrays, known as hypermatrices or tensors, are decomposed into low-rank representations. In particular, the real symmetric CANDECOMP/PARAFAC (CP) decomposition [13] is often appropriate:
T = K∑ i=1 λia ⊗m i . (1)
Here, we are given T , a real symmetric tensor of size D×m. The goal is to expand T as a sum of K rank-1 terms, coming from scalar/vector pairs (λi, ai) ∈ R× RD. Importantly, the number of terms K must be minimal possible for the given tensor, in which case K is called the rank of the input T .
When m > 2 and K = O(Dm), CP decompositions are generically unique (up to permutation and scaling) by fundamental results in algebraic geometry [10]. An actionable interpretation [3] is that CP decomposition infers well-defined latent variables {(λi, ai) : i ∈ [K]} encoded by T . Indeed in learning applications, where symmetric tensors are formed from statistical moments (higher-order covariances) or multivariate derivatives (higher-order Hessians), CP decomposition has enabled parameter estimation for mixtures of Gaussians [20, 35], generalized linear models [34], shallow neural networks [19, 24, 42], deeper networks [17, 18, 30], hidden Markov models [5], among others.
Unfortunately, CP decomposition is NP-hard in the worst case [23]. In fact, it is believed to possess a computational-to-statistical gap [7], with efficient algorithms expected to exist for random tensors
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
only up to rank K = O(Dm/2). For this reason, and due to its non-convexity and high-dimensionality, CP decomposition has served theoretically as a key testing ground for better understanding mysteries of non-convex optimization landscapes. To date, this focus has been on the non-convex program
max x∈RD:∥x∥2=1
〈 T, x⊗m 〉 . (PM-P)
We label the problem PM-P, standing for Power Method Program, because projected gradient ascent applied to PM-P corresponds to the Shifted Symmetric Higher-Order Power Method of Kolda and Mayo [28]. Important analyses of PM-P include Ge and Ma’s [21] and the earlier [4] on overcomplete tensors, as well as [33] for low-rank tensors, and [31] which studied PM-P assuming the tensor components form a unit-norm tight frame with low incoherence.
In this paper, we perform an analysis of the non-convex optimization landscape associated with the Subspace Power Method (SPM) for computing symmetric tensor decompositions. The first and third authors introduced SPM in [25]. This method is based on the following non-convex program:
max x∈RD:∥x∥2=1 FA(x), (SPM-P)
where FA(x) := ∥PA(x⊗n)∥2F , n := ⌈m/2⌉, A := Span{a⊗n1 , . . . , a ⊗n K },
and PA : (RD)⊗n → A is orthogonal projection with respect to Frobenius inner product. Note that SPM-P is a particular polynomial optimization problem of degree 2n on the unit sphere.
There are at least two motivating reasons to study the optimization landscape of SPM-P. Firstly, it was observed in numerical experiments in [25] that the SPM algorithm is competitive within its applicable rank range of K = O(D⌊m/2⌋). It gave a roughly one-order of magnitude speed-up over the decomposition methods in [26] as implemented in Tensorlab [40], while matching the numerical stability of FOOBI [14]. Thus SPM is a practical algorithm. Secondly, from a theory standpoint, the program SPM-P has certain desirable properties which PM-P lacks. Specifically for an input tensor T = ∑K i=1 λia ⊗m i with rank K ≲ D ⌊m/2⌋ and Zariski-generic1 {(λi, ai)}Ki=1, SPM-P is such that:
• Each component ±ai is exactly a global maximum, and there are no other global maxima. • The globally maximal value is known in advance to be exactly 1. So the objective value
gives a certificate for global optimality, and non-global critical points can be discarded.
These properties were shown for SPM-P in [25], but both fail for PM-P (see Figure 1). Thus, SPM-P is more relevant theoretically than PM-P as a test problem for non-convex CP decomposition.
Prior theory. In [25], it is proven that projected gradient ascent applied to SPM-P, initialized at almost all starting points with a constant explicitly-bounded step size, must converge to a secondorder critical point of SPM-P at a power rate or faster. However this left open the possible existence
1Zariski-generic means that the failure set can be described by the vanishing non-zero polynomials [22], so in particular, has Lebesgue measure 0.
of spurious second-order critical points, i.e., second-order points with reasonably high objective value that are not global maxima (unequal to, and possibly distant from, each CP component ±ai). Such critical points could pose trouble for the successful optimization of SPM-P. Furthermore all theory for SPM-P in [25] was restricted to the clean case: that is, when the input tensor T is exactly of a sufficiently low rank K. The analysis for PM-P in [21, 31, 33] also assume noiseless inputs. However, tensors arising in practice are invariably noisy due, e.g., to sampling or measurement errors.
Main contributions. We perform a landscape analysis of SPM-P by characterizing all second-order critical points, using suitable assumptions on a1, . . . , aK . Under deterministic frame conditions on a1, . . . aK , which are satisfied by mutually incoherent vectors, near-orthonormal systems, and random ensembles of size K = O(D), Theorem 7 shows that all second-order critical points of SPM-P coincide exactly with±ai. Theorem 16 shows the same result for overcomplete random ensembles of size K = õ(D⌊m/2⌋), however requiring an additional superlevel set condition to exclude maximizers with vanishing objective values. Both results extend to noisy tensor decomposition, where SPM is applied to a perturbation T̂ ≈ T . In this setting, second-order critical points with objective values exceeding O(∥T̂ − T∥F ) are O(∥T̂ − T∥F )-near to one of the components ±ai. We also show in Lemma 9 that spurious local maximizers (with lower objective values) do exist in the noisy case.
The results imply a clear separation of the functional landscape between near-global maximizers, with objective values close to 1, and spurious local maximizers with small objective value. Hence, the SPM objective can be used to validate a final iterate of projected gradient ascent in the noisy case. In Theorem 18, we combine our landscape analysis with bounds on error propagation incurred during SPM’s deflation steps. This gives guarantees for end-to-end tensor decomposition using SPM.
Lastly, we expose the relation between PM-P and SPM-P. Specifically, SPM-P can be expressed as PM-P with the appropriate insertion of the inverse of the Grammian (Gn)ij := ⟨ai, aj⟩n. The resulting de-biasing effect on the local maximizers with respect to the components ±ai (cf. Figure 1) is responsible for many advantages of SPM-P. Along the way, we state a conjecture about the minimal eigenvalue of Gn when ai are i.i.d. uniform on the sphere, which may be of independent interest.

heading: 2 Notation
text: Vectors and matrices. When x is a vector, ∥x∥p is the ℓp-norm (p ∈ R≥1∪{∞}). For x, y ∈ RD, the entrywise (or Hadamard) product is x⊙y ∈ RD, and the entrywise power is x⊙s := x⊙. . .⊙x ∈ RD (s ∈ N). When M is a matrix, ∥M∥2 is the spectral norm. If M is real symmetric, µj(M) is the eigenvalue of M that is the j-th largest in absolute value. We denote the identity by IdD ∈ RD×D.
Tensors. A real tensor of length D and order m is an array of size D ×D × . . .×D (m times) of real numbers. Write T mD := ( RD )⊗m ∼= RDm for the space of tensors of size D×m. Meanwhile, Sym(T mD ) ⊆ T mD is the subspace of symmetric tensors (i.e., tensors unchanged by any permutation of indices). The Frobenius inner product and norm are denoted by ⟨·, ·⟩ and ∥·∥F , respectively. Given any linear subspace of tensors A ⊆ T mD , let PA : T mD → A denote the orthogonal projector onto A with respect to ⟨·, ·⟩. In the case A = Sym(T mD ), the projector PSym(T mD ) is the symmetrization operator, Sym : T mD → Sym(T mD ). Given T ∈ T m1 D and S ∈ T m2 D , the tensor (or outer) product is T ⊗ S ∈ T m1+m2D , defined by (T ⊗ S)i1,...,im1+m2 := Ti1,...,im1Sim1+1,...,im2 . For T ∈ T mD and s ∈ N, the tensor power is T⊗s := T ⊗ . . . ⊗ T ∈ T smD . For T ∈ T m1 D , S ∈ T m2D with m1 ≥ m2, the contraction T · S ∈ T m1−m2 D is defined by (T · S)i1,...,im1−m2 :=∑
j1,...,jm2 Ti1,...,im1−m2 ,j1,...,jm2Sj1,...,jm2 . Let Reshape(T, [d1, . . . , dℓ]) be the function that re-
shapes the tensor T to have dimensions d1, . . . , dℓ, as in corresponding Matlab/NumPy commands.
Other. The unit sphere in RD is SD−1, and Unif(SD−1) is the associated uniform probability distribution. Given a function f : RD → R, the Euclidean gradient and Hessian matrix at x ∈ RD are ∇f(x) ∈ RD and ∇2f(x) ∈ Sym(T 2D). The Riemannian gradient and Hessian with respect to SD−1 at x ∈ SD−1 are ∇SD−1f(x) and ∇2SD−1f(x) (see [1]). Write Span for linear span, [K] := {1, . . . ,K}, and |A| for the cardinality of a finite set A. Lastly, we use asymptotic notation freely.

heading: 3 Symmetric tensor decomposition via Subspace Power Method
text: In this section, we outline the tensor decomposition method SPM of [25], and provide basic insights on the program SPM-P. Throughout we assume that m ≥ 3 is an integer and define n := ⌈m/2⌉.
SPM algorithm. The input is a tensor T̂ ∈ Sym(T mD ), with the promise that T̂ ≈ T = ∑K i=1 λia ⊗m i
for {(λi, ai)}Ki=1 Zariski-generic and K ≤ ( D+n−1
n
) −D (if m is even) and K ≤ Dn (if m is odd).
As a first step, SPM obtains the orthogonal projector PÂ : T mD → Â that projects onto the column span of Reshape(T̂ , [Dn, Dm−n]), by using matrix singular value decomposition. Provided that T̂ ≈ T , the associated subspace approximation error Â ≈ A = span{a⊗n1 , . . . , a ⊗n K } defined by
∆A := ∥∥PA − PÂ∥∥F→F = sup
T∈T nD , ∥T∥F=1 ∥∥PA(T )− PÂ(T )∥∥F , (2) can be bounded as follows. (Note that by [9, Lem. 2.3], we know ∆A ≤ 1 a priori.) Lemma 1 (Error in subspace). Let m ≥ 3, n = ⌈m2 ⌉, T ∈ Sym(T m D ) and assume that M := Reshape(T, [Dn, Dm−n]) has exactly K nonzero singular values σ1(M) ≥ . . . ≥ σK(M) > 0. Let T̂ ∈ Sym(T mD ), M̂ := Reshape(T, [Dn, Dm−n]). Assume ∆M := ∥M − M̂∥2 < σK(M). Then∥∥∥PIm(M) − PImK(M̂)∥∥∥2 ≤ ∆MσK(M)−∆M , (3) where Im(M) ⊆ RDn denotes the image of M and ImK(M̂) ⊆ RD n
denotes the subspace spanned by the K leading left singular vectors of M̂ . In particular, if T = ∑K i=1 λia ⊗2n i , A = span{a⊗ni : i ∈ [K]}, dim(A) = K, and Â is the subspace spanned by K leading tensorized left singular vectors of M̂ , the right-hand side of (3) upper-bounds ∆A.
Remark 2. If T = ∑K
i=1 λia ⊗m i , one coefficient λi is small and the vectors {ai : i ∈ [K]} are
not too correlated, then the flattened tensor Reshape(T, [Dn, Dm−n]) has a small eigenvalue. This makes estimating the corresponding eigenvector sensitive to noise. See Remark S.41 in the appendix.
Given Â, SPM seeks one tensor component ai by solving the noisy variant of SPM-P defined by
max x∈SD−1
FÂ(x), where FÂ(x) := ∥∥PÂ(x⊗n)∥∥2F . (nSPM-P)
Starting from a random initial point x0 ∼ Unif(SD−1), the projected gradient ascent iteration
x← x+ γPÂ(x
⊗n) · x⊗(n−1)
∥x+ γPÂ(x⊗n) · x⊗(n−1)∥2 , (4)
with a constant step-size γ, is guaranteed to converge to a second-order critical point of nSPM-P almost surely by [25]. Here we require that the step-size γ is less than an explicit upper bound given in [25]. Denoting by âi the final iterate obtained by SPM, we accept the candidate approximate tensor component âi if FÂ(âi) is large enough; otherwise we draw a new starting point x0 and re-run (4).
Next given âi, SPM evaluates a deflation formula based on Wedderburn rank reduction [12] from matrix algebra to compute the corresponding weight λ̂i. Then, we update the tensor T̂ ← T̂ − λ̂iâ⊗mi . To finish the tensor decomposition, SPM performs the projected gradient ascent and deflation steps K times to compute all of the tensor components and weights {(λ̂i, âi)}Ki=1.
Preparatory material about nSPM-P. The goal of this paper is to show that second-order critical points of nSPM-P with reasonable function value must be near the global maximizers±a1, . . . ,±aK of SPM-P, under suitable incoherence assumptions on the rank-one components a1, . . . , aK . Naturally, the optimality conditions for nSPM-P play an important part in this analysis. Proposition 3 (Optimality conditions). Let x ∈ SD−1 be first and second-order critical for nSPM-P. Then for each z ∈ SD−1 with z ⊥ x, we have
PÂ(x n) · xn−1 = FÂ(x)x, (5) FÂ(x) ≥ n∥PÂ(x n−1z)∥2F + (n− 1)⟨PÂ(x
n), xn−2z2⟩. (6) Furthermore, for any y ∈ SD−1 we have
FÂ(x) ≥ n ∥∥PÂ(xn−1y)∥∥2F + (n− 1)⟨PÂ(xn), xn−2y2⟩ − 2(n− 1)FÂ(x)⟨x, y⟩2. (7)
In the analysis later, we make frequent use of expressing the objective FA(x) using the Gram matrix
Gn ∈ Sym(T 2D) defined by (Gn)ij := ⟨a⊗ni , a ⊗n j ⟩ = ⟨ai, aj⟩ n. (8)
Under linear independence of the tensors a⊗n1 , . . . , a ⊗n K , which is implied by our assumptions made later, the inverse G−1n exists and the noiseless program SPM-P can be expressed as follows. Lemma 4. Let A := [a1| . . . |aK ] ∈ RD×K and {a⊗ni : i ∈ [K]} be linearly independent. We have
FA(x) = ∥∥PA(x⊗n)∥∥2F = ((A⊤x)⊙n)⊤ G−1n ((A⊤x)⊙n) . (9)
Lemma 4 exposes the relation between SPM-P and PM-P. While PM-P can be rewritten ⟨T, x⊗m⟩ = ⟨(A⊤x)⊙n, (A⊤x)⊙n⟩ if m is even, SPM-P takes into account correlations among the tensors a⊗n1 , . . . , a ⊗n K and inserts the Grammian G −1 n into (9). Consequently, correlations among the tensor components are considered in SPM-P, without any a priori knowledge of the tensors a⊗n1 , . . . , a ⊗n K .
In the special case of orthonormal systems, or more generally systems that resemble equiangular tight frames [16, 37], the SPM-P and PM-P objectives coincide up to shift and scaling. Lemma 5. Assume there exist ρ ∈ (−1, 1) \ { −1K−1} and M ∈ R such that ⟨ai, aj⟩
n = ρ for all i ̸= j and ∑ i∈[K]⟨x, ai⟩n = M for all x ∈ SD−1. Denote A = [a1| . . . |aK ] ∈ RD×K . Then
FA(x) = (1− ρ)−1∥A⊤x∥2n2n − ( (1− ρ)2 +Kρ(1− ρ) )−1 ρM2. (10)

heading: 4 Main results
text: In this section, we present the main results about local maximizers of the nSPM-P program. Section 4.1 is tailored to low-rank tensor models with K = O(D) components that satisfy certain deterministic frame conditions. Section 4.2 then considers the overcomplete case K = õ(D⌊m/2⌋) in an average case scenario, where a1, . . . , aK are modeled as independent copies of an isotropic random vector.

heading: 4.1 Low-rank tensors under deterministic frame conditions
text: Motivated by frame constants in frame theory [11], we measure the incoherence of the ensemble a1, . . . , aK by scalars ρs ∈ R≥0, which are defined via
ρs := sup x∈SD−1 K∑ i=1 |⟨x, ai⟩|s − 1. (11)
They satisfy the order relation ρs ≤ ρs′ for s′ ≤ s, due to ∥ai∥2 = 1, and can be related to extremal eigenvalues of Grammians Gs and G⌊s/2⌋ as shown in the following result.
Lemma 6. Let {ai : i ∈ [K]} ⊆ SD−1 and (Gs)ij := ⟨ai, aj⟩s for s ∈ N. Then
1− ρs ≤ µK(Gs) ≤ µ1(Gs) ≤ 1 + ρs ≤ µ1(G⌊s/2⌋). (12)
The characterization in Lemma 6 allows to compute bounds for ρs for low-rank tensors with mutually incoherent components or rank-O(D) tensor with random components. We provide details on this in Remark 8 below, but first state the main guarantee about local maximizers using ρ2 and ρn. Theorem 7 (Main deterministic result). Let {ai : i ∈ [K]} ⊆ SD−1 and A = Span{ani : i ∈ [K]}. Let Â ⊆ Sym(T nD ) be a perturbation of A with ∆A = ∥PA − PÂ∥F→F . Let
τ := 1
6 − n2ρ2 − (n2 + n)ρn and ∆0 :=
2τ
2 + 4τ + 3n2 . (13)
Then, if ∆A < ∆0, the program nSPM-P has exactly 2K second-order critical points in the superlevel set where
FÂ(x) ≥ 2 + 2τ + 3n2
2τ ∆A. (14)
Each of these critical points is a strict local maximizer for FÂ. Further for each such point x ∗, there exists unique i ∈ [K] and s ∈ {−1, 1} such that
∥x∗ − sai∥22 ≤ 2∆A n . (15)
In the noiseless case (∆A = 0), if τ ≥ 0 then there are precisely 2K second-order critical points of SPM-P with positive functional value, and they are the global maximizers sai (i ∈ [K], s ∈ {−1, 1}). Remark 8. Using Lemma 6, we identify two situations where ρ2 and ρn can be bounded from above.
1. Mutually incoherent ensembles. Let a1, . . . , aK have mutual incoherence ρ := max i ̸=j |⟨ai, aj⟩|. Using Gershgorin’s circle theorem, we obtain
ρs ≤ µ1(G⌊s/2⌋)− 1 ≤ max i∈[K] ∑ j ̸=i |⟨ai, aj⟩|⌊s/2⌋ ≤ (K − 1)ρ⌊s/2⌋ for any s ∈ N≥2,
which implies that the conditions of Theorem 7 are saisfied if Kρ is sufficiently small. This setting is comparable to the analysis for PM-P in [33]. Moreover, if K = D and a1, . . . , aK are mutually orthogonal, then ρs = 0 for each s ≥ 2. Therefore Theorem 7 holds with τ = 1/6 for orthogonally decomposable tensors [6].
2. Low-rank random ensembles. Let a1, . . . , aK be independent copies of an isotropic unitnorm random vector a with sub-Gaussian norm O(1/ √ D) (e.g., ai ∼ Unif(SD−1)). With
high probability, the ensemble can achieve arbitrarily small ρ2, provided that K ≤ CD for a sufficiently small constant C > 0. The proof of this fact relies on A = [a1| . . . |aK ] satisfying, with high probability, the so-called (K, δ) restricted isometry property [38, Thm. 5.65] as defined in Definition 10 in the next section. We note that despite requiring conditions milder than those in the unit-norm tight frame analysis for PM-P in [31], we achieve a comparable scaling of K = O(D).
In the noiseless case where ∆A = 0, Theorem 7 shows that all local maximizers coincide with global maximizers ±a1, . . . ,±aK , provided the tensor components are sufficiently incoherent to ensure τ ≥ 0. In the noisy case, all local maximizers with objective values FÂ(x) ≥ C(n, τ)∆A are O(∆A)-close to the global optimizers of the noiseless objective FA, where the constant C(n, τ) increases as the incoherence of the vectors a1, . . . , aK shrinks. Unfortunately, the presence of spurious local maximizers of FÂ with small objective values cannot be avoided under a deterministic noise model, as the next result shows. Lemma 9. Let δ ∈ (0, 1), a ∈ SD−1 and A = span{a⊗n}. Then there exists a subspace Â ⊆ Sym(T nD ) with dim(Â) = 1 and ∥PA − PÂ∥F→F = δ such that nSPM-P possesses a strict local maximizer of objective value exactly δ2.

heading: 4.2 Average case analysis of overcomplete tensors
text: The overcomplete case with K = õ(D⌊m/2⌋) falls outside the range of Theorem 7, because τ in (13) becomes negative when K ≫ D. Instead, our analysis for the overcomplete case relies on A = [a1| . . . |aK ] ∈ RD×K obeying the (p, δ)-restricted isometry property (RIP) for p = O(D/ log(K)). Definition 10. Let A ∈ RD×K , 1 ≤ p ≤ K be an integer, and δ ∈ (0, 1). We say that A is (p, δ)-RIP if every D × p submatrix Ap of A satisfies ∥A⊤p Ap − Idp∥2 ≤ δ.
A consequence of the RIP, which is particularly useful in the analysis of overcomplete tensor models, is that the correlation coefficients {⟨ai, x⟩ : i ∈ [K]} can naturally be split into two groups. Lemma 11 (RIP-induce partitioning of correlation coefficients). Suppose that A = [a1| . . . |aK ] ∈ RD×K satisfies the (p, δ)-RIP for p = ⌈cδD/ log(K)⌉. Let c̃δ := (1+ δ)/cδ . Then for all x ∈ SD−1 there is a subset of indices I(x) ⊆ [K] with cardinality p such that
1− δ ≤ ∑
i∈I(x)
⟨ai, x⟩2 ≤ 1 + δ and ⟨ai, x⟩2 ≤ c̃δ log(K)
D for i ̸∈ I(x), (16)
Let us now collect all assumptions needed to analyze the overcomplete case.
Definition 12. Let A := [a1| . . . |aK ]. We require the following assumptions:
A1 There exists cδ > 0, depending only on δ, such that A is (⌈cδD/ log(K)⌉, δ)-RIP.
A2 There exists c1 > 0, independent of K,D, such that maxi,j:i ̸=j⟨ai, aj⟩2 ≤ c1 log(K)/D. A3 There exists c2 > 0, independent of K,D, such that ∥∥G−1n ∥∥2 ≤ c2.
To the best of our knowledge, there are no known deterministically constructed systems of size K = õ(D⌊m/2⌋) which satisfy A1 - A3 for small δ ≪ 1 and constants cδ, c1, c2 that do not depend on K,D. However, by modeling components via a sufficiently-spread random vector, i.e., by considering an average case scenario with ai ∼ Unif(SD−1), they hold with high probability. Proposition 13 (A1 - A3 for random ensembles). Let a ∼ Unif(SD−1) and assume a1, . . . , aK are K independent copies of a, where log(K) = o(D). Fix an arbitrary constant δ ∈ (0, 1). There exist constants C > 0 and D0 ∈ N depending only on δ such that for all D ≥ D0, and with probability at least 1−K−1 − 2 exp(−Cδ2D), conditions A1 and A2 (with c1 ≤ C) hold. Furthermore, if n = 2 and K = o(D2), A3 holds with probability 1− C(eD/ √ K)−C √ K for c2 ≤ C.
Remark 14 (A3 when n > 2). Following [2, 15], we give a self-contained proof for A3 when n = 2 in the appendix. We are currently not able to extend the proof to n > 2, because some technical tools such as the Hanson-Wright inequality [32] and an extension of [2, Thm. 3.3], which ensures the RIP for matrices whose columns have sub-exponential tails, have not yet been fully developed for random vectors with dependent entries and heavier tails (so-called α-sub-exponential tails or sub-Weibull tails). However we strongly believe that A3 holds for n > 2 and K = o(Dn). For now, we formulate this as a conjecture, supplemented with numerical evidence presented in Figure 2. We also add that the conjecture, once proven, would complement recent advances on the well-posedness of random tensors with fewer statistical dependencies among the components [8, 39].
Conjecture 15 (Grammians of independent symmetric rank-one tensors). Let a1, . . . , aK be independent copies of the random vector a ∼ Unif(SD−1) and fix ϵ > 0 arbitrarily. Then there exists some constant κn > 0 and an increasing function γn : (0, κn)→ (0, 1), both depending only on n, such that γn(κ)→ 1 as κ→ 0, and if K ≤ κDn for some κ < κn we have
P( ∥∥G−1n ∥∥2 ≥ γn(κ)− ϵ)→ 1 as D →∞. (17)
In particular, if K = o(Dn) then ∥G−1n ∥2 → 1 as D →∞.
We now present our main theorem about local maximizers of nSPM-P in the overcomplete case. Theorem 16 (Main random overcomplete result). Let K,D ∈ N, define εK := K logn(K)/Dn, and suppose that limD→∞ εK = 0. Assume a1, . . . , aK ∈ RD satisfy A1 - A3 for some δ, c1, c2 > 0. Then there exist δ0, depending only on n, c1 and c2, and D0,∆0, C, which depend additionally on cδ , such that if δ < δ0, D > D0, and ∆A ≤ ∆0, the program nSPM-P has exactly 2K second-order critical points in the superlevel set{
x ∈ SD−1 : FÂ(x) ≥ CεK + 5∆A } . (18)
Each of these critical points is a strict local maximizer for FÂ. Further for each such point x ∗, there exists unique i ∈ [K] and s ∈ {−1, 1} such that
∥x∗ − sai∥22 ≤ 2∆A n . (19)
In particular, in the noiseless case (∆A = 0), the second-order critical points in the superlevel set (18) are exactly the global maximizers sai (i ∈ [K], s ∈ {−1, 1}).
By Theorem 16, all non-degenerate local maximizers in the superlevel set (18) are close to global maximizers of the noiseless objective FA. In the noiseless case ∆A = 0, they coincide with global maximizers and the right hand side in (18) tends to 0 as D →∞. Hence, we have a clear separation between the global maximizers ±a1, . . . ,±aK and degenerate local maximizers with vanishing objective value, so that the objective value acts as a certificate for the validity of an identified maximizer.
Theorem 16 is not a global guarantee because a random starting point x0 ∼ Unif(SD−1), used for starting the projected gradient ascent iteration (4), is, with high probability, not contained in (18).
Remark 17 (Objective value at random initialization). For a random sample x0 ∼ Unif(SD−1) in the clean case ∆A = 0, we empirically observe the objective value FA(x0) ≈ CK/Dn as illustrated in Figure 3. As such, we conjecture Ex0,a1,...,ak∼Unif(SD−1)[FA(x0)] = O( K Dn ). Comparing the objective to the level set condition in Theorem 16,{ x ∈ SD−1 : FA(x) ≥ C K logn(K)
Dn + 5∆A
} , (20)
a random starting point x0 therefore falls short of satisfying (20) by a logn(K)-factor only. In this sense, Theorem 16 furnishes a “near-global" guarantee for nSPM-P in the random overcomplete case.

heading: 4.3 End-to-end tensor decomposition
text: By combining our landscape analyses with bounds for error propagation during deflation, we obtain a theorem about end-to-end tensor decomposition using the SPM algorithm. That is, under the conditions of Theorem 7 or 16, a tweaking of SPM (Algorithm 1 in the appendix) recovers the entire CP decomposition exactly in the noiseless case (T̂ = T ). In the noisy regime, it obtains an approximate CP decomposition, and we bound the error in terms of ∆A. Due to space constraints, we leave precise descriptions of Algorithm 1 and our deflation bounds to the supplementary material.
Theorem 18 (Main result on end-to-end tensor decomposition). Let T = ∑K
i=1 λia ⊗m i ∈ Sym(T mD )
and M := Reshape(T, [Dn, Dm−n]). Let σ1(M) ≥ . . . ≥ σK(M) be the singular values of M , and assume σK(M) > 0. For other tensor T̂ ∈ Sym(T mD ), let M̂ = Reshape(T̂ , [Dn, Dm−n]), assume ∆M := ∥M − M̂∥2 < 12σK(M) and let ∆̂A = ∆M σK(M)−∆M . Suppose that T satisfies the assumptions of either Theorem 7 or Theorem 16, define ∆0 as in the corresponding theorem statement and let ℓ(∆A) be the corresponding level set threshold.2 Then there exist constants C1, C2, not depending on T̂ or ∆M , such that if we define ∆̃A := C1∆̂A+C2 √ ∆̂A the following holds. Assume
2In both theorem statements, the level set threshold depends on ∆A.
∆̃A < ∆0 and when we apply Algorithm 1 to T̂ each run of projected gradient ascent converges to a point with functional value at least ℓ(∆̃A). Then, letting (âi, λ̂i)i∈[K] be the output of Algorithm 1 applied to T̂ , there exist a permutation π ∈ Perm(K) and signs s1, . . . , sK ∈ {−1,+1} with
∥siaπ(i) − âi∥2 ≤ √ 2∆̂A n and ∣∣∣∣ smiλπ(i) − 1λ̂i ∣∣∣∣ ≤ 2 √ m/n σK(M) √ ∆̂A +
4
σK(M) ∆̂A ∀ i ∈ [K].
In particular, in the noiseless case (∆M = 0), Algorithm 1 returns the exact CP decomposition of T .
We conclude that the Subspace Power Method, with an initialization scheme for nSPM-P that gives x ∈ RD where FÂ(x) > ℓ(∆A), is a guaranteed algorithm for low-rank tensor decomposition.

heading: 5 Numerical experiments
text: Here we present numerical experiments that corroborate the theoretical findings of Section 4. We illustrate that SPM identifies exact tensor components in the noiseless case (in contrast to PM), and that SPM behaves robustly in noisy tensor decomposition. We use the implementation of SPM of the first and third authors, available at https://github.com/joaompereira/SPM, which is licensed under the MIT license. All of our experiments presented below may be conducted on a standard laptop computer within a few hours. For further numerical experiments, we refer the reader to [25], where SPM was tested in a variety of other scenarios, justifying it as a possible replacement for state-of-the-art symmetric tensor decomposition methods such as FOOBI [14] or Tensorlab [40].
Global optimizers and noise robustness. In the first set of experiments, we are interested in the recovery of individual tensor components ai for different D, K and m from noisy approximately low rank tensors. With m = 2n as the tensor order, we create noiseless tensors with ai ∼ Unif(SD−1)) as the tensor components and λi = √ Dm/K λ̃i as the tensor weights, where λ̃i ∼ Unif([1/2, 2]). This way, the variance of each entry of the tensor is about 1. We construct noisy tensors T̂ ≈ T by adding independent copies of ϵ ∼ N (0,m!σ2) to each entry of the tensor and then project onto
Sym(T mD ). The entries with all distint indices of the projected noise tensor have variance σ2; the variance of the remaining entries is also a multiple of σ2 and the number of such entries is o(Dm).
After constructing the tensor, we sample 10 points from Unif(SD−1) as initializations for projected gradient ascent applied to PM-P and nSPM-P. We iterate until convergence and record the distance between the final iterate and the closest ±ai among the tensor components. Averages and standard deviations (represented by shaded areas) over 1000 distances (10 distances per tensor, 100 tensors) are depicted in Figures 4a and 4b. In Figure 4a, we use noiseless tensors with σ = 0, while in Figure 4b σ ranges from 10−5 to 1, and the plots for 4th and 6th order tensors are superposed. We set D = 20,K = 100 for the 4th order tensor (n = 2) and D = 10,K = 100 for the 6th order (n = 3).
Figure 4a illustrates that projected gradient ascent applied to nSPM-P always converges to the global maxima in the noiseless case (up to numerical error), provided the scaling is K = O(D⌊m/2⌋) = O(D2) and the constant adheres to the rank constraints pointed out in [25, Prop. 4.2]. This is in agreement with Theorem 16, which shows local maxima with sufficiently large objective have to coincide with global maximizers. In the noisy case, Figure 4b illustrates the robustness of nSPM-P, giving an error of O(σ). In contrast to nSPM-P, the PM-P objective suffers from a large bias and does not recover exact tensor components. The effect of the bias even dominates errors induced by moderately-sized entrywise noise.
Deflation with SPM. In the second experiment, we recover complete tensors by using the deflation procedure described in [25]. We are mostly interested in the noisy case, since deflation with exact ai’s, as identified by nSPM-P in the noiseless case, does not induce any additional error.
We test 4th order tensors (n = 2) with D = 40,K = 200 and 6th order tensors (n = 3) with D = 15,K = 200, vary the noise level σ, and construct random tensors as in the previous experiment. Figure 4c plots the average error over 100 repetitions of successively recovered tensor components, where the x-axis ranges from the first recovered component at 1 to the last component at index 300. The figure illustrates that nSPM-P combined with modified deflation allows for recovering all tensor components up to an error of O(σ). Surprisingly, we do not observe error propagation, despite the fact that noisy recovered tensor components are being used within each deflation step.

heading: 6 Conclusion
text: We presented a quantitative picture for the optimization landscape of a recent formulation [25] of the non-convex, high-dimensional problem of symmetric tensor decomposition. We identified different assumptions on the tensor components a1, . . . , aK and bounds on the rank K so that all second-order critical points of the optimization problem SPM-P with sufficiently high functional value must equal one of the input tensor’s CP components. In Theorem 7 the assumptions were deterministic frame and low-rank conditions, while in Theorem 16 the hypotheses were random components and an overcomplete rank. Our proofs accommodated noise in the input tensor’s entries, and we obtained robust results for only by analyzing the program nSPM-P. Our analysis has algorithmic implications. As the Subspace Power Method of [25] is guaranteed to converge to second-order critical points, by combining with analysis of deflation, it follows that SPM (with sufficient initialization) is provable under our assumptions. In Theorem 18 we gave guarantees for end-to-end decomposition using SPM.
Compared to the usual power method functional, the novelty of the SPM functional is the de-biasing role played by the inverse of a Grammian matrix recording correlations between rank-1 tensors (recall Eq. (8)). This Grammian matrix is responsible for many of the SPM functional’s desirable properties, but it also complicated our analysis. We showed that there are theoretical and numerical advantages in using the SPM functional for tensor decomposition over its usual power method counterpart.
This paper suggests several directions for future research:
• What are the average-case properties of the Grammian matrix Gn? We formulated Conjecture 15 about the minimal eigenvalue of Gn. How about the other eigenvalues?
• In the random overcomplete setting, if we assume no noise (∆A = 0) then can we dispense with the superlevel condition in Theorem 16? This would give a fully global guarantee.
• Why do we see no error propagation when using deflation and sequential solves of nSPM-P for CP decomposition? Can errors accumulate if we choose the noise deterministically?
Acknowledgements and funding disclosure
We thank Massimo Fornasier for suggesting this collaboration in the beginning, after a fortuitous exchange in Oaxaca, Mexico at the Computational Harmonic Analysis and Data Science workshop (October 2019). J.K. acknowledges partial support from start-up grants provided by the College of Natural Sciences and Oden Institute for Computational Engineering and Sciences at UT Austin.

Review summaries:
This paper considers the subspace power method (SPM) for symmetric tensors. Specifically, this paper analyzes the non-convex optimization landscape associated with the SPM objective. The main results are two theorems, theorem 6 and theorem 15. Under some assumptions, it is shown in theorem 6  that for low rank case, nSPM-P has 2K local maximizers which are good approximations of all components. For the overcomplete tensors, similar conclusion is given. Numerical experiments is presented to validate the theoretical findings.

This paper studies the non-convex landscape of the SPM-P optimization objective (in page 2) for tensor decomposition. Under both deterministic and random models of the components, the work shows that all second-order critical points coincide with the ground truth vectors (up to a sign difference) with different levels of overcompleteness. The result can also be generalized to the noisy cases under certain assumptions. 

In this paper, the authors consider the optimization formulation for symmetric tensor decomposition introduced in the Subspace Power Method (SPM). They analyze the non-convex optimization landscape associated with the SPM objective. The analysis focuses on the noisy tensors, and implies SPM with suitable initialization is a provable, efficient, robust algorithm for low-rank symmetric tensor decomposition. Numerical results coincide with the theoretical analysis.  The contributions is to analyze the landscape of SPM-P by characterizing all second-order critical points.