A statistical model for tensor PCA

Introduction

Given a data matrix $\mathbf{X}$, Principal Component Analysis (PCA) can be regarded as a ‘denoising’ technique that replaces $\mathbf{X}$ by its closest rank-one approximation. This optimization problem can be solved efficiently, and its statistical properties are well-understood. The generalization of PCA to tensors is motivated by problems in which it is important to exploit higher order moments, or data elements are naturally given more than two indices. Examples include topic modeling [AGH+12], video processing, collaborative filtering in presence of temporal/context information, community detection [AGHK13], spectral hypergraph theory and hyper-graph matching [DBKP09]. Further, finding a rank-one approximation to a tensor is a bottleneck for tensor-valued optimization algorithms using conditional gradient type of schemes. While tensor factorization is NP-hard [HL13], this does not necessarily imply intractability for natural statistical models. Over the last ten years, it was repeatedly observed that either convex optimization or greedy methods yield optimal solutions to statistical problems that are intractable from a worst case perspective (well-known examples include sparse regression [DE03, Tro04, CT07] and low-rank matrix completion [CR09, KMO10]).

In order to investigate the fundamental tradeoffs between computational resources and statistical power in tensor PCA, we consider the simplest possible model where this arises, whereby an unknown unit vector $\mathbf{v_{0}}$ is to be inferred from noisy multilinear measurements. Namely, for each unordered $k$-uple $\{i_{1},i_{2},\dots,i_{k}\}\subseteq[n]$, we measure

with $\mathbf{Z}$ Gaussian noise (see below for a precise definition) and wish to reconstruct $\mathbf{v_{0}}$. In tensor notation, the observation model reads (see the end of this section for notations)

This is analogous to the so called ‘spiked covariance model’ used to study matrix PCA in high dimensions [JL09].

It is immediate to see that maximum-likelihood estimator $\mathbf{v}^{\mbox{\tiny\sc ML}}$ is given by a solution of the following problem

Solving it exactly is –in general– NP hard [HL13].

Assuming unbounded computational resources, we can solve the Tensor PCA optimization problem and hence implement the maximum likelihood estimator $\mathbf{\widehat{v}}^{\mbox{\tiny\sc ML}}$. We use recent results in probability theory to show that this approach is successful for $\beta\geq\mu_{k}$ (here $\mu_{k}$ is a constant given explicitly below, with $\mu_{k}=\sqrt{k\log k}(1+o_{k}(1))$). In particular, above this thresholdNote that, for $k$ even, $\mathbf{v_{0}}$ can only be recovered modulo sign. For the sake of simplicity, we assume here that this ambiguity is correctly resolved. we have, with high probability,

We use an information-theoretic argument to show that no approach can do significantly better, namely no procedure can estimate $\mathbf{v_{0}}$ accurately for $\beta\leq c\sqrt{k}$ (for $c$ a universal constant).

A heuristics argument suggests that the necessary and sufficient condition for tensor unfolding to succeed is indeed $\beta\gtrsim n^{(k-2)/4}$ (which is below the rigorous bound by a factor $n^{1/4}$ for $k$ odd). We can indeed confirm this conjecture for $k$ even and under an asymmetric noise model. Numerical simulations confirm the conjecture for $k=3$.

We then consider a simple tensor power iteration method, that proceeds by repeatedly applying the tensor to a vector. We prove that, initializing this iteration uniformly at random, it converges very rapidly to an accurate estimate provided $\beta\gtrsim n^{(k-1)/2}$. A heuristic argument suggests that the correct necessary and sufficient threshold is given by $\beta\gtrsim n^{(k-2)/2}$. In other words, power iteration is substantially less powerful than unfolding.

Motivated by the last observation, we consider a ‘warm-start’ power iteration algorithm, in which we initialize power iteration with the output of tensor unfolding. This approach appears to have the same threshold signal-to-noise ratio as simple unfolding, but significantly better accuracy above that threshold.

We also study a number of variations on this, with improved unfolding methods.

Finally we consider an approximate message passing (AMP) algorithm [DMM09, BM11]. Such algorithms proved effective in compressed sensing and several other estimation problems. We show that the behavior of AMP is qualitatively similar to the one of naive power iteration. In particular, AMP fails for any $\beta$ bounded as $n\to\infty$.

Given the above computational complexity barrier, it is natural to study weaker version of the original problem. Here we assume that extra information about $\mathbf{v_{0}}$ is available. This can be provided by additional measurements or by approximately solving a related problem, for instance a matrix PCA problem as in [AGH+12]. We model this additional information as ${\mathbf{y}}=\gamma\mathbf{v_{0}}+\mathbf{g}$ (with $\mathbf{g}$ an independent Gaussian noise vector), and incorporate it in the initial condition of AMP algorithm. We characterize exactly the threshold value $\gamma_{*}=\gamma_{*}(\beta)$ above which AMP converges to an accurate estimator.

The thresholds for various classes of algorithms are summarized below.

We will conclude the paper with some insights that we believe provide useful guidance for tensor factorization heuristics. We illustrate these insights through simulations.

Throughout the paper, proofs will be deferred to the Appendices.

We define the Frobenius (Euclidean) norm of a tensor $\mathbf{X}$ by $\|\mathbf{X}\|_{F}=\sqrt{\langle\mathbf{X},\mathbf{X}\rangle}$, and its operator norm by

For a permutation $\pi\in\mathfrak{S}_{k}$, we will denote by $\mathbf{X}^{\pi}$ the tensor with permuted indices $\mathbf{X}^{\pi}_{i_{1},\cdots,i_{k}}=\mathbf{X}_{\pi(i_{1}),\cdots,\pi(i_{k})}$. We call the tensor $\mathbf{X}$ symmetric if, for any permutation $\pi\in\mathfrak{S}_{k}$, $\mathbf{X}^{\pi}=\mathbf{X}$. It is proved [Wat90] that, for symmetric tensors, the value of problem Tensor PCA coincides with $\|\mathbf{X}\|_{op}$ up to a sign. More precisely, for symmetric tensors we have the equivalent representation

where $\mathbf{g}\sim{\sf N}(0,{\rm I}_{n})$, and $\mathbf{o}(1)$ is a vector with $\|\mathbf{o}(1)\|_{2}\to 0$ in probability as $n\to\infty$. We further have,

with $g\sim{\sf N}(0,1)$. Finally notice that, for $k$ even, in Spiked Tensor Model, the vector $\mathbf{v_{0}}$ can always be recovered up to a sign flip. This suggest the use of the loss function

Ideal estimation

In this section we consider the problem of estimating $\mathbf{v_{0}}$ under the Spiked Tensor Model, when no constraint is imposed on the complexity of the estimator. Our first result is a lower bound on the loss of any estimator.

In order to establish a matching upper bound on the loss, we consider the maximum likelihood estimator $\mathbf{\widehat{v}}^{\mbox{\tiny\sc ML}}$, obtained by solving the Tensor PCA problem. As in the case of matrix denoising, we expect the properties of this estimator to depend on signal to noise ratio $\beta$, and on the ‘norm’ of the noise $\|\mathbf{Z}\|_{op}$ (i.e. on the value of the optimization problem Tensor PCA in the case $\beta=0$). For the matrix case $k=2$, this coincides with the largest eigenvalue of $\mathbf{Z}$. Classical random matrix theory shows that –in this case– $\|\mathbf{Z}\|_{op}$ concentrates tightly around $2$ [Gem80, DS01a, BS10].

It turns out that tight results for $k\geq 3$ follow immediately from a technically sophisticated analysis of the stationary points of random Morse functions by Auffinger, Ben Arous and Cerny [ABAC13]. (See Appendix B.1 for further background.)

There exists a sequence of real numbers $\{\mu_{k}\}_{k\geq 2}$, such that

Further $\|\mathbf{Z}\|_{op}$ concentrates tightly around its expectation. Namely, for any $n,k$

Finally $\mu_{k}=\sqrt{k\log k}(1+o_{k}(1))$ for large $k$.

An explicit expression for the quantity $\mu_{k}$ is given in Appendix B (which also contains a proof, that uses [ABAC13]). Evaluating this expression for small values of $k$, we get the following explicit values, that we also compare with the large-$k$ asymptotics $\sqrt{k\log k}$. (It is not hard to increase the number of digits in these evaluations, using the expressions in Appendix.)

For instance, this table indicates that a large order-$3$ Gaussian tensor should have $\|\mathbf{Z}\|_{op}\approx 2.87$, while a large order $10$ tensor has $\|\mathbf{Z}\|_{op}\approx 6.75$. As a simple consequence of Lemma 2.1, we establish an upper bound on the error incurred by the maximum likelihood estimator, see Section B.2 for a proof.

Let $\mu_{k}$ be the sequence of real numbers introduced above. Letting $\mathbf{\widehat{v}}^{\mbox{\tiny\sc ML}}$ denote the maximum likelihood estimator (i.e. the solution of Tensor PCA), we have for $n$ large enough, and all $s>0$

with probability at least $1-2e^{-ns^{2}/(16k)}$.

The following upper bound on the value of the problem Tensor PCA is proved using Sudakov-Fernique inequality. While it is looser than Lemma 2.1 (corresponding to the case $\beta=0$), we expect it to become sharp for $\beta\geq\beta_{k}$ a suitably large constant. We refer to Appendix B.3 for its proof.

The most striking prediction from statistical physics is that the function $H_{\mathbf{Z}}(\mathbf{v})$ has an exponential number of local maxima on the unit sphere [CS95]. Furthermore, there exists $\eta_{k}<\mu_{k}$ such that, for each $x\in[\eta_{k},\mu_{k})$ the number of local maxima with value $H_{\mathbf{Z}}(\mathbf{v})\approx x$ is $\exp\{\Theta(n)\}$. In [ABAC13] rigorous evidence is developed to this support picture.

In the Spiked Tensor Model these local maxima translate into undesired local maxima of $\langle\mathbf{X},\mathbf{v}{{}^{\otimes k}}\rangle$. It is natural to guess that these local maxima affect local iterative algorithms, and that these do not converge to a good estimate of $\mathbf{v_{0}}$ unless they are initialized within thee ‘basin of attraction’ of $\mathbf{v_{0}}$. The analysis in the next sections confirms this intuition.

Tensor Unfolding

A simple and popular heuristics to obtain tractable estimators of $\mathbf{v_{0}}$ consists in constructing a suitable matrix with the entries of $\mathbf{X}$, and performing principal component analysis on this matrix. Since the number of distinct entries of $\mathbf{X}$ is of order $n^{k}$, the resulting matrix ${\sf Mat}_{q}(\mathbf{X})$ has dimension $\Theta(n^{q})\times\Theta(n^{k-q})$. This operation is variously referred as matricization, unfolding, flattening. While the details of this construction can vary, we do not expect them to affect qualitatively our results, that we summarize for the sake of convenience:

The best way to unfold $\mathbf{X}$ amounts to form a matrix as square as possible.

Setting $b=(\lceil k/2\rceil-1)/2$ (in particular $b=1/2$ for $k\in\{3,4\}$), the unfolding approach succeeds when $\beta$ is larger than $n^{b}$. This is to be compared with $\beta=\Theta(1)$ that is sufficient for the maximum likelihood estimator (see previous section).

Based on heuristic arguments, we believe that the tight threshold is $\beta\gtrsim n^{(k-2)/4}$ (i.e. that $\beta\gtrsim n^{(k-2)/4}$ is both necessary and sufficient –modulo constants).

A sharper analysis is possible when the symmetric noise tensor $\mathbf{Z}$ in our Spiked Tensor Model is replaced by non-symmetric Gaussian noise, and $k$ is even. In particular, we can confirm the above conjecture in this case. (As mentioned, we expect similar results to hold more generally.)

In this case, if $\beta\leq(1-\varepsilon)n^{b}$, then the estimator from unfolding is essentially orthogonal to the signal $\mathbf{v_{0}}$. On the other hand, if $\beta\geq(1+\varepsilon)n^{b}$, we construct an estimator with $|\langle\mathbf{\widehat{v}},\mathbf{v_{0}}\rangle|\to 1$.

We achieves the remarkable behavior at the last point by a recursive unfolding method. In a nutshell we perform principal component analysis on ${\sf Mat}_{q}(\mathbf{X})$, construct a matrix out of the principal vector, and then perform again principal component analysis.

Standard convex relaxations of low-rank tensor estimation problem compute factorizations of ${\sf Mat}_{q}(\mathbf{X})$[TSHK11, LMWY13, MHG13, RPP13]. Not all unfoldings (choices of $q$) are equivalent. It is natural to expect that this approach will be successful only if the signal-to-noise ratio exceeds the operator norm of the unfolded noise $\|{\sf Mat}_{q}(\mathbf{Z})\|_{op}$. The next lemma suggests that the latter is minimal when ${\sf Mat}_{q}(\mathbf{Z})$ is ‘as square as possible’ . A similar phenomenon was observed in a different context in [MHG13].

For any integer $0\leq q\leq k$ we have, for some universal constant $C_{k}$,

For all $n$ large enough, both bounds are minimized for $q=\lceil k/2\rceil$. Further

is a Lipschitz function of the Gaussian vector $\mathbf{G}$ with modulus at most $\sqrt{k/n}$. Hence the same holds for $\|{\sf Mat}_{q}(\mathbf{Z})\|_{op}=\max_{\mathbf{u},\mathbf{v}}\langle\mathbf{u},{\sf Mat}_{q}(\mathbf{Z})\mathbf{v}\rangle$, and the claim follows from Gaussian concentration of measure.

For the upper bound in Eq. (21), note that

Let us recall the following standard result derived directly from Wedin perturbation Theorem [Wed72], and stated in the context of the spiked model.

Note $\beta>0$ is the only singular value of $\beta\mathbf{u_{0}}\mathbf{w_{0}}^{\sf T}$, while the second singular value of $(\beta\mathbf{u_{0}}\mathbf{w_{0}}^{\sf T}+{\bf\Xi})$ is at most $\|{\bf\Xi}\|_{op}$. Wedin Theorem states that, for all $\beta>\|{\bf\Xi}\|_{op}$, we have

In particular $|\sin(\widehat{\mathbf{w}},\mathbf{w_{0}})|\leq 2\|{\bf\Xi}\|_{op}/\beta$ for $\beta\geq 2\|{\bf\Xi}\|_{op}$. Hence the claim (26) follows from

Letting $\mathbf{w}=\mathbf{w}(\mathbf{X})$ denote the top right singular vector of ${\sf Mat}(\mathbf{X})$, we have the following, for some universal constant $C=C_{k}>0$, and $b\equiv(1/2)(\lceil k/2\rceil-1)$.

If $\beta\geq 5\,k^{1/2}\,n^{b}$ then, with probability at least $1-n^{-2}$, we have

where $\mathbf{u_{0}}={\sf{vec}}(\mathbf{v_{0}}^{\otimes\lfloor k/2\rfloor})$, $\mathbf{w_{0}}={\sf{vec}}(\mathbf{v_{0}}^{\otimes\lceil k/2\rceil})$. We know by Lemma 3.1 that $\|{\sf Mat}(\mathbf{Z})\|_{op}\leq(5/2)\sqrt{k}\,n^{b}$ with the claimed probability. The loss upper bound (29) follows immediately from this upper bound and Wedin’s theorem Eq. (26). ∎

2 Asymmetric noise and recursive unfolding

A technical complication in analyzing the random matrix ${\sf Mat}_{q}(\mathbf{X})$ lies in the fact that its entries are not independent, because the noise tensor $\mathbf{Z}$ is assumed to be symmetric. In the next theorem we consider the case of non-symmetric noise and even $k$. This allows us to leverage upon known results in random matrix theory [Pau07, FP09, BGN12] to obtain: $(i)$ Asymptotically sharp estimates on the critical signal-to-noise ratio; $(ii)$ A lower bound on the loss below the critical signal-to-noise ratio. Namely, we consider observations

In other words $\mathbf{w}(\mathbf{\widetilde{X}})$ is a good estimate of $\mathbf{v_{0}}^{\otimes(k/2)}$ if and only if $\beta$ is larger than $n^{b}$.

we then let $\mathbf{\widehat{v}}$ to be the left principal vector of ${\sf Mat}_{1}(\mathbf{X})$. We refer to this algorithmIn practice int might be more effective to use a balanced matricization at the second step. For instance if $k$ is a power of two one could construct a square matricization and repeat the same process. For analysis purposes, we prefer the version described here. as to recursive unfolding.

Let $\mathbf{\widetilde{X}}$ be distributed according to the non-symmetric model (31) with $k\geq 4$ even, define $b\equiv(k-2)/4$. and let $\mathbf{\widehat{v}}$ be the estimate obtained by two-steps recursive unfolding.

If $\beta\geq(1+\varepsilon)n^{b}$ then, almost surely

For the sake of simplicity, we assume $\beta/n^{b}\to\varepsilon$. The limit along other sequences follows from a standard subsequence argument.

It follows from the invariance of the noise distribution in Eq. (31) that

where $\mathbf{g}\sim{\sf N}(0,{\rm I}_{n^{k/2}})$. It follows from Eq. (33), together with the almost sure limits $\lim_{n\to\infty}\|\mathbf{g}\|_{2}/n^{k/4}=1$ and $\lim_{n\to\infty}\langle\mathbf{g},{\sf{vec}}(\mathbf{v_{0}}^{\otimes(k/2)})\rangle/n^{k/4}=1$ that (almost surely)

Using the definition (34), we then have (recall that $b=(k-2)/4$)

Since $\rho_{n}$ is bounded away from zero as $n\to\infty$, Wedin’s theorem implies $\lim_{n\to\infty}|\langle\mathbf{\widehat{v}},\mathbf{v_{0}}\rangle|=1$, and therefore the claim (35). ∎

We conjecture that the weaker condition $n\gtrsim n^{(k-2)/4}$ is indeed sufficient also for our original symmetric noise model model, both for $k$ even and for $k$ odd.

Power Iteration

Iterating over (multi-) linear maps induced by a (tensor) matrix is a standard method for finding leading eigenpairs, see [KM11] and references therein for tensor-related results. In this section we will consider a simple power iteration, and then its possible uses in conjunction with tensor unfolding. Finally, we will compare our analysis with results available in the literature.

Approximate Message Passing (AMP) provides a different iterative strategy and will be discussed in Section 5. While the qualitative behavior is the same as for naive power iteration, a sharper asymptotic analysis is possible for AMP.

The simplest iterative approach is defined by the following recursion

The following result establishes convergence criteria for this iteration, first for generic noise $\mathbf{Z}$ and then for standard normal noise (using Lemma 2.1).

Then for all $t\geq t_{0}(k)$, the power iteration estimator satisfies

If $\mathbf{Z}$ is a standard normal noise tensor, then conditions (41), (41) are satisfied with high probability provided

We next discuss two aspects of this result: $(i)$ The requirement of a positive correlation between initialization and ground truth ; $(ii)$ Possible scenarios under which the assumptions of Theorem 6 are satisfied.

Notice that we require a positive correlation of the initialization ${\mathbf{y}}$ with the ground truth $\mathbf{v_{0}}$. The underlying reason is that, if $\langle\mathbf{v}^{0},\mathbf{v_{0}}\rangle$ is small, then $\langle\mathbf{v}^{t},\mathbf{v_{0}}\rangle$ remains small at all subsequent iterations. In order to clarify this point, it is instructive to compute the distribution of $\mathbf{v}^{1}$ for standard Gaussian noise $\mathbf{Z}$. We let

Using Eq. (9) and the fact that $\mathbf{v}^{0}$ is independent of $\mathbf{Z}$, we get

where $\mathbf{g}\sim{\sf N}(0,{\rm I}_{n})$, and $\mathbf{o}(1)$ is a vector with $\|\mathbf{o}(1)\|_{2}\to 0$ in probability as $n\to\infty$. In particular

In particular $\tau^{1}\lesssim\tau^{0}$ only if $\beta\tau_{0}^{k-2}\gtrsim 1$, or, equivalently, $\langle{\mathbf{y}},\mathbf{v_{0}}\rangle/\|{\mathbf{y}}\|_{2}\gtrsim\beta^{-1/(k-2)}$. This suggest that the condition in Eq. (45) is not too far from being tight (in the sense that the exponent $-1/(k-1)$ can at best replaced by $-1/(k-2)$).

In general we cannot assume that an initialization satisfying the conditions of Theorem 6. Hence, unlike for ordinary matrix factorization, power iteration is not a practical solution to the tensor principal component problem. There are however circumstances under which a sufficiently good initialization exists.

If $y$ is a uniformly random vector on the unit sphere, then $\langle\mathbf{v_{0}},{\mathbf{y}}\rangle$ is approximately normal with mean zero and variance $1/n$. For instance $|\langle\mathbf{v_{0}},{\mathbf{y}}\rangle|\geq 1/\sqrt{n}$ with probability roughly $0.32$.

Comparing this with condition (42), we obtain that a random initialization succeed with positive probability if

For standard Gaussian noise, this amounts to requiring $\beta\geq(2n)^{(k-1)/2}\mu_{k}$. The above heuristic analysis suggests that the correct condition should be $\beta\gtrsim n^{(k-2)/2}$.

Additional information might be available about the vector $\mathbf{v_{0}}$. This information can be used for initiating the power iteration. In the next section we consider the special case in which tensor unfolding is used for initializing power iteration.

2 Comparison with Tensor Unfolding

It is instructive to compare the result of the previous section with the ones for tensor unfolding, cf. Section 3. Summarizing, for standard Gaussian noise

Tensor unfolding is guaranteed to succeed provided $\beta\gtrsim n^{b}$, with $b=(\lceil k/2\rceil-1)/2$. We conjecture that a necessary and sufficient condition is in fact $\beta\gtrsim n^{(k-2)/4}$ (e.g. $\beta\gtrsim n^{1/4}$ for order $3$ tensors).

Power iteration, with random initialization requires $\beta\gtrsim n^{(k-1)/2}$. Our heuristic calculation suggests that a necessary and sufficient condition is in fact $\beta\gtrsim n^{(k-2)/2}$ (e.g. $\beta\gtrsim n^{1/2}$ for order $3$ tensors)..

In other words, tensor unfolding is successful under a signal-to-noise ratio that is order of magnitudes smaller than power iteration. This suggests the following warm start procedure: $(i)$ Compute a first estimate $\mathbf{\widehat{v}}^{\rm Unfold}$ of $\mathbf{v_{0}}$ using tensor unfolding; $(ii)$ Use this as initialization for the power iteration, hence setting $\mathbf{v}^{0}=\mathbf{\widehat{v}}^{\rm Unfold}$. We will explore this approach numerically in Section 6.

3 Related work

As mentioned above, power iteration is a natural approach to tensor factorization and was studied in several earlier papers. Most recently, interest within machine learning was spurred by [AGH+12].

Our Theorem 6 is analogous to the main result of [AGH+12] although incomparable:

In [AGH+12] the ‘signal’ part of the tensor $\mathbf{X}$ is assumed to have an orthogonal decomposition $\sum_{i=1}^{n}\lambda_{i}\mathbf{v}_{i}^{\otimes k}$ with $\min_{i}(\lambda_{i})$ bounded away from zero. Here, the signal part has rank one (equivalently, all the $\lambda_{i}$’s but one vanish).

In [AGH+12] only the case of third order tensors ($k=3$) is considered. We characterize power iteration for general $k$.

We establish convergence in a number of iterations $t$ that is independent of the dimensions $n$. In [AGH+12] the number of iterations is bounded by a polynomial in $n$.

We evaluate our bounds in the case of Gaussian noise. This allows a comparison with other methods, such as tensor unfolding.

Asymptotics via Approximate Message Passing

Approximate message passing (AMP) algorithms [DMM09, BM11] proved successful in several high-dimensional estimation problems including compressed sensing, low rank matrix reconstruction, and phase retrieval [FRVB11, KRFU12, SC11, SR12]. An appealing feature of this class of algorithms is that their high-dimensional limit can be characterized exactly through a technique known as ‘state evolution.’ Here we develop an AMP algorithm for tensor data, and its state evolution analysis focusing on the fixed $\beta$, $n\to\infty$ limit. Proofs follows the approach of [BM11] and will be presented in a journal publication.

(Note that, unlike in power iteration, we normalize $\mathbf{v}^{t}$ ‘before’ multiplying it by $\mathbf{X}$. This choice is equivalent but yields slightly simpler expression.)

Our main conclusion is that the behavior of AMP is qualitatively similar to the one of power iteration. However, we can establish stronger results in two respects:

We can prove that, unless side information is provided about the signal $\mathbf{v_{0}}$, the AMP estimates remains essentially orthogonal to $\mathbf{v_{0}}$, for any fixed number of iterations. This corresponds to a converse to Theorem 6.

Since state evolution is asymptotically exact, we can prove sharp phase transition results with explicit characterization of their locations.

We assume that the additional information takes the form of a noisy observation ${\mathbf{y}}=\gamma\,\mathbf{v_{0}}+\mathbf{z}$, where $\mathbf{z}\sim{\sf N}(0,{\rm I}_{n}/n)$. Our next results summarizes the state evolution analysis. Its proof is deferred to a journal publication.

where $\mathbf{v}^{t}_{\parallel}$ is proportional to $\mathbf{v_{0}}$, and $\mathbf{v}^{t}_{\perp}$ is perpendicular. Then $\mathbf{v}^{t}_{\perp}$ is uniformly random, conditional on its norm. Further, almost surely

where $\tau_{t}$ is given recursively by letting $\tau_{0}=\gamma$ and, for $t\geq 0$ (we refer to this as to state evolution):

Note that state evolution coincides with the equation that we derived for the first iteration of power iteration, cf. Eq. (48) (apart from the different scaling). It is important to notice that for subsequent iterations $t\geq 1$, state evolution (53) does not correctly describe naive power iteration. The reason is that $\mathbf{v}^{t}$ depends on $\mathbf{X}$, and hence the same argument does not apply. The AMP iteration differ from naive power iteration because of the ‘memory term’, $-{\sf b}_{t}\,f(\mathbf{v}^{t-1})$. As shown in [BM11], this memory term approximately cancels dependencies. As a consequence, the resulting algorithm obeys state evolution.

The following result characterizes the minimum required additional information $\gamma$ to allow AMP to escape from those undesired local optima. We will say that $\{\mathbf{v}^{t}\}_{t}$ converges almost surely to a desired local optimum if, almost surely,

Consider the Tensor PCA problem with $k\geq 3$ and

Then AMP converges almost surely to a desired local optimum if and only if $\gamma>\sqrt{1/\epsilon_{k}(\beta)-1}$ where $\epsilon_{k}(\beta)$ is the largest solution of $(1-\epsilon)^{(k-2)}\epsilon=\beta^{-2}$,

In the special case $k=3$, and $\beta>2$, assuming $\gamma>\beta(1/2-\sqrt{1/4-1/\beta^{2}})$, AMP tends to a desired local optimum. Numerically $\beta>2.69$ is enough for AMP to achieve $\langle\mathbf{v_{0}},\mathbf{\widehat{v}}\rangle\geq 0.9$ if $\gamma>0.45$.

As a final remark, we note that the methods of [MR14] can be used to show that, under the assumptions of Theorem 7, for $\beta>\beta_{k}$ a sufficiently large constant, AMP asymptotically solves the optimization problem Tensor PCA. Formally, we have, almost surely,

Numerical experiments

Let us emphasize two practical suggestions that arise from our work:

Tensor unfolding is superior to tensor power iteration under our spiked model. For instance, for $k=3$, we expect tensor power iteration to require $\beta\gtrsim n^{1/4}$ and unfolding to require $\beta\gtrsim n^{1/2}$.

For smaller values of $\beta$, iterative methods (tensor power iteration or approximate message passing) only produce a good estimate if the initialization has a scalar product with the ground truth $\mathbf{v_{0}}$ that is bounded away from zero.

As a consequence of the above, side information about the unknown vector $\mathbf{v_{0}}$ can greatly improve performances.

A special case, we will study the behavior of warm start algorithms that first perform a singular value decomposition of ${\sf Mat}(\mathbf{X})$, and then apply an iterative method (tensor power iteration or approximate message passing).

In this section we will illustrate these suggestions through numerical simulations.

Section 6.1 describes a refinement of tensor unfolding that provides a tighter relaxation. Section 6.2 compares different algorithms. Finally, Section 6.3 provides additional illustration of how side information can dramatically simplify the estimation problem.

This optimization problem is NP hard, since it includes copositive programming as a special case. However [DMR14] provides rigorous and empirical evidence that problems of this type can be solved efficiently by a projected power iteration, under statistical model dor $\mathbf{X}$.

2 Comparison of different algorithms

In Fig. 1 we compare different algorithms on data generated following Spiked Tensor Model with $k=3$, and $n\in\{25,50,100,200,400,800\}$ and for a range of values of $\beta\in$. The plots represent measured values of the absolute correlation $|\langle\mathbf{\widehat{v}},\mathbf{v_{0}}\rangle|$ versus $\beta$, averaged over $50$ samples (except for $n=800$, where we used $8$ samples).

The main findings are consistent with the theory developed above:

Tensor power iteration (with random initialization) performs poorly with respect to other approaches that use some form of tensor unfolding. The gap widens as the dimension $n$ increases.

PSD-constrained principal component analysis (described in the last section) is slightly superior to plain unfolding.

All algorithms based on initial unfolding have essentially the same threshold. Above that threshold, those that process the singular component (either by recursive unfolding or by tensor power iteration) have superior performances over simpler one-step algorithms.

In addition, we noted that the two iterative algorithms (Power Iteration and AMP) show very close behaviors in our experiments.

In Figure 2 we compare the scaling with $n$ of the threshold signal-to-noise ratio for different type of algorithms. Our heuristic arguments suggest that tensor power iteration with random initialization will work for $\beta\gtrsim n^{1/2}$, while unfolding only requires $\beta\gtrsim n^{1/4}$ (our theorems guarantee this for, respectively, $\beta\gtrsim n$ and $\beta\gtrsim n^{1/2}$). We plot the average correlation $|\langle\mathbf{\widehat{v}},\mathbf{v_{0}}\rangle|$ versus (respectively) $\beta/n^{1/2}$ and $\beta/n^{1/4}$. The curve superposition confirms that our prediction captures the correct behavior already for $n$ of the order of $50$.

3 The value of side information

The analysis in previous sections suggest to use the leading eigenvector of $\mathbf{M}$ as the initial point of AMP algorithm for tensor PCA on $\mathbf{X}$. We performed the experiments on $100$ randomly generated instances with $n=50,200,500$ and report in Figure 3 the mean values of $|\langle\mathbf{v_{0}},\mathbf{\widehat{v}}(\mathbf{X})\rangle|$ with confidence intervals.

Random matrix theory predicts $\lim_{n\to\infty}\langle\mathbf{\widehat{v}}_{1}(M),\mathbf{v_{0}}\rangle=\sqrt{1-\lambda^{-2}}$ [FP09]. Thus we can set $\gamma=\sqrt{1-\lambda^{-2}}$ and apply the theory of the previous section. In particular, Proposition 5.1 implies

and $\lim_{n\to\infty}\langle\mathbf{\widehat{v}}(\mathbf{X}),\mathbf{v_{0}}\rangle=0$ otherwise Simultaneous PCA appears vastly superior to simple PCA. Our theory captures this difference quantitatively already for $n=500$.

Acknowledgements

This work was partially supported by the NSF grant CCF-1319979 and the grants AFOSR/DARPA FA9550-12-1-0411 and FA9550-13-1-0036.

Appendix A Information theoretic bound: Proof of Theorem 1

where for indices $i_{1}<i_{2}<\cdots<i_{k}$, we have ${\sf{U}}(\mathbf{X})_{a(i_{1},\cdots,i_{k})}=\mathbf{X}_{i_{1},\cdots,i_{k}}$ with $a(i_{1},\cdots,i_{k})=1+\sum_{j=1}^{k}n^{j-1}(i_{j}-1)$. Let $D(\cdot\|\cdot)$ denote the Kullback-Leiber divergence where $P_{\mathbf{w}}$ is the law of ${\sf{U}}(\mathbf{X})$ conditional on $\mathbf{v}_{0}=\mathbf{w}$.

We are now in position to prove Theorem 1. Let $\mathcal{V}$ denote the class of estimators $\mathbf{\widehat{v}}$ with unit norm:

(Here ${\rm Vol}_{n-1}(\,\cdot\,)$ denotes the $(n-1)$-dimensional volume, and $B(\mathbf{x},\varepsilon)$ the ball of radius $\varepsilon$ centered at $\mathbf{x}$.)

Let $\mathcal{N}$ denote an $\varepsilon$-packing with cardinality $|\mathcal{N}|\geq N_{n}(\varepsilon)$. Let $\mathbf{v_{0}}$ be uniformly distributed in the set $\mathcal{N}$. For an estimator $\mathbf{\widehat{v}}\in\mathcal{V}$, we define $G(\mathbf{\widehat{v}}(\mathbf{X}))=\arg\min_{\mathbf{w}\in\mathcal{N}}\|\mathbf{\widehat{v}}(\mathbf{X})-\mathbf{w}\|_{2}$. Consider the error event $\{G(\mathbf{\widehat{v}}(\mathbf{X}))\neq\mathbf{v_{0}}\}$. By definition of $G(\mathbf{\widehat{v}}(\mathbf{X}))$, the event $G(\mathbf{\widehat{v}}(\mathbf{X}))\neq\mathbf{v_{0}}$ implies $(\|\mathbf{\widehat{v}}(\mathbf{X})-\mathbf{v}_{0}\|_{2}\wedge\|\mathbf{\widehat{v}}(\mathbf{X})+\mathbf{v}_{0}\|_{2})\geq\varepsilon/2$. By Markov inequality we have:

By Fano’s inequality [CT12] we have that:

where $\Delta=\max_{\mathbf{w}\neq\mathbf{w}^{\prime}\in\mathcal{N}}D(P_{\mathbf{w}}\|P_{\mathbf{w}^{\prime}})$, and in the second inequality we used [HV94]

Using Eq. (59) and Lemma A.1, in Eq. (61), we get

Appendix B Maximum likelihood: Proof Theorem 2

While the function $H_{\mathbf{Z}}(\,\cdot\,)$ is obviously non-convex, it turns out that –for random data $\mathbf{Z}$– it is dramatically so. Namely, it has an exponential number of local maximum, whose value is –typically– only a fraction of the value of the global maximum.

where, for $x\geq\eta_{k}\equiv 2\sqrt{k-1}$

Further, for $x<\eta_{k}$, $g_{k}(x)=g_{k}(\eta_{k})$.

The function $g_{k}(x)$ is monotone decreasing for $x\geq\eta_{k}$, and non-negative if and only if $x\in[\eta_{k},\mu_{k}]$ for some $\mu_{k}>0$ (strictly positive for $x\in[\eta_{k},\mu_{k})$). In Figure 4, we plot $g_{k}(x)$ for $k\in\{3,4,5\}$. Informally, this means that the function $H_{\mathbf{Z}}(\mathbf{v})$ has exponentially many local maxima with value $H_{\mathbf{Z}}(\mathbf{v})\approx x$ for any $x\in[\eta_{k},\mu_{k})$. To leading exponential order, the number of such maxima is given by $\exp\{n\,g_{k}(x)\}$.

The value $\mu_{k}$ can be determined as the unique solution to the equation $g(x)=0$. It is immediately to do this numerically, obtaining the values in Section 2.

The last result implies that the global maximum of $H_{\mathbf{Z}}(\mathbf{v})$ is (asymptotically) at least $\mu_{k}$. Indeed the global maximum is necessarily a local maximum as well. The next result implies that indeed the global maximum converges to $\mu_{k}$.

Let $\mu_{k}$ denote the unique non-negative root of the equation $g_{k}(x)=0$, for $x\geq\eta_{k}\equiv 2\sqrt{k-1}$. Then

In order to derive the large-$k$ asymptotics of $\mu_{k}$, we rewrite Eq. (69) in terms of the variable $y\equiv k^{2}z^{2}/2$. We get $g_{k}(x)=f_{k}(y(x))/2$, where

Further $y\in(0,k/(k-1)]$. The claimed asymptotics follows by showing that the only solution of $f_{k}(y)=0$ in this interval obeys $y_{k}=(\log k)^{-1}(1+o_{k}(1))$. This in turns can be showed by using the bounds

and showing that the solution of $y^{-1}+\log(y)=\log(a)$ for large $a$ is $y^{-1}=a+\Theta(\log(a))$.

Finally, the norm $\|\mathbf{Z}\|_{op}$ concentrates tightly around its expectation.

is a Lipschitz function with Lipschitz modulus $\sqrt{k/n}$ (with respect to Euclidean norm) of the Gaussian vector (tensor) $\mathbf{G}$. Hence $\|\mathbf{Z}\|_{op}$ is Lipchitz continuous with the same modulus. The claim follows from Gaussian isoperimetry [Led01]. ∎

Note that to make the connection with the notations used in [ABAC13], one has to use the proper scaling $H_{n,k}(\mathbf{v})=\frac{n}{\sqrt{k}}\,L_{\mathbf{Z}}(\mathbf{v}/\sqrt{n})$ ($H_{n,k}(\mathbf{v})$is the objective function considered in [ABAC13]).

The upper bound on the tensor operator norm obtained from Sudakov-Fernique inequality is not tight. In fact taking $\beta=0$ in Lemma 2.2 gives the loose upper bound $\|\mathbf{Z}\|_{op}\leq k$. Except in the case of random matrices ($k=2$), this is loose roughly by a factor $\sqrt{k}$.

B.2 Proof of Theorem 2

By optimality of $\mathbf{\widehat{v}}$, we have

Using $(1-\alpha)^{1/k}\geq(1-\alpha)$ which holds for $\alpha\in$, and rescaling $s$, we get $|\langle\mathbf{v_{0}},\hat{\mathbf{v}}\rangle|\geq 1-(\mu_{k}+s)/\beta$ with probability at least $1-2e^{-ns^{2}/(16k)}$ for all $n$ large enough.

B.3 Proof of Lemma 2.2

Since $x\mapsto\sqrt{x}$ is uniformly continuous on bounded intervals $[0,M]$, it is sufficient to prove

for all $n\geq n_{0}$, and an eventually different sequence $\delta_{n}$. By triangular inequality and using $\|\mathbf{v_{0}}\|_{2}=1$,

Next we have $\lim_{n\to\infty}\|\mathbf{g}\|^{2}=1$ almost surely by the strong law of large numbers, and $\langle\mathbf{v_{0}},\mathbf{g}\rangle\sim{\sf N}(0,1/n)$ whence $\lim_{n\to\infty}\langle\mathbf{v_{0}},\mathbf{g}\rangle=0$ by Borel-Cantelli. ∎

The function $\mathbf{X}\mapsto M_{\mathbf{X}}(\kappa)$ is a Lipschitz continuous function with Lipschitz constant $\sqrt{k/n}$ of the standard Gaussian tensor $\mathbf{G}$ (namely $|M_{\mathbf{X}}(\kappa)-M_{\mathbf{X}^{\prime}}(\kappa)|\leq(k/n)^{1/2}\|\mathbf{G}-\mathbf{G}^{\prime}\|_{F}$). Hence, by Gaussian isoperimetry, we have

Further we claim that $\kappa\mapsto M_{\mathbf{X}}(\kappa)$ is uniformly continuous for $\kappa\in$. In order to prove this, let

where $\langle\mathbf{v}^{\perp},\mathbf{v_{0}}\rangle=0$. We have, for $\kappa_{1},\kappa_{2}\in$, and by letting $\mathbf{v}^{\perp}$ and $\mathbf{w}^{\perp}$ denote the perpendicular components of $\mathbf{v}(\kappa_{1})$ and $\mathbf{v}(\kappa_{2})$, we have for some constant $c>0$

where Eq. (88) was obtained by exploiting the symmetry of the tensor $\mathbf{X}$ and Eq. (89) was derived using the norm of the vector $\left\{(\kappa_{1}-\kappa_{2})\mathbf{v_{0}}+(\sqrt{1-\kappa_{1}^{2}}-\sqrt{1-\kappa_{2}^{2}})\mathbf{w}^{\perp}\right\}$. Using Eq. (86) over a grid $\kappa\in\{0,1/n,2/n,\dots,1-1/n,1\}$, and the factThis follows from Lemma 2.1 and triangular inequality, or from a standard $\varepsilon$-net argument. that $\|\mathbf{X}\|_{op}\leq C$ for some constant $C>0$ with probability $1-e^{-\Theta(n)}$, we have for all $t>0$ and some constant $c>0$

In particular, by Borel-Cantelli we have, almost surely,

In order to upper bound $\overline{M}(\kappa)$, we apply Sudakov-Fernique inequality for non-centered Gaussian processes [Vit00, Theorem 1] to the two processes $\{{\cal X}_{\mathbf{v}}\}$, $\{{\cal Y}_{\mathbf{v}}\}$ indexed by $\mathbf{v}\in\mathcal{W}_{\kappa}$ defined as follows:

where $\tau=\kappa/\sqrt{1-\kappa^{2}}$ and $\delta_{n}$ satisfies $\lim_{n\to\infty}\delta_{n}=0$ uniformly over $\kappa\in$, by Lemma B.4. We finally conclude that

Concentration around the expectation follows by Gaussian isoperimetry as in the proof of Lemma 2.1. ∎

Appendix C Power Iteration: Proof of Theorem 6

Let $\tau_{t}\equiv\langle\mathbf{v_{0}},\mathbf{v}^{t}\rangle$ and $\xi\equiv\|\mathbf{Z}\|_{op}/\beta$. Let $\tau_{\rm min}$, $\tau_{*}\in$ be the two solutions of

We will show below that our assumptions imply $\tau_{0}>\tau_{\rm min}$. Further $\tau\geq\tau_{\rm min}$ implies $\tau^{k-1}-\xi\geq 0$.

We will prove the first inequality $\tau_{t}\geq\tau_{\rm min}$ by induction. It is true at $t=0$ by assumption. Assume it is true at $t$. Then $\tau_{t+1}\geq 0$ using Eq. (101).

Hence we can divide the two inequalities above obtaining $\tau_{t+1}\geq(\tau_{t}^{k-1}-\xi)/(\tau_{t}^{k-1}+\xi)$ which implies

To conclude the proof of Eq. (43), we notice that, for $\xi\leq 1/(2e(k-1))$

where we recall that $\tau_{\rm min}$, $\tau_{*}$ are the two solutions of $g_{k}(x)\equiv x^{k-1}(1-x)=\xi$ in the interval $$. For the first inequality, note that, in the interval$[e^{-1/(k-1)},1]$,$g_{k}(x)$is decreasing with$g_{k}(x)\geq e^{-1}(1-x)$. This implies

i.e. $\tau_{*}\geq 1-e\,\xi$ as long as $1-e~{}\xi\geq e^{-1/(k-1)}$, which is implied by $\xi\leq 1/(2e(k-1))$.

For the second inequality, note that, in the interval $[0,1-(k-1)^{-1}]$, we have $g_{k}(x)$ increasing with $g_{k}(x)\geq x^{k-1}/(k-1)$. This implies

as long as $[(k-1)\xi]^{1/k-1}\leq 1-(k-1)^{-1}$, which follows, again, by our assumptions.

Finally, conditions (44), (45) follow directly by applying Lemma 2.1.

Appendix D Approximate Message Passing: Proof of Theorem 7

Let us recall the state evolution recursion (53)

The fixed point equation $\tau^{2}=f(\tau^{2};\beta)$ has two strictly positive solutions $\tau_{1}^{2}(\beta)<\tau_{2}^{2}(\beta)$.

The smallest fixed point is given by $\tau_{1}(\beta)=\sqrt{1/\epsilon_{k}(\beta)-1}$ as in the statement.

The largest fixed point satisfies $\tau_{2}(\beta)>1-(2/\beta^{2})$.

The behavior of the function $f(\tau^{2};\beta)$ is illustrated in Fig. 5.

Now, the function $x\mapsto h_{k}(x)$ is continuously differentiable and strictly positive in the interval $(0,1)$, with $h_{k}(0)=h_{k}(1)=0$. Further, simple calculus shows it has a unique stationary point (a maximum) at $x_{*}=(k-2)/(k-1)$ with $h_{k}(x_{*})=1/\omega_{k}^{2}$. This implies that, for $\beta>\omega_{k}$, Eq. (111) has two fixed points $0<x_{2}(\beta)<x_{*}<x_{1}(\beta)<1$ thus implying points 1 and 2 above (the latter immediately follows from inverting the re-parametrization).

where the second inequality follows since $x_{2}>x_{*}$. By state evolution (Proposition 5.1), together with the fact that $\tau_{t}\to\tau_{2}$, we have

References