Tensor principal component analysis via sum-of-squares proofs

Samuel B. Hopkins, Jonathan Shi, David Steurer

Introduction

Principal component analysis (pca), the process of identifying a direction of largest possible variance from a matrix of pairwise correlations, is among the most basic tools for data analysis in a wide range of disciplines. In recent years, variants of pca have been proposed that promise to give better statistical guarantees for many applications. These variants include restricting directions to the nonnegative orthant (nonnegative matrix factorization) or to directions that are sparse linear combinations of a fixed basis (sparse pca). Often we have access to not only pairwise but also higher-order correlations. In this case, an analog of pca is to find a direction with largest possible third moment or other higher-order moment (higher-order pca or tensor pca).

All of these variants of pca share that the underlying optimization problem is NP-hard for general instances (often even if we allow approximation), whereas vanilla pca boils down to an efficient eigenvector computation for the input matrix. However, these hardness result are not predictive in statistical settings where inputs are drawn from particular families of distributions. Here efficient algorithm can often achieve much stronger guarantees than for general instances. Understanding the power and limitations of efficient algorithms for statistical models of NP-hard optimization problems is typically very challenging: it is not clear what kind of algorithms can exploit the additional structure afforded by statistical instances, but, at the same time, there are very few tools for reasoning about the computational complexity of statistical / average-case problems. (See [BR13] and [BKS13] for discussions about the computational complexity of statistical models for sparse pca and random constraint satisfaction problems.)

We study a statistical model for the tensor principal component analysis problem introduced by [MR14] through the lens of a meta-algorithm called the sum-of-squares method, based on semidefinite programming. This method can capture a wide range of algorithmic techniques including linear programming and spectral algorithms. We show that this method can exploit the structure of statistical tensor pca instances in non-trivial ways and achieves guarantees that improve over the previous ones. On the other hand, we show that those guarantees are nearly tight if we restrict the complexity of the sum-of-squares meta-algorithm at a particular level. This result rules out better guarantees for a fairly wide range of potential algorithms. Finally, we develop techniques to turn algorithms based on the sum-of-squares meta-algorithm into algorithms that are truly efficient (and even easy to implement).

Montanari and Richard propose the following statistical modelMontanari and Richard use a different normalization for the signal-to-noise ratio. Using their notation, β≈τ/n\beta\approx\tau/\sqrt{n}. for tensor pca.

Montanari and Richard show that when τ⩽o(n)\tau\leqslant o(\sqrt{n}) Problem 1.1 becomes information-theoretically unsolvable, while for τ⩾ω(n)\tau\geqslant\omega(\sqrt{n}) the maximum likelihood estimator (MLE) recovers v′v^{\prime} with ⟨v,v′⟩⩾1−o(1)\langle v,v^{\prime}\rangle\geqslant 1-o(1).

The maximum-likelihood-estimator (MLE) problem for Problem 1.1 is an instance of the following meta-problem for k=3k=3 and f ⁣:x↦∑ijkTijkxixjxkf\colon x\mapsto\sum_{ijk}\mathbf{T}_{ijk}x_{i}x_{j}x_{k} [MR14].

For k=2k=2, this problem is just an eigenvector computation. Already for k=3k=3, it is NP-hard. Our algorithms proceed by relaxing Problem 1.2 to a convex problem. The latter can be solved either exactly or approximately (as will be the case of our faster algorithms). Under the Gaussian assumption on the noise in Problem 1.1, we show that for τ⩾ω(n3/4log⁡(n)1/4)\tau\geqslant\omega(n^{3/4}\log(n)^{1/4}) the relaxation does not substantially change the global optimum.

Montanari and Richard actually consider two variants of this model. The first we have already described. In the second, the noise is symmetrized, (to match the symmetry of potential signal tensors v⊗3v^{\otimes 3}).

It turns out that for our algorithms based on the sum-of-squares method, this kind of symmetrization is already built-in. Hence there is no difference between Problem 1.1 and Problem 1.3 for those algorithms. For our faster algorithms, such symmetrization is not built in. Nonetheless, we show that a variant of our nearly-linear-time algorithm for Problem 1.1 also solves Problem 1.3 with matching guarantees.

1 Results

We consider the degree-44 sum-of-squares relaxation for the MLE problem. (See Section 1.2 for a brief discussion about sum-of-squares. All necessary definitions are in Section 2. See [BS14] for more detailed discussion.) Note that the planted vector vv has objective value (1−o(1))τ(1-o(1))\tau for the MLE problem with high probability (assuming τ=Ω(n)\tau=\Omega(\sqrt{n}) which will always be the case for us).

There exists a polynomial-time algorithm based on the degree-44 sum-of-squares relaxation for the MLE problem that given an instance of Problem 1.1 or Problem 1.3 with τ⩾n3/4(log⁡n)1/4/ε\tau\geqslant n^{3/4}(\log n)^{1/4}/\varepsilon outputs a unit vector v′v^{\prime} with ⟨v,v′⟩⩾1−O(ε)\langle v,v^{\prime}\rangle\geqslant 1-O(\varepsilon) with probability 1−O(n−10)1-O(n^{-10}) over the randomness in the input. Furthermore, the algorithm works by rounding any solution to the relaxation with objective value at least (1−o(1))τ(1-o(1))\tau. Finally, the algorithm also certifies that all unit vectors bounded away from v′v^{\prime} have objective value significantly smaller than τ\tau for the MLE problem Problem 1.2.

We complement the above algorithmic result by the following lower bound.

We interpret a tensor-unfolding algorithm studied by Montanari and Richard as a spectral relaxation of the degree-4 sum-of-squares program for the MLE problem. This interpretation leads to an analysis that gives better guarantees in terms of signal-to-noise ratio τ\tau and also informs a more efficient implementation based on shifted matrix power iteration.

Our algorithmic results also extend in a straightforward way to tensors of order higher than 33. (See Section 7 for some details.) For simplicity we give some of these results only for the higher-order analogue of Problem 1.1; we conjecture however that all our results for Problem 1.3 generalize in similar fashion.

There is a polynomial-time algorithm, based on semidefinite programming, which on input T(x)=τ⋅⟨v0,x⟩k+A(x)\mathbf{T}(x)=\tau\cdot\langle v_{0},x\rangle^{k}+\mathbf{A}(x) returns a unit vector vv with ⟨v0,v⟩⩾1−O(ε)\langle v_{0},v\rangle\geqslant 1-O(\varepsilon) with probability 1−O(n−10)1-O(n^{-10}) over random choice of A\mathbf{A}.

There is a polynomial-time algorithm, based on semidefinite programming, which on input T(x)=τ⋅⟨v0,x⟩k+A(x)\mathbf{T}(x)=\tau\cdot\langle v_{0},x\rangle^{k}+\mathbf{A}(x) certifies that T(x)⩽τ⋅⟨v,x⟩k+O(nk/4log⁡(n)1/4)\mathbf{T}(x)\leqslant\tau\cdot\langle v,x\rangle^{k}+O(n^{k/4}\log(n)^{1/4}) for some unit vv with probability 1−O(n−10)1-O(n^{-10}) over random choice of A\mathbf{A}. This guarantees in particular that vv is close to a maximum likelihood estimator for the problem of recovering the signal v0v_{0} from the input τ⋅v0⊗k+A\tau\cdot v_{0}^{\otimes k}+\mathbf{A}.

For even kk, the above all hold, except now we recover vv with ⟨v0,v⟩2⩾1−O(ε)\langle v_{0},v\rangle^{2}\geqslant 1-O(\varepsilon), and the algorithms can be implemented in nearly linear time.

When A\mathbf{A} is a symmetric noise tensor (the higher-order analogue of Problem 1.3), (1–2) above hold. We conjecture that (3) does as well.

The last theorem, the higher-order generalization of Theorem 1.6, almost completely resolves a conjecture of Montanari and Richard regarding tensor unfolding algorithms for odd kk. We are able to prove their conjectured signal-to-noise ratio τ\tau for an algorithm that works mainly by using an unfolding of the input tensor, but our algorithm includes an extra random-rotation step to handle sparse signals. We conjecture but cannot prove that the necessity of this step is an artifact of the analysis.

2 Techniques

To maximize gg, we apply the Sum-of-Squares meta-algorithm (SoS). SoS provides a hierarchy of strong convex relaxations of Problem 1.2. Using convex duality, we can recast the optimization problem as one of efficiently certifying the upper bound on hh which shows that optima of gg are dominated by the signal. SoS efficiently finds boundedness certificates for hh of the form

Finally, we analyze a third algorithm for TPCA which simply computes the highest singular vector of a matrix unfolding of the input tensor. This algorithm was considered in depth by Montanari and Richard, who fully characterized its behavior in the case of even-order tensors (corresponding to k=4,6,8,…k=4,6,8,\ldots in Problem 1.2). They conjectured that this algorithm successfully recovers the signal vv at the signal-to-noise ratio τ\tau of Theorem 1.7 for Problem 1.1 and Problem 1.3. Up to an extra random rotations step before the tensor unfolding in the case that the input comes from Problem 1.3 (and up to logarithmic factors in τ\tau) we confirm their conjecture. We observe that their algorithm can be viewed as a method of rounding a non-optimal solution to the SoS relaxation to find the signal. We show, also, that for k=4k=4, the degree-44 SoS relaxation does no better than the simpler tensor unfolding algorithm as far as signal-to-noise ratio is concerned. However, for odd-order tensors this unfolding algorithm does not certify its own success in the way our other algorithms do.

In Theorem 1.5, we show that degree-44 SoS cannot certify that the noise polynomial A(x)=∑ijkaijkxixjxk\mathbf{A}(x)=\sum_{ijk}a_{ijk}x_{i}x_{j}x_{k} for aijka_{ijk} iid standard Gaussians satisfies A(x)⩽o(n3/4)\mathbf{A}(x)\leqslant o(n^{3/4}).

3 Related Work

There is a vast literature on tensor analogues of linear algebra problems—too vast to attempt any survey here. Tensor methods for machine learning, in particular for learning latent variable models, have garnered recent attention, e.g., with works of Anandkumar et al. [AGH+14, AGHK13]. These approaches generally involve decomposing a tensor which captures some aggregate statistics of input data into rank-one components. A recent series of papers analyzes the tensor power method, a direct analogue of the matrix power method, as a way to find rank-one components of random-case tensors [AGJ14b, AGJ14a].

Another recent line of work applies the Sum of Squares (a.k.a. Lasserre or Lasserre/Parrilo) hierarchy of convex relaxations to learning problems. See the survey of Barak and Steurer for references and discussion of these relaxations [BS14]. Barak, Kelner, and Steurer show how to use SoS to efficiently find sparse vectors planted in random linear subspaces, and the same authors give an algorithm for dictionary learning with strong provable statistical guarantees [BKS14b, BKS14a]. These algorithms, too, proceed by decomposition of an underlying random tensor; they exploit the strong (in many cases, the strongest-known) algorithmic guarantees offered by SoS for this problem in a variety of average-case settings.

Concurrently and independently of us, and also inspired by the recently-discovered applicability of tensor and sum-of-squares methods to machine learning, Barak and Moitra use SoS techniques formally related to ours to address the tensor prediction problem: given a low-rank tensor (perhaps measured with noise) only a subset of whose entries are revealed, predict the rest of the tensor entries [BM15]. They work with worst-case noise and study the number of revealed entries necessary for the SoS hierarchy to successfully predict the tensor. By constrast, in our setting, the entire tensor is revealed, and we study the signal-to-noise threshold necessary for SoS to recover its principal component under distributional assumptions on the noise that allow us to avoid worst-case hardness behavior.

Since Barak and Moitra work in a setting where few tensor entries are revealed, they are able to use algorithmic techniques and lower bounds from the study of sparse random constraint satisfaction problems (CSPs), in particular random 3XOR [GK01, FGK05, FO07, FKO06]. The tensors we study are much denser. In spite of the density (and even though our setting is real-valued), our algorithmic techniques are related to the same spectral refutations of random CSPs. However our lower bound techniques do not seem to be related to the proof-complexity techniques that go into sum-of-squares lower bound results for random CSPs.

The analysis of tractable tensor decomposition in the rank one plus noise model that we consider here (the spiked tensor model) was initiated by Montanari and Richard, whose work inspired the current paper [MR14]. They analyze a number of natural algorithms and find that tensor unfolding algorithms, which use the spectrum of a matrix unfolding of the input tensor, are most robust to noise. Here we consider more powerful convex relaxations, and in the process we tighten Montanari and Richard’s analysis of tensor unfolding in the case of odd-order tensors. In concurrent and independent work, Zheng and Tomioka also give a tight analysis of tensor unfolding for the asymmetric version of the spiked model of tensor pca (Problem 1.1) [ZT15, Theorem 1].

Related to our lower bound, Montanari, Reichman, and Zeitouni (MRZ) prove strong impossibility results for the problem of detecting rank-one perturbations of Gaussian matrices and tensors using any eigenvalue of the matrix or unfolded tensor; they are able to characterize the precise threshold below which the entire spectrum of a perturbed noise matrix or unfolded tensor becomes indistinguishable from pure noise [MRZ14]. This lower bound is incomparable to our lower bound for the degree-4 SoS relaxation. The MRZ lower bound considers fine-grained information about the spectrum of a single matrix associated with the detection problem. Our lower bound considers coarser information (just the top eigenvalue) but it applies to a wide range of matrices associated with the problem (all matrices generated via the degree-4 sum-of-squares proof system).

Preliminaries

We employ the usual Loewner (a.k.a. positive semi-definite) ordering ⪰\succeq on Hermitian matrices.

We will be heavily concerned with tensors and matrix flattenings thereof. In general, boldface capital letters T\mathbf{T} denote tensors and ordinary capital letters denote matrices AA. We adopt the convention that unless otherwise noted for a tensor T\mathbf{T} the matrix TT is the squarest-possible unfolding of T\mathbf{T}. If T\mathbf{T} has even order kk then TT has dimensions nk/2×nk/2n^{k/2}\times n^{k/2}. For odd kk it has dimensions n⌊k/2⌋×n⌈k/2⌉n^{\lfloor k/2\rfloor}\times n^{\lceil k/2\rceil}. All tensors, matrices, vectors, and scalars in this paper are real.

For a kk-tensor T\mathbf{T}, we write T(v){\mathbf{T}}(v) for ⟨v⊗k,T⟩\langle v^{\otimes k},\mathbf{T}\rangle. Thus, T(x)\mathbf{T}(x) is a homogeneous real polynomial of degree kk.

We use Sk\mathcal{S}_{k} to denote the symmetric group on kk elements. For a kk-tensor T\mathbf{T} and π∈Sk\pi\in\mathcal{S}_{k}, we denote by Tπ{\mathbf{T}}^{\pi} the kk-tensor with indices permuted according to π\pi, so that Tαπ=Tπ−1(α){\mathbf{T}}^{\pi}_{\alpha}={\mathbf{T}}_{\pi^{-1}(\alpha)}. A tensor T\mathbf{T} is symmetric if for all π∈Sk\pi\in\mathcal{S}_{k} it is the case that Tπ=T{\mathbf{T}}^{\pi}=\mathbf{T}. (Such tensors are sometimes called “supersymmetric.”)

For clarity, most of our presentation focuses on 33-tensors. For an n×nn\times n 33-tensor T\mathbf{T}, we use TiT_{i} to denote its n×nn\times n matrix slices along the first mode, i.e., (Ti)j,k=Ti,j,k(T_{i})_{j,k}=\mathbf{T}_{i,j,k}.

2 Polynomials and Matrices

3 The Sum of Squares (SoS) Algorithm

Pseudo-distributions were first introduced in [BBH+12] and are surveyed in [BS14].

We employ the standard result that, up to negligible issues of numerical accuracy, if there exists a degree-dd pseudo-distribution satisfying constraints {p0(x)=0,…,pm(x)=0}\{p_{0}(x)=0,\ldots,p_{m}(x)=0\}, then it can be found in time nO(d)n^{O(d)} by solving a semidefinite program of size nO(d)n^{O(d)}. (See [BS14] for references.)

Certifying Bounds on Random Polynomials

To better exploit the benefits of square matrices, we bound the maxima of degree-33 homogeneous ff by a degree-44 polynomial. In the case that ff is multi-linear, we have the polynomial identity f(x)=13⟨x,∇f(x)⟩f(x)=\frac{1}{3}\langle x,\nabla f(x)\rangle. Using Cauchy-Schwarz, we then get f(x)⩽13∥x∥∥∇f(x)∥f(x)\leqslant\frac{1}{3}\|x\|\|\nabla f(x)\|. This inequality suggests using the degree-44 polynomial ∥∇f(x)∥2\|\nabla f(x)\|^{2} as a bound on ff. Note that local optima of ff on the sphere occur where ∇f(v)∝v\nabla f(v)\propto v, and so this bound is tight at local maxima. Given a random homogeneous ff, we will associate a degree-44 polynomial related to ∥∇f∥2\|\nabla f\|^{2} and show that this polynomial yields the best possible degree-44 SoS-certifiable bound on max⁡∥v∥=1f(v)\max_{\|v\|=1}f(v).

We observe that for ff multi-linear in the coordinates xix_{i} of xx, up to a constant factor we may take the matrices AiA_{i} to be matrix representations of ∂if\partial_{i}f, so that ∑iAi⊗Ai\sum_{i}A_{i}\otimes A_{i} is a matrix representation of the polynomial ∥∇f∥2\|\nabla f\|^{2}. This choice of AiA_{i} may not, however, yield the optimal spectral bound λ2\lambda^{2}.

The following theorem is the reason for our definition of λ\lambda-boundedness.

We now state the degree-33 case of a general λ\lambda-boundedness fact for homogeneous polynomials with random coefficients. The SoS-certifiable bound for a random degree-33 polynomial this provides is the backbone of our SoS algorithm for tensor PCA in the spiked tensor model.

Let A\mathbf{A} be a 33-tensor with independent entries from N(0,1)\mathcal{N}(0,1). Then A(x)\mathbf{A}(x) is λ\lambda-bounded with λ=O(n3/4log⁡(n)1/4)\lambda=O(n^{3/4}\log(n)^{1/4}), with high probability.

The full statement and proof of this theorem, generalized to arbitrary-degree homogeneous polynomials, may be found as Theorem B.5; we prove the statement above as a corollary in Section B. Here provide a proof sketch.

Immediate by combining Theorem 3.3 with Theorem 3.2. ∎

Polynomial-Time Recovery via Sum of Squares

The following theorem characterizes the success of Algorithm 4.1 and Algorithm 4.2

if there exists a sufficiently good upper bound on A(x)\mathbf{A}(x) (or in the case of the symmetric noise input, on Aπ(x)\mathbf{A}^{\pi}(x) for every π∈S3\pi\in\mathcal{S}_{3}) which is degree-4 SoS certifiable, then the vector recovered by the algorithm will be very close to vv, and that

in the case of A\mathbf{A} with independent entries from N(0,1)\mathcal{N}(0,1), such a bound exists with high probability.

Conveniently, Item 2 is precisely the content of Corollary 3.4. The following lemma expresses Item 1.

Rewriting ⟨v0⊗3,x⊗3⟩\langle v_{0}^{\otimes 3},x^{\otimes 3}\rangle as ⟨v0,x⟩3\langle v_{0},x\rangle^{3}, we obtain

We discuss here a modified TPCA model, which will illustrate the qualitative differences between the new tensor PCA algorithms we propose in this paper and previously-known algorithms. The model is semi-random and semi-adversarial. Such models are often used in average-case complexity theory to distinguish between algorithms which work by solving robust maximum-likelihood-style problems and those which work by exploiting some more fragile property of a particular choice of input distribution.

Here we show that Algorithm 4.1 succeeds in recovering vv in the semi-random model.

Let T′\mathbf{T}^{\prime} be the semi-random-model tensor PCA input, with τ⩾n3/4log⁡(n)1/4/ε\tau\geqslant n^{3/4}\log(n)^{1/4}/\varepsilon. With high probability over randomness in T′\mathbf{T}^{\prime}, Algorithm 4.1 outputs a vector vv with ⟨v,v0⟩⩾1−O(ε)\langle v,v_{0}\rangle\geqslant 1-O(\varepsilon).

Linear Time Recovery via Further Relaxation

We now attack the problem of speeding up the algorithm from the preceding section. We would like to avoid solving a large semidefinite program to optimality: our goal is to instead use much faster linear-algebraic computations—in particular, we will recover the tensor PCA signal vector by performing a single singular vector computation on a relatively small matrix. This will complete the proofs of Theorem 1.7 and Theorem 1.6, yielding the desired running time.

Our SoS algorithm in the preceding section turned on the existence of the λ\lambda-boundedness certificate ∑iAi⊗Ai\sum_{i}A_{i}\otimes A_{i}, where AiA_{i} are the slices of a random tensor A\mathbf{A}. Let T=τ⋅v0⊗3+A\mathbf{T}=\tau\cdot v_{0}^{\otimes 3}+\mathbf{A} be the spiked-tensor input to tensor PCA. We could look at the matrix ∑iTi⊗Ti\sum_{i}T_{i}\otimes T_{i} as a candidate λ\lambda-boundedness certificate for T(x)\mathbf{T}(x). The spectrum of this matrix must not admit the spectral bound that ∑iAi⊗Ai\sum_{i}A_{i}\otimes A_{i} does, because T(x)\mathbf{T}(x) is not globally bounded: it has a large global maximum near the signal vv. This maximum plants a single large singular value in the spectrum of ∑iTi⊗Ti\sum_{i}T_{i}\otimes T_{i}. The associated singular vector is readily decoded to recover the signal.

Before stating and analyzing this fast linear-algebraic algorithm, we situate it more firmly in the SoS framework. In the following, we discuss spectral SoS, a convex relaxation of Problem 1.2 obtained by weakening the full-power SoS relaxation. We show that the spectrum of the aforementioned ∑iTi⊗Ti\sum_{i}T_{i}\otimes T_{i} can be viewed as approximately solving the spectral SoS relaxation. This gives the fast, certifying algorithm of Theorem 1.7. We also interpret the tensor unfolding algorithm given by Montanari and Richard for TPCA in the spiked tensor model as giving a more subtle approximate solution to the spectral SoS relaxation. We prove a conjecture by those authors that the algorithm successfully recovers the TPCA signal at the same signal-to-noise ratio as our other algorithms, up to a small pre-processing step in the algorithm; this proves Theorem 1.6 [MR14]. This last algorithm, however, succeeds for somewhat different reasons than the others, and we will show that it consequently fails to certify its own success and that it is not robust to a certain kind of semi-adversarial choice of noise.

To obtain spectral SoS, the convex relaxation of Problem 1.2 which we will be able to (approximately) solve quickly in the random case, we first need to return to the full-strength SoS relaxation and examine it from a more linear-algebraic standpoint.

A polynomial may have many matrix representations, but a pseudo-distribution has just one: a matrix representation of a pseudo-distribution must obey strong symmetry conditions in order to assign the same pseudo-expectation to every representation of the same polynomial. We will have much more to say about constructing matrices satisfying these symmetry conditions when we state and prove our lower bounds, but here we will in fact profit from relaxing these symmetry constraints.

It may not be immediately obvious why this program optimizes only over MM which are matrix representations of pseudo-distributions. If, however, some MM does not obey the requisite symmetries, then min⁡Mp∈Mp⟨M,Mp⟩=−∞\min_{M_{p}\in\mathcal{M}_{p}}\langle M,M_{p}\rangle=-\infty, since the asymmetry may be exploited by careful choice of Mp∈MpM_{p}\in\mathcal{M}_{p}. Thus, at optimality this program yields MM which is the matrix representation of a pseudo-distribution {x}\{x\} satisfying {∥x∥2−1=0}\{\|x\|^{2}-1=0\}.

1.2 Relaxing to the Degree-444 Dual

By weak duality, we can interchange the min⁡\min and the max⁡\max in (5.2) to obtain the dual program:

We call this dual program the spectral SoS relaxation of max⁡∥x∥=1p(x)\max_{\|x\|=1}p(x). If p=∑i⟨x,Aix⟩p=\sum_{i}\langle x,A_{i}x\rangle for A\mathbf{A} with independent entries from N(0,1)\mathcal{N}(0,1), the spectral SoS relaxation achieves the same bound as our analysis of the full-strength SoS relaxation: for such pp, the spectral SoS relaxation is at most O(n3/2log⁡(n)1/2)O(n^{3/2}\log(n)^{1/2}) with high probability. The reason is exactly the same as in our analysis of the full-strength SoS relaxation: the matrix ∑iAi⊗Ai\sum_{i}A_{i}\otimes A_{i}, whose spectrum we used before to bound the full-strength SoS relaxation, is still a feasible dual solution.

Let T=τ⋅v0⊗3+A\mathbf{T}=\tau\cdot v_{0}^{\otimes 3}+\mathbf{A} be the spiked-tensor input to tensor PCA. We know from our initial characterization of SoS proofs of boundedness for degree-33 polynomials that the polynomial T′(x):=(x⊗x)T(∑iTi⊗Ti)(x⊗x)\mathbf{T}^{\prime}(x)\mathrel{\mathop{:}}=(x\otimes x)^{T}(\sum_{i}T_{i}\otimes T_{i})(x\otimes x) gives SoS-certifiable upper bounds on T(x)\mathbf{T}(x) on the unit sphere. We consider the spectral SoS relaxation of max⁡∥x∥=1T′(x)\max_{\|x\|=1}\mathbf{T}^{\prime}(x),

The following theorem describes the behavior of Algorithm 5.1 and Algorithm 5.2 and gives a proof of Theorem 1.7 and Corollary 1.7.

With high probability, Algorithm 5.1 returns vv with ⟨v,v0⟩2⩾1−O(ε)\langle v,v_{0}\rangle^{2}\geqslant 1-O(\varepsilon).

If Algorithm 5.2 outputs certify then T(x)⩽τ⋅⟨v,x⟩3+O(n3/4log⁡(n)1/4)\mathbf{T}(x)\leqslant\tau\cdot\langle v,x\rangle^{3}+O(n^{3/4}\log(n)^{1/4}) (regardless of the distribution of A\mathbf{A}). If A\mathbf{A} is distributed as above, then Algorithm 5.2 outputs certify with high probability.

Both Algorithm 5.1 and Algorithm 5.2 can be implemented in time O(n4log⁡(1/ε))O(n^{4}\log(1/\varepsilon)).

A similar fact to Lemma 5.6 appears in [MR14].

The proofs of Lemma 5.4 and Lemma 5.6 follow here. The proof of Lemma 5.5 uses only standard concentration of measure arguments; we defer it to Section B.

Let u,wu,w be the top left and right singular vectors of MM. We have

Let v0,v′,V′,vv_{0},v^{\prime},V^{\prime},v, be as in the lemma statement. We know vv is the maximizer of max⁡∥w∥,∥w′∥=1wTV′w′\max_{\|w\|,\|w^{\prime}\|=1}w^{T}V^{\prime}w^{\prime}. By assumption,

Thus, the top singular value of V′V^{\prime} is at least 1−O(ε)1-O(\varepsilon), and since ∥v′∥\|v^{\prime}\| is a unit vector, the Frobenius norm of V′V^{\prime} is 11 and so all the rest of the singular values are O(ε)O(\varepsilon). Expressing v0v_{0} in the right singular basis of V′V^{\prime} and examining the norm of V′v0V^{\prime}v_{0} completes the proof. ∎

In the first case, we start by observing that it is enough to find a vector ww which has ⟨w,v′⟩⩾1−ε\langle w,v^{\prime}\rangle\geqslant 1-\varepsilon, where v′v^{\prime} is a top singular vector of MM. Let λ1,λ2\lambda_{1},\lambda_{2} be the top two singular values of MM. The analysis of the algorithm already showed that λ1/λ2⩾Ω(1/ε)\lambda_{1}/\lambda_{2}\geqslant\Omega(1/\varepsilon). Standard analysis of the matrix power method now yields that O(log⁡(1/ε))O(\log(1/\varepsilon)) iterations will suffice.

3 Nearly-Linear-Time Recovery via Tensor Unfolding and Spectral SoS

Despite its a priori simplicity, the analysis of Algorithm 5.7 is more subtle than for any of our other algorithms. This would not be true for even-order tensors, for which the square matrix unfolding tensor has one singular value asymptotically larger than all the rest, and indeed the corresponding singular vector is well-correlated with v0v_{0}. However, in the case of odd-order tensors the unfolding has no spectral gap. Instead, the signal v0v_{0} has some second-order effect on the spectrum of the matrix unfolding, which is enough to recover it.

Again by triangle inequality, uTMu⩾v0TMv=τ2−O(ετ2)u^{T}Mu\geqslant v_{0}^{T}Mv=\tau^{2}-O(\varepsilon\tau^{2}). So rearranging we get ⟨u,v0⟩2⩾1−O(ε)\langle u,v_{0}\rangle^{2}\geqslant 1-O(\varepsilon) as desired. ∎

The following lemma is a consequence of standard matrix concentration inequalities; we defer its proof to Section B, Lemma B.10.

Immediate from Lemma 5.9, Lemma 5.10, and Lemma 5.11. ∎

4 Fast Recovery in the Semi-Random Model

There is a qualitative difference between the aggregate matrix statistics needed by our certifying algorithms (Algorithm 4.1, Algorithm 4.2, Algorithm 5.1, Algorithm 5.2) and those needed by rounding the tensor unfolding solution spectral SoS Algorithm 5.7. In a precise sense, the needs of the latter are greater. The former algorithms rely only on first-order behavior of the spectra of a tensor unfolding, while the latter relies on second-order spectral behavior. Since it uses second-order properties of the randomness, Algorithm 5.7 fails in the semi-random model.

The argument that that Algorithm 5.1 and Algorithm 5.2 still succeed in the semi-random model is routine; for completeness we discuss here the necessary changes to the proof of Theorem 5.3. The non-probabilistic certification claims made in Theorem 5.3 are independent of the input model, so we show that Algorithm 5.1 still finds the signal with high probability and that Algorithm 5.2 still fails only with only a small probability.

In the semi-random model, ε⩾n−1/4\varepsilon\geqslant n^{-1/4} and τ⩾n3/4log⁡(n)1/4/ε\tau\geqslant n^{3/4}\log(n)^{1/4}/\varepsilon, with high probability, Algorithm 5.1 returns vv with ⟨v,v0⟩2⩾1−O(ε)\langle v,v_{0}\rangle^{2}\geqslant 1-O(\varepsilon) and Algorithm 5.2 outputs certify.

5 Fast Recovery with Symmetric Noise

We suppose now that A\mathbf{A} is a symmetric Gaussian noise tensor; that is, that A\mathbf{A} is the average of A0π\mathbf{A}_{0}^{\pi} over all π∈S3\pi\in\mathcal{S}_{3}, for some order-33 tensor A0\mathbf{A}_{0} with iid standard Gaussian entries.

Our previous techniques fail in this symmetric noise scenario due to lack of independence between the entries of the noise tensor. However, we sidestep that issue here by restricting our attention to an asymmetric block of the input tensor.

The resulting algorithm is not precisely identical to the tensor unfolding algorithm investigated by Montanari and Richard, but is based on tensor unfolding with only superficial modifications.

It is possible to implement each iteration of the matrix power method in Algorithm 5.14 in linear time. We focus on multiplying a vector by MXM_{X} in linear time; the other cases follow similarly.

To accomplish this, we simply reflatten our tensors. Let VV be the nn-by-nn matrix flattening of vv. Then we compute the matrix RTPYR⋅V⋅RTPZTRR^{T}P_{Y}R\cdot V\cdot R^{T}{P_{Z}}^{T}R, and return its flattening back into an n2n^{2}-dimensional vector, and this will be equal to (R⊗R)T(PY⊗PZ)(R⊗R) v(R\otimes R)^{T}(P_{Y}\otimes P_{Z})(R\otimes R)\,v. This equivalence follows by taking the singular value decomposition V=∑iλiuiwiTV=\sum_{i}\lambda_{i}u_{i}w_{i}^{T}, and noting that v=∑iλiui⊗wiv=\sum_{i}\lambda_{i}u_{i}\otimes w_{i}.

So γ2∥PRu∥2\gamma^{2}\|PRu\|^{2} is the sum of the squares of mm independent variables drawn from N(0,1/n)\mathcal{N}(0,1/n). By a Bernstein inequality, \big{|}\gamma^{2}\|PRu\|^{2}-m/n\big{|}\leqslant O(\sqrt{m/n^{2}}\log m) with high probability. Also by a Bernstein inequality, γ2−1<O(1/nlog⁡n)\gamma^{2}-1<O(\sqrt{1/n}\log n) with high probability. ∎

For τ⩾n3/4/ε\tau\geqslant n^{3/4}/\varepsilon, with high probability, Algorithm 5.14 recovers a vector vv with ⟨v,v0⟩⩾1−O(ε)\langle v,v_{0}\rangle\geqslant 1-O(\varepsilon) when A\mathbf{A} is a symmetric Gaussian noise tensor (as in Problem 1.3) and ε⩾log⁡(n)/n\varepsilon\geqslant\log(n)/\sqrt{n}.

Name the projections UX:=(PY⊗PZ)UPXU_{X}\mathrel{\mathop{:}}=(P_{Y}\otimes P_{Z})UP_{X}, UY:=(PZ⊗PX)UPYU_{Y}\mathrel{\mathop{:}}=(P_{Z}\otimes P_{X})UP_{Y}, and UZ:=(PX⊗PY)UPZU_{Z}\mathrel{\mathop{:}}=(P_{X}\otimes P_{Y})UP_{Z}.

First off, U=τ(Rv0)⊗3+A′{\mathbf{U}}=\tau(Rv_{0})^{\otimes 3}+\mathbf{A}^{\prime} where A′\mathbf{A}^{\prime} is a symmetric Gaussian tensor (distributed identically to A\mathbf{A}). This follows by noting that multiplication by R⊗3R^{\otimes 3} commutes with permutation of indices, so that (R⊗3B)π=R⊗3Bπ(R^{\otimes 3}{\mathbf{B}})^{\pi}=R^{\otimes 3}{\mathbf{B}}^{\pi}, where we let B{\mathbf{B}} be the asymmetric Gaussian tensor so that A=∑π∈S3Bπ\mathbf{A}=\sum_{\pi\in\mathcal{S}_{3}}{\mathbf{B}}^{\pi}. Then A′=R⊗3∑π∈S3Bπ=∑π∈S3(R⊗3B)π\mathbf{A}^{\prime}=R^{\otimes 3}\sum_{\pi\in\mathcal{S}_{3}}{\mathbf{B}}^{\pi}=\sum_{\pi\in\mathcal{S}_{3}}(R^{\otimes 3}{\mathbf{B}})^{\pi}. This is identically distributed with A\mathbf{A}, as follows from the rotational symmetry of B{\mathbf{B}}.

Thus UX=τ(PY⊗PZ)(R⊗R)(v0⊗v0)(PXRv0)T+(PY⊗PZ)A′PXU_{X}=\tau(P_{Y}\otimes P_{Z})(R\otimes R)(v_{0}\otimes v_{0})(P_{X}Rv_{0})^{T}+(P_{Y}\otimes P_{Z})A^{\prime}P_{X}, and

Let SS refer to Expression 5.5. By Lemma 5.16, \big{|}\|PRv_{0}\|^{2}-\tfrac{1}{3}\big{|}<O(\sqrt{1/n}\log n) with high probability for P∈{PX,PY,PZ}P\in\{P_{X},P_{Y},P_{Z}\}. Hence S=(19±O(1/nlog⁡n))τ2(PXRv0)(PXRv0)TS=(\tfrac{1}{9}\pm O(\sqrt{1/n}\log n))\tau^{2}(P_{X}Rv_{0})(P_{X}Rv_{0})^{T} and ∥S∥=(127±O(1/nlog⁡n))τ2\|S\|=(\tfrac{1}{27}\pm O(\sqrt{1/n}\log n))\tau^{2}.

Let CC refer to Expression 5.6 so that Expression 5.7 is CTC^{T}. Let also A′′=(PY⊗PZ)A′PXA^{\prime\prime}=(P_{Y}\otimes P_{Z})A^{\prime}P_{X}. Note that, once the identically-zero rows and columns of A′′A^{\prime\prime} are removed, A′′A^{\prime\prime} is a matrix of iid standard Gaussian entries. Finally, let v′′=PYRv0⊗PZRv0v^{\prime\prime}=P_{Y}Rv_{0}\otimes P_{Z}Rv_{0}. By some substitution and by noting that ∥PXR∥⩽1\|P_{X}R\|\leqslant 1, we have that ∥C∥⩽τ ∥v0v′′TA′′∥\|C\|\leqslant\tau\,\|v_{0}{v^{\prime\prime}}^{T}A^{\prime\prime}\|. Hence by Lemma B.10, ∥C∥⩽O(ετ2)\|C\|\leqslant O(\varepsilon\tau^{2}).

The recovered eigenvector vXv_{X} satisfies ⟨vX,MXvX⟩⩾Ω(τ2)\langle v_{X},M_{X}v_{X}\rangle\geqslant\Omega(\tau^{2}) and ⟨vX,(MX−S)vX⟩⩽O(ετ2)\langle v_{X},(M_{X}-S)v_{X}\rangle\leqslant O(\varepsilon\tau^{2}) and therefore ⟨vX,SvX⟩=(127±O(ε+1/nlog⁡n))τ2\langle v_{X},Sv_{X}\rangle=(\tfrac{1}{27}\pm O(\varepsilon+\sqrt{1/n}\log n))\tau^{2}. Substituting in the expression for SS, we conclude that ⟨PXRv0,vX⟩=(13±O(ε+1/nlog⁡n))\langle P_{X}Rv_{0},v_{X}\rangle=(\tfrac{1}{\sqrt{3}}\pm O(\varepsilon+\sqrt{1/n}\log n)).

The analyses for vYv_{Y} and vZv_{Z} follow in the same way. Hence

At the same time, since vXv_{X}, vYv_{Y}, and vZv_{Z} are each orthogonal to each other, ∥vX+vY+vZ∥=3\|v_{X}+v_{Y}+v_{Z}\|=\sqrt{3}. Hence with the output vector being v:=R−1(vX+vY+vZ)/∥vX+vY+vZ∥v\mathrel{\mathop{:}}=R^{-1}(v_{X}+v_{Y}+v_{Z})/\|v_{X}+v_{Y}+v_{Z}\|, we have

6 Numerical Simulations

We report now the results of some basic numerical simulations of the algorithms from this section. In particular, we show that the asymptotic running time differences among Algorithm 5.1, Algorithm 5.7 implemented naïvely, and the linear-time implementation of Algorithm 5.7 are apparent at reasonable values of nn, e.g. n=200n=200.

Specifics of our experiments are given in Figure 1. We find pronounced differences between all three algorithms. The naïve implementation of Algorithm 5.7 is markedly slower than the linear implementation, as measured either by number of matrix-vector multiplies or processor time. Algorithm 5.1 suffers greatly from the need to construct an n2×n2n^{2}\times n^{2} matrix; although we do not count the time to construct this matrix against its reported running time, the memory requirements are so punishing that we were unable to collect data beyond n=100n=100 for this algorithm.

Lower Bounds

We will now prove lower bounds on the performance of degree-44 SoS on random instances of the degree-44 and degree-33 homogeneous polynomial maximization problems. As an application, we show that our analysis of degree-44 for Tensor PCA is tight up to a small logarithmic factor in the signal-to-noise ratio.

The existence of the maps η\eta depending only on the random part A\mathbf{A} of the tensor PCA input v0⊗3+Av_{0}^{\otimes 3}+\mathbf{A} formalizes the claim from Theorem 1.5 that no algorithm can reliably recover v0v_{0} from the pseudo-distribution η(A)\eta(\mathbf{A}).

Additionally, the lower-bound construction holds for the symmetric noise model also: the input tensor A\mathbf{A} is symmetrized wherever it occurs in the construction, so it does not matter if it had already been symmetrized beforehand.

The rest of this section is devoted to proving these theorems, which we eventually accomplish in Section 6.2.

The general outline of the proof will be as follows:

But before we can state a formal version of our theorem, we will need a few facts about polynomials, pseudo-distributions, matrices, vectors, and how they are related by symmetries under actions of permutation groups.

1 Polynomials, Vectors, Matrices, and Symmetries, Redux

Here we further develop the matrix view of SoS presented in Section 5.1.1.

in order that they assign consistent values to each representation of the same polynomial. We call such matrices maximally symmetric (following Doherty and Wehner [DW12]).

The degree dd will always be clear from context.

1.2 The Monomial-Indexed (i.e. Symmetric) Subspace

We let Π\Pi be the projector to this subspace. For any maximally-symmetric MM we have ΠMΠ=M\Pi M\Pi=M, but the reverse implication is not true (for readers familiar with quantum information: any MM which has M=ΠMΠM=\Pi M\Pi is Bose-symmetric, but may not be PPT-symmetric; maximally symmetric matrices are both. See [DW12] for further discussion.)

1.3 Maximally-Symmetric Matrices from Tensors

2 Formal Statement of the Lower Bound

We will warm up with the degree-44 lower bound, which is conceptually somewhat simpler.

Let A\mathbf{A} be a 44-tensor and let λ>0\lambda>0 be a function of nn. Suppose the following conditions hold:

A\mathbf{A} is significantly correlated with ∑π∈S4Aπ\sum_{\pi\in\mathcal{S}_{4}}\mathbf{A}^{\pi}. ⟨A,∑π∈S4Aπ⟩⩾Ω(n4)\langle\mathbf{A},\sum_{\pi\in\mathcal{S}_{4}}\mathbf{A}^{\pi}\rangle\geqslant\Omega(n^{4}).

Permutations have lower-bounded spectrum. For every π∈S4\pi\in\mathcal{S}_{4}, the Hermitian n2×n2n^{2}\times n^{2} unfolding 12(Aπ+(Aπ)T)\frac{1}{2}(A^{\pi}+(A^{\pi})^{T}) of Aπ\mathbf{A}^{\pi} has no eigenvalues smaller than −λ2-\lambda^{2}.

Then n3/2δ2′+n2δ2⩽O(1)n^{3/2}\delta_{2}^{\prime}+n^{2}\delta_{2}\leqslant O(1).

The degree-33 version of our lower bound requires bounds on the spectra of the flattenings not just of the 33-tensor A\mathbf{A} itself but also of the flattenings of an associated 44-tensor, which represents the polynomial ⟨x⊗2,(∑iAi⊗Ai)x⊗2⟩\langle x^{\otimes 2},(\sum_{i}A_{i}\otimes A_{i})x^{\otimes 2}\rangle.

Let A\mathbf{A} be a 33-tensor and let λ>0\lambda>0 be a function of nn. Suppose the following conditions hold:

A\mathbf{A} is significantly correlated with ∑π∈S3Aπ\sum_{\pi\in\mathcal{S}_{3}}\mathbf{A}^{\pi}. ⟨A,∑π∈S3Aπ⟩⩾Ω(n3)\langle\mathbf{A},\sum_{\pi\in\mathcal{S}_{3}}\mathbf{A}^{\pi}\rangle\geqslant\Omega(n^{3}).

Permutations have lower-bounded spectrum. For every π∈S3\pi\in\mathcal{S}_{3}, we have

Then nδ1+n3/2δ2′+n2δ2⩽O(1)n\delta_{1}+n^{3/2}\delta_{2}^{\prime}+n^{2}\delta_{2}\leqslant O(1).

Then there is a degree-44 pseudo-distribution {x}\{x\} satisfying {∥x∥22=1}\{\|x\|_{2}^{2}=1\} so that

We prove the degree-33 corollary; the degree-44 case is almost identical using Theorem 6.3 and Lemma B.12 in place of their degree-33 counterparts.

Let A\mathbf{A} be a 33-tensor. If A\mathbf{A} satisfies the conditions of Theorem 6.4 with λ=O(n3/4log⁡(n)1/4)\lambda=O(n^{3/4}\log(n)^{1/4}), we let η(A)\eta(\mathbf{A}) be the pseudo-distribution described there, with

3 In-depth Preliminaries for Pseudo-Expectation Symmetries

Let D8<S4\mathcal{D}_{8}<\mathcal{S}_{4} be given by D8=⟨(12),(34),(13)(24)⟩\mathcal{D}_{8}=\langle(12),(34),(13)(24)\rangle. Let C3={(),σ,σ2}=⟨σ⟩\mathcal{C}_{3}=\{(),\sigma,\sigma^{2}\}=\langle\sigma\rangle, where ()() denotes the identity in S4\mathcal{S}_{4}. Then {gh:g∈D8,h∈C3}=S4\{gh:g\in\mathcal{D}_{8},h\in\mathcal{C}_{3}\}=\mathcal{S}_{4}.

The proof is routine; we provide it here for completeness. Note that C3\mathcal{C}_{3} is a subgroup of order 33 in the alternating group A4\mathcal{A}_{4}. This alternating group can be decomposed as A4=K4⋅C3\mathcal{A}_{4}=\mathcal{K}_{4}\cdot\mathcal{C}_{3}, where K4=⟨(12)(34),(13)(24)⟩\mathcal{K}_{4}=\langle(12)(34),(13)(24)\rangle is a normal subgroup of A4\mathcal{A}_{4}. We can also decompose S4=C2⋅A4\mathcal{S}_{4}=\mathcal{C}_{2}\cdot\mathcal{A}_{4} where C2=⟨(12)⟩\mathcal{C}_{2}=\langle(12)\rangle and A4\mathcal{A}_{4} is a normal subgroup of S4\mathcal{S}_{4}. Finally, D8=C2⋅K4\mathcal{D}_{8}=\mathcal{C}_{2}\cdot\mathcal{K}_{4} so by associativity, S4=C2⋅A4=C2⋅K4⋅C3=D8⋅C3\mathcal{S}_{4}=\mathcal{C}_{2}\cdot\mathcal{A}_{4}=\mathcal{C}_{2}\cdot\mathcal{K}_{4}\cdot\mathcal{C}_{3}=\mathcal{D}_{8}\cdot\mathcal{C}_{3}. ∎

For any subset S⊆S4S\subseteq\mathcal{S}_{4}, we have {ghs:g∈D8,h∈C3,s∈S}=S4\{ghs:g\in\mathcal{D}_{8},h\in\mathcal{C}_{3},s\in S\}=\mathcal{S}_{4}.

and so M′∈Sym⁡MM^{\prime}\in\operatorname{Sym}M. ∎

We make an useful observation about the nontrivial permutations of MM, in the special case that M=AATM=AA^{T} for some 33-tensor A\mathbf{A}.

We observe that AAT[(j1,j2),(j3,j4)]=∑iAij1j2Aij3j4AA^{T}[(j_{1},j_{2}),(j_{3},j_{4})]=\sum_{i}A_{ij_{1}j_{2}}A_{ij_{3}j_{4}} and that (∑iAi⊗Ai)[(j1,j2),(j3,j4)]=∑iAij1j3Aij2j4(\sum_{i}A_{i}\otimes A_{i})[(j_{1},j_{2}),(j_{3},j_{4})]=\sum_{i}A_{ij_{1}j_{3}}A_{ij_{2}j_{4}}. Multiplication by PP on the right has the effect of switching the order of the second indexing pair, so [(∑iAi⊗Ai)P][(j1,j2),(j3,j4)]=∑iAij1j4Aij2j3[(\sum_{i}A_{i}\otimes A_{i})P][(j_{1},j_{2}),(j_{3},j_{4})]=\sum_{i}A_{ij_{1}j_{4}}A_{ij_{2}j_{3}}. From this it is easy to see that σ⋅AAT=(234)⋅AAT=(∑iAi⊗Ai)P\sigma\cdot AA^{T}=(234)\cdot AA^{T}=(\sum_{i}A_{i}\otimes A_{i})P.

from which we see that σ2⋅AAT=∑iAi⊗AiT\sigma^{2}\cdot AA^{T}=\sum_{i}A_{i}\otimes A_{i}^{T}. ∎

4 Construction of Initial Pseudo-Distributions

We begin by discussing how to create an initial guess at a pseudo-distribution whose third moments are highly correlated with the polynomial A(x)\mathbf{A}(x). This initial guess will be a valid pseudo-distribution, but will fail to be on the unit sphere, and so will require some repairing later on. For now, the method of creating this initial pseudo-distribution involves using a combination of symmetrization techniques to ensure that the matrices we construct are well defined as linear functionals over polynomials, and spectral techniques to establish positive-semidefiniteness of these matrices.

where B⪰0B\succeq 0 and is full rank. Then M⪰0M\succeq 0 if and only if D⪰CB−1CTD\succeq CB^{-1}C^{T}.

We would ideally take DD to be the spectrally-least maximally-symmetric matrix so that D⪰CB−1CTD\succeq CB^{-1}C^{T}. But this object might not be well defined, so we instead take the following substitute.

4.2 Symmetries at Degree Three

Each matrix in the sum defining MM is positive-semidefinite, so M⪰0M\succeq 0. Each DiD_{i} is maximally symmetric and therefore so is ∑i=1kDi\sum_{i=1}^{k}D_{i}. We know that ML⁡∣3=∑i=1kML⁡∣3iM_{\operatorname{\mathcal{L}}|_{3}}=\sum_{i=1}^{k}M_{\operatorname{\mathcal{L}}|_{3}}^{i} is maximally-symmetric, so it follows that MM is the matrix representation of a valid pseudo-expectation. ∎

5 Getting to the Unit Sphere

In particular, since (∥x∥22−1)(\|x\|_{2}^{2}-1) is in the kernel of L⁡′\operatorname{\mathcal{L}}^{\prime}, either λmin⁡L⁡′=0\operatorname{\lambda_{min}}\operatorname{\mathcal{L}}^{\prime}=0 or

The condition p⊥(∥x∥22−1)p\perp(\|x\|_{2}^{2}-1) yields p0=∑ipiip_{0}=\sum_{i}p_{ii}. Substituting into the above, we obtain the sum of squares

and note that this is maximized in absolute value when all the signs line up:

where we have used Cauchy-Schwarz and the fact max⁡0⩽α⩽1α(1−α)=(1/2)2\max_{0\leqslant\alpha\leqslant 1}\alpha(1-\alpha)=(1/2)^{2}. The other terms are all similar:

6 Repairing Almost-Pseudo-Distributions

λmin⁡L⁡=−ε\operatorname{\lambda_{min}}\operatorname{\mathcal{L}}=-\varepsilon.

is a valid pseudo-expectation satisfying {∥x∥2=1}\{\|x\|^{2}=1\}.

7 Putting Everything Together

We are ready to prove Theorem 6.3 and Theorem 6.4. The proof of Theorem 6.3 is somewhat simpler and contains many of the ideas of the proof of Theorem 6.4, so we start there.

where we recall δ2\delta_{2} and δ2′\delta_{2}^{\prime} defined in the theorem statement. Finally, for ξ2′\xi_{2}^{\prime}, we have

7.2 The Degree-3 Lower Bound

The functional L⁡\operatorname{\mathcal{L}} contains our current best guess at the degree 1 and 2 moments of a pseudo-distribution whose degree-3 moments are ε\varepsilon-correlated with A(x)\mathbf{A}(x).

The next step is to use symmetric Schur complement to extend L⁡\operatorname{\mathcal{L}} to a degree-44 pseudo-expectation. Note that ML⁡∣3M_{\operatorname{\mathcal{L}}|_{3}} decomposes as

Since we have the same assumptions on AπA^{\pi} for all π∈S3\pi\in\mathcal{S}_{3}, without loss of generality we analyze just the case that π\pi is the identity permutation, in which case Aπ=AA^{\pi}=A.

Here we have used Corollary 6.7 and Corollary 6.6 to express a general element of Sym⁡(ε2n2ΠAATΠ)\operatorname{Sym}(\frac{\varepsilon^{2}}{n^{2}}\Pi AA^{T}\Pi) in terms of Π,AAT,σ⋅AAT\Pi,AA^{T},\sigma\cdot AA^{T}, and σ2⋅AAT\sigma^{2}\cdot AA^{T}.

Finally, our assumptions on ⟨A,∑π∈S3Aπ⟩\langle A,\sum_{\pi\in\mathcal{S}_{3}}A^{\pi}\rangle yield

where δ1\delta_{1} is as defined in the theorem statement.

Now using Lemma 6.15, we can correct the negative eigenvalue of L⁡1\operatorname{\mathcal{L}}^{1} to get a pseudo-expectation

Higher-Order Tensors

We have heretofore restricted ourselves to the case k=3k=3 in our algorithms for the sake of readability. In this section we state versions of our main results for general kk and indicate how the proofs from the 33-tensor case may be generalized to handle arbitrary kk. Our policy is to continue to treat kk as constant with respect to nn, hiding multiplicative losses in kk in our asymptotic notation.

For even kk, the degree-kk SoS approach does not improve on the tensor unfolding algorithms of Montanari and Richard [MR14]. Indeed, by performing a similar variable substitution, yβ=xβy_{\beta}=x^{\beta} for all ∣β∣=k/2|\beta|=k/2, the SoS algorithm reduces exactly to the eigenvalue/eigenvector computation from tensor unfolding. If we perform instead the substitution yβ=xβy_{\beta}=x^{\beta} for ∣β∣=k/2−1|\beta|=k/2-1, it becomes possible to extract v0v_{0} directly from the degree-22 pseudo-moments of an (approximately) optimal degree-44 pseudo-distribution, rather than performing an extra step to recover v0v_{0} from vv well-correlated with v0⊗k/2v_{0}^{\otimes k/2}. Either approach recovers v0v_{0} only up to sign, since the input is unchanged under the transformation v0↦−v0v_{0}\mapsto-v_{0}.

We now state analogues of all our results for general kk. Except for the above noted differences from the k=3k=3 case, the proofs are all easy transformations of the proofs of their degree-33 counterparts.

There is an algorithm, based on semidefinite programming, which on input T(x)=τ⋅⟨v0,x⟩k+A(x)\mathbf{T}(x)=\tau\cdot\langle v_{0},x\rangle^{k}+\mathbf{A}(x) returns a unit vector vv with ⟨v0,v⟩⩾1−ε\langle v_{0},v\rangle\geqslant 1-\varepsilon with high probability over random choice of A\mathbf{A}.

There is an algorithm, based on semidefinite programming, which on input T(x)=τ⋅⟨v0,x⟩k+A(x)\mathbf{T}(x)=\tau\cdot\langle v_{0},x\rangle^{k}+\mathbf{A}(x) certifies that T(x)⩽τ⋅⟨v,x⟩k+O(nk/4log⁡(n)1/4)\mathbf{T}(x)\leqslant\tau\cdot\langle v,x\rangle^{k}+O(n^{k/4}\log(n)^{1/4}) for some unit vv with high probability over random choice of A\mathbf{A}. This guarantees in particular that vv is close to a maximum likelihood estimator for the problem of recovering the signal v0v_{0} from the input τ⋅v0⊗k+A\tau\cdot v_{0}^{\otimes k}+\mathbf{A}.

For even kk, the above all hold, except now we recover vv with ⟨v0,v⟩2⩾1−ε\langle v_{0},v\rangle^{2}\geqslant 1-\varepsilon, and the algorithms can be implemented in nearly-linear time.

The next theorem partially resolves a conjecture of Montanari and Richard regarding tensor unfolding algorithms for odd kk. We are able to prove their conjectured signal-to-noise ratio τ\tau, but under an asymmetric noise model. They conjecture that the following holds when A\mathbf{A} is symmetric with unit Gaussian entries.

Conclusion

One theme in this work has been efficiently certifying upper bounds on homogeneous polynomials with random coefficients. It is an interesting question to see whether one can (perhaps with the degree d>4d>4 SoS meta-algorithm) give an algorithm certifying a bound of n3/4−δn^{3/4-\delta} over the unit sphere on a degree 33 polynomial with standard Gaussian coefficients. Such an algorithm would likely yield improved signal-to-noise guarantees for tensor PCA, and would be of interest in its own right.

Conversely, another problem is to extend our lower bound to handle degree d>4d>4 SoS. Together, these two problems suggest (as was independently suggested to us by Boaz Barak) the problem of characterizing the SoS degree required to certify a bound of n3/4−δn^{3/4-\delta} as above.

Another problem is to simplify the linear time algorithm we give for tensor PCA under symmetric noise. Montanari and Richard’s conjecture can be interpreted to say that the random rotations and decomposition into submatrices involved in our algorithm are unnecessary, and that in fact our linear time algorithm for recovery under asymmetric noise actually succeeds in the symmetric case.

Acknowledgments

We thank Moses Charikar for bringing to our attention the work of Montanari and Richard. We would like to thank Boaz Barak, Rong Ge, and Ankur Moitra for enlightening conversations. S. B. H. acknowledges the support of an NSF Graduate Research Fellowship under award no. 1144153. D. S. acknowledges support from the Simons Foundation, the National Science Foundation, an Alfred P. Sloan Fellowship, and a Microsoft Research Faculty Fellowship, A large portion of this work was completed while the authors were long-term visitors to the Simons Institute for the Theory of Computing (Berkeley) for the program on Algorithmic Spectral Graph Theory.

References

Appendix A Pseudo-Distribution Facts

Let x,yx,y be vector-valued polynomials. Then

Let x,yx,y be vector-valued polynomials and d>0d>0 an integer. Then

Note that ⟨x,y⟩d=⟨x⊗d,y⊗d⟩\langle x,y\rangle^{d}=\langle x^{\otimes d},y^{\otimes d}\rangle and apply Lemma A.2. ∎

Yet another version of pseudo-Cauchy-Schwarz will be useful:

Let {x,y}\{x,y\} be a degree dd pseudo-distribution over a pair of vectors, d⩾2d\geqslant 2. Then

Again, see [BKS14b] for the cleanest proof.

Let p(u)p(u) be the univariate polynomial p(u)=1−2u3+up(u)=1-2u^{3}+u. It is easy to check that p(u)⩾0p(u)\geqslant 0 for u∈u\in. It follows from classical results about univariate polynomials that p(u)p(u) then can be written as

for some SoS polynomials s0,s1,s2s_{0},s_{1},s_{2} of degrees at most 22. (See [OZ13], fact 3.2 for a precise statement and attributions.)

We have by Lemma A.2 that ⟨x,v0⟩⪯12(∥x∥2+1)\langle x,v_{0}\rangle\preceq\frac{1}{2}(\|x\|^{2}+1) and also that ⟨x,v0⟩⪰−12(∥x∥2+1)\langle x,v_{0}\rangle\succeq-\frac{1}{2}(\|x\|^{2}+1). Multiplying the latter SoS relation by the SoS polynomial s1(⟨x,v0⟩)s_{1}(\langle x,v_{0}\rangle) and the former by s2(⟨x,v0⟩)s_{2}(\langle x,v_{0}\rangle), we get that

where in the second-to-last step we have used the assumption that {x}\{x\} satisfies {∥x∥2=1}\{\|x\|^{2}=1\}. A similar analysis yields

We will need a bound on the pseudo-expectation of a degree-33 polynomial in terms of the operator norm of its coefficient matrix.

We begin by expanding in the monomial basis and using pseudo-Cauchy-Schwarz:

Appendix B Concentration bounds

We will be extensively concerned with various real random matrices. A great deal is known about natural classes of such matrices; see the excellent book of Tao [Tao12] and the notes by Vershynin and Tropp [Ver11, Tro12].

It will be convenient to use the following standard result on the concentration of empirical covariance matrices. This statement is borrowed from [Ver11], Corollary 5.50.

We will also need the matrix Bernstein inequality. This statement is borrowed from Theorem 1.6.2 of Tropp [Tro12].

We will need bounds on the operator norm of random square rectangular matrices, both of which are special cases of Theorem 5.39 in [Ver11].

Let AA be an n×nn\times n matrix with independent entries from N(0,1)\mathcal{N}(0,1). Then with probability 1−n−ω(1)1-n^{-\omega(1)}, the operator norm ∥A∥\|A\| satisfies ∥A∥⩽O(n)\|A\|\leqslant O(\sqrt{n}).

Let AA be an n2×nn^{2}\times n matrix with independent entries from N(0,1)\mathcal{N}(0,1). Then with probability 1−n−ω(1)1-n^{-\omega(1)}, the operator norm ∥A∥\|A\| satisfies ∥A∥⩽O(n)\|A\|\leqslant O(n).

Our first concentration theorem provides control over the nontrivial permutations of the matrix AATAA^{T} under the action of S4\mathcal{S}_{4} for a tensor A\mathbf{A} with independent entries.

Let c∈{1,2}c\in\{1,2\} and d⩾1d\geqslant 1 an integer. Let A1,…,AncA_{1},\ldots,A_{n^{c}} be iid random matrices in {±1}nd×nd\{\pm 1\}^{n^{d}\times n^{d}} or with independent entries from N(0,1)\mathcal{N}(0,1). Then, with probability 1−O(n−100)1-O(n^{-100}),

We can prove Theorem 3.3 as a corollary of the above.

Now by Theorem B.5, we know that for AiA_{i} the slices of the tensor A\mathbf{A} from the statement of Theorem 3.3,

Now we prove Theorem B.5. We will prove only the statement about ∑iAi⊗Ai\sum_{i}A_{i}\otimes A_{i}, as the case of ∑iAi⊗AiT\sum_{i}A_{i}\otimes A_{i}^{T} is similar.

Let A1,…,AncA_{1},\ldots,A_{n^{c}} be as in Theorem B.5. We first need to get a handle on their norms individually, for which we need the following lemma.

Let AA be a random matrix in {±1}nd×nd\{\pm 1\}^{n^{d}\times n^{d}} or with independent entries from N(0,1)\mathcal{N}(0,1). For all t⩾1t\geqslant 1, the probability of the event {∥A∥>tnd/2}\{\lVert A\rVert>tn^{d/2}\} is at most 2−t2nd/K2^{-t^{2}{n^{d}}/K} for some absolute constant KK.

The subgaussian norm of the rows of AA is constant and they are identically and isotropically distributed. Hence Theorem 5.39 of [Ver11] applies to give the result. ∎

Since the norms of the matrices A1,…,AncA_{1},\ldots,A_{n^{c}} are concentrated around nd/2n^{d/2} (by Lemma B.6), it will be enough to prove Theorem B.5 after truncating the matrices A1,…,AncA_{1},\ldots,A_{n^{c}}. For t⩾1t\geqslant 1, define iid random matrices A1′,…,Anc′A^{\prime}_{1},\ldots,A^{\prime}_{n^{c}} such that

for some tt to be chosen later. Lemma B.6 allows us to show that the random matrices Ai⊗AiA_{i}\otimes A_{i} and Ai′⊗Ai′A^{\prime}_{i}\otimes A^{\prime}_{i} have almost the same expectation. For the remainder of this section, let KK be the absolute constant from Lemma B.6.

For every i∈[nc]i\in[n^{c}] and all t⩾1t\geqslant 1, the expectations of Ai⊗AiA_{i}\otimes A_{i} and Ai′⊗Ai′A^{\prime}_{i}\otimes A^{\prime}_{i} satisfy

Using Jensen’s inequality and that Ai=Ai′A_{i}=A_{i}^{\prime} unless ∥Ai∥>tnd/2\|A_{i}\|>tn^{d/2}, we have

For R=2t2ndR=2t^{2}n^{d}, the random matrices B1′,…,Bnc′B^{\prime}_{1},\ldots,B^{\prime}_{n^{c}} satisfy {∥Bi′∥⩽R}\{\lVert B^{\prime}_{i}\rVert\leqslant R\} with probability 11. Therefore, by the Bernstein bound for non-symmetric matrices [Tro12, Theorem 1.6],

Since our parameters satisfy t2C⋅n(4d+c)/2/3⩽t4n(2d+c)t^{2}C\cdot n^{(4d+c)/2}/3\leqslant t^{4}n^{(2d+c)}, this probability is bounded by

At this point, we have all components of the proof of Theorem B.5.

At the same time, by Lemma B.6 and a union bound,

We choose t=1t=1 and C=1002Kdlog⁡nC=100\sqrt{2Kd\log n} and assume that nn is large enough so that C⋅n(2d+c)/2⩾nc⋅2−tnd/KC\cdot n^{(2d+c)/2}\geqslant n^{c}\cdot 2^{-t{n^{d}}/K} and 2n2d⋅exp⁡(−C2Kt4)⩾nc⋅2−t2nd/K2n^{2d}\cdot\exp\left(\frac{-C^{2}}{Kt^{4}}\right)\geqslant n^{c}\cdot 2^{-t^{2}{n^{d}}/K}. Then the probability satisfies

B.3 Concentration for Spectral SoS Analyses

The first claim is immediate from Theorem B.5. For the second, we note that since v0v_{0} is a unit vector, the matrix ∑iv0(i)Ai\sum_{i}v_{0}(i)A_{i} has independent entries from N(0,1)\mathcal{N}(0,1). Thus, by Lemma B.3, ∥∑iv0(i)Ai∥⩽O(n)\|\sum_{i}v_{0}(i)A_{i}\|\leqslant O(\sqrt{n}) with probability 1−O(n−100)1-O(n^{-100}), as desired. ∎

With δ=O(1/n)\delta=O(1/\sqrt{n}) and t=1t=1, our parameters will satisfy n(k+1)/2⩾(t/δ)2n(k−1)/2n^{(k+1)/2}\geqslant(t/\delta)^{2}n^{(k-1)/2}. Hence, by Lemma B.1,

with probability at least 1−2exp⁡(−n(k+1)/2)⩾1−O(n−100)1-2\exp(-n^{(k+1)/2})\geqslant 1-O(n^{-100}).

B.4 Concentration for Lower Bounds

The next theorems collects the concentration results necessary to apply our lower bounds Theorem 6.3 and Theorem 6.4 to random polynomials.

For (B.11), from Theorem B.5, Lemma 6.8, the observation that multiplication by an orthogonal operator cannot increase the operator norm, a union bound over all π\pi, and the triangle inequality, it follows that:

with probability 1−n−1001-n^{-100}. By the definition of the operator norm and another application of triangle inequality, this implies

We turn to (B.2). By a Chernoff bound, ⟨A,A⟩=Ω(n3)\langle\mathbf{A},\mathbf{A}\rangle=\Omega(n^{3}) with probability 1−n−1001-n^{-100}. Let π∈S3\pi\in\mathcal{S}_{3} be a nontrivial permutation. To each multi-index α\alpha with ∣α∣=3|\alpha|=3 we associate its orbit Oα\mathcal{O}_{\alpha} under ⟨π⟩\langle\pi\rangle. If α\alpha has three distinct indices, then ∣Oα∣>1|\mathcal{O}_{\alpha}|>1 and ∑β∈OαAβAβπ\sum_{\beta\in\mathcal{O}_{\alpha}}A_{\beta}A^{\pi}_{\beta} is a random variable XαX_{\alpha} with the following properties:

∣Xα∣<O(log⁡n)|X_{\alpha}|<O(\log n) with probability 1−n−ω(1)1-n^{-\omega(1)}.

XαX_{\alpha} and −Xα-X_{\alpha} are identically distributed.

We have Πw=w\Pi w=w and we let eij:=Π(ei⊗ej)=12(ei⊗ej+ej⊗ei)e_{ij}\mathrel{\mathop{:}}=\Pi(e_{i}\otimes e_{j})=\frac{1}{2}(e_{i}\otimes e_{j}+e_{j}\otimes e_{i}). So using Lemma 6.8,

For i≠ji\neq j, each term wT(Akej⊗Akei)w^{T}(A_{k}e_{j}\otimes A_{k}e_{i}) (or similar, with various transposes) is the sum of nn independent products of pairs of independent unit Gaussians, so by a Chernoff bound followed by a union bound, with probability 1−n−ω(1)1-n^{-\omega(1)} all of them are O(nlog⁡n)O(\sqrt{n}\log n). There are O(n)O(n) such terms, for an upper bound of O(n3/2(log⁡n))O(n^{3/2}(\log n)) on the contribution from the tensored parts.

At the same time, wTAw^{T}A is a sum ∑kakk\sum_{k}a_{kk} of nn rows of AA and AeijAe_{ij} is the average of two rows of AA; since i≠ji\neq j these rows are independent from wTAw^{T}A. Writing this out, wTAATeij=12∑k⟨akk,aij+aji⟩w^{T}AA^{T}e_{ij}=\frac{1}{2}\sum_{k}\langle a_{kk},a_{ij}+a_{ji}\rangle. Again by a standard Chernoff and union bound argument this is in absolute value at most O(n3/2(log⁡n))O(n^{3/2}(\log n)) with probability 1−n−ω(1)1-n^{-\omega(1)}. In sum, when i≠ji\neq j, with probability at least 1−n−ω(1)1-n^{-\omega(1)}, we get ∣L⁡∥x∥2xixj∣=O(1/n2log⁡n)|\operatorname{\mathcal{L}}\|x\|^{2}x_{i}x_{j}|=O(1/n^{2}\log n). After a union bound, the maximum over all i,ji,j is O(1/n2)O(1/n^{2}). This concludes (B.5).

In the i=ji=j case, since ∑k⟨w,Akei⊗Akei⟩=∑j,k⟨ej,Akei⟩2\sum_{k}\langle w,A_{k}e_{i}\otimes A_{k}e_{i}\rangle=\sum_{j,k}\langle e_{j},A_{k}e_{i}\rangle^{2} is a sum of n2n^{2} independent square Gaussians, by a Bernstein inequality, ∣∑k⟨w,Akei⊗Akei⟩−n2∣⩽O(nlog⁡1/2n)|\sum_{k}\langle w,A_{k}e_{i}\otimes A_{k}e_{i}\rangle-n^{2}|\leqslant O(n\log^{1/2}n) with probability 1−n−ω(1)1-n^{-\omega(1)}. The same holds for the other tensored terms, and for wTAATeiiw^{T}AA^{T}e_{ii}, so when i=ji=j we get that ∣O(λ2)L⁡∥x∥2xi2−5∣⩽O((log⁡1/2n)/n)|O(\lambda^{2})\operatorname{\mathcal{L}}\|x\|^{2}x_{i}^{2}-5|\leqslant O((\log^{1/2}n)/n) with probability 1−n−ω(1)1-n^{-\omega(1)}. Summing over all ii, we find that ∣O(λ2)L⁡∥x∥4−5n∣⩽O(log⁡1/2n)|O(\lambda^{2})\operatorname{\mathcal{L}}\|x\|^{4}-5n|\leqslant O(\log^{1/2}n), so that O(λ2)∣L⁡∥x∥2xi2−1nL⁡∥x∥4∣⩽O((log⁡1/2n)/n)O(\lambda^{2})|\operatorname{\mathcal{L}}\|x\|^{2}x_{i}^{2}-\tfrac{1}{n}\operatorname{\mathcal{L}}\|x\|^{4}|\leqslant O((\log^{1/2}n)/n) with probability 1−n−ω(1)1-n^{-\omega(1)}. A union bound over ii completes the argument. ∎