Guaranteed Non-Orthogonal Tensor Decomposition via Alternating Rank-$1$ Updates

Tensor decomposition, alternating minimization, overcomplete representation, latent variable models.

Introduction

Tensor decompositions have been recently popular for unsupervised learning of a wide range of latent variable models such as independent component analysis (De Lathauwer et al., 2007), topic models, Gaussian mixtures, hidden Markov models (Anandkumar et al., 2014a), network community models (Anandkumar et al., 2013a), and so on. The decomposition of a certain low order multivariate moment tensor (typically up to fourth order) in these models is guaranteed to provide a consistent estimate of the model parameters. Moreover, the sample and computational requirements are only a low order polynomial in the rank of the tensor (Anandkumar et al., 2014a; Song et al., 2013). In practice, the tensor decomposition techniques have been shown to be effective in a number of applications such as blind source separation (Comon, 2002), computer vision (Vasilescu and Terzopoulos, 2003), contrastive topic modeling (Zou et al., 2013), and community detection (Huang et al., 2013). In many cases, the tensor approach is shown to be orders of magnitude faster than existing techniques such as the stochastic variational approach.

The state of art for guaranteed tensor decomposition involves two steps: converting the input tensor to an orthogonal symmetric form, and then solving the orthogonal decomposition through tensor eigen decomposition (Comon, 1994; Kolda and Mayo, 2011; Zhang and Golub, 2001; Anandkumar et al., 2014a). The first step of converting the input tensor to an orthogonal symmetric form is known as whitening. For the second step, the tensor eigen pairs can be found through a simple tensor power iteration procedure.

While having efficient guarantees, the above procedure suffers from a number of theoretical and practical limitations. For instance, in practice, the learning performance is especially sensitive to whitening (Le et al., 2011). Moreover, whitening is computationally the most expensive step in deployments (Huang et al., 2013), and it can suffer from numerical instability in high-dimensions due to ill-conditioning. Lastly, the above approach is unable to learn overcomplete representations (this is the case when number of features/components is much larger than the dimension) due to the orthogonality constraint, which is especially limiting, given the recent popularity of overcomplete feature learning in many domains (Bengio et al., 2012; Lewicki and Sejnowski, 2000).

The current practice for tensor decomposition is the alternating least squares (ALS) procedure, which has been described as the “workhorse” of tensor decomposition (Kolda and Bader, 2009). This involves solving the least squares problem on a mode of the tensor, while keeping the other modes fixed, and alternating between the tensor modes. The method is extremely fast since it involves calculating linear updates, but is not guaranteed to converge to the global optimum in general (Kolda and Bader, 2009).

In this paper, we provide local and global convergence guarantees for a modified alternating method, for which the main step is making rank-$1$ updates along different modes of the tensor. This update is basically a rank-1 ALS update. This method is extremely fast to deploy, trivially parallelizable, and does not suffer from ill-conditioning issues faced by both ALS (Kolda and Bader, 2009) and whitening approaches (Le et al., 2011). Our analysis assumes the presence of incoherent tensor components, which can be viewed as a soft-orthogonality constraint. Incoherent representations have been extensively considered in literature in a number of contexts, e.g., compressed sensing (Donoho, 2006) and sparse coding (Arora et al., 2013; Agarwal et al., 2013). Incoherent representations provide flexible modeling, can handle overcomplete signals, and are robust to noise (Lewicki and Sejnowski, 2000). Moreover, when the latent variable model parameters are generic or when we have randomly constructed (multiview) features (McWilliams et al., 2013), the moment tensors have incoherent components, as assumed here. In this work, we establish that incoherence leads to efficient guarantees for tensor decomposition. The guarantees also include a tight perturbation analysis. In a subsequent work (Anandkumar et al., 2014b), we apply the tensor decomposition guarantees of this paper to various learning settings, and derive sample complexity bounds through novel covering arguments.

In this paper, we propose and analyze an algorithm for non-orthogonal CP (Candecomp/Parafac) tensor decomposition; see Figure 1 for the details of the algorithm. The main step of the algorithm is a simple alternating rank-$1$ update which is the alternating version of the tensor power iteration adapted for asymmetric tensors. In each iteration, one of the tensor modes is updated by projecting the other modes along their estimated directions, and the process is alternated between all the modes of the tensor; see (5) for this update.

In the undercomplete or mildly overcomplete settings $(k=O(d))$, a simple initialization procedure (see Procedure 2) based on rank-$1$ SVD of random tensor slices is provided. This initialization procedure lands the estimate in the basin of attraction for the alternating update procedure in polynomial number of trials (in the tensor rank $k$). This leads to global convergence guarantees for tensor decomposition.

Greedy or rank-$1$ updates are perhaps the most natural procedure for CP tensor decomposition. For orthogonal tensors, they lead to guaranteed recovery (Zhang and Golub, 2001). However, when the tensor is non-orthogonal, greedy procedure is not optimal in general (Kolda, 2001). Finding tensor decomposition in general is NP-hard (Hillar and Lim, 2009). We circumvent this obstacle by limiting ourselves to tensors with incoherent components. We exploit incoherence to prove error contraction under each step of the alternating update procedure with an approximation error, which is decaying, when $k=o(d^{1.5})$. To this end, we require tools from random matrix theory, bounds on $2\to p$ norm for random matrices (Guédon and Rudelson, 2007; Adamczak et al., 2011) for some $p<3$, and matrix perturbation results to provide tight bounds on error contraction.

2 Related work

CP tensor decomposition (Carroll and Chang, 1970), also known as PARAFAC decomposition (Harshman, 1970; Harshman and Lundy, 1994) is a classical definition for tensor decomposition with many applications. The most commonly used algorithm for CP decomposition is Alternating Least Squares (ALS) (Comon et al., 2009), which has no convergence guarantees in general. Kolda (2001) and Zhang and Golub (2001) analyze the greedy or the rank-1 updates in the orthogonal setting. In the noisy setting, Anandkumar et al. (2014a) analyze deflation procedure for orthogonal decomposition, and Song et al. (2013) extend the analysis to the nonparametric setting. For the non-orthogonal tensors, a common strategy is to first apply a procedure called whitening to reduce it to the orthogonal case. But as discussed earlier, the whitening procedure can lead to poor performance and bad sample complexity. Moreover, it requires the tensor factors to have full column rank, which rules out overcomplete tensors.

Learning overcomplete tensors is challenging, and they may not even be identifiable in general. Kruskal (1976, 1977) provided an identifiability result based on the Kruskal rank of the factor matrices of the tensor. However, this result is limiting since it requires $k=O(d)$, where $k$ is the tensor rank and $d$ is the dimension. The FOOBI procedure by De Lathauwer et al. (2007) overcomes this limitation by assuming generic factors, and shows that a polynomial-time procedure can recover the tensor components when $k=O(d^{2})$, and the tensor is fourth order. However, the procedure does not work for third-order overcomplete tensors, and has no polynomial sample complexity bounds. Simple procedures can recover overcomplete tensors for higher order tensors (five or higher). For instance, for the fifth order tensor, when $k=O(d^{2})$, we can utilize random slices along a mode of the tensor, and perform simultaneous diagonalization on the matricized versions. Note that this procedure cannot handle the same level of overcompleteness as FOOBI, since an additional dimension is required for obtaining two (or more) fourth order tensor slices. The simultaneous diagonalization procedure entails careful perturbation analysis, carried out by (Goyal et al., 2013; Bhaskara et al., 2013). In addition, Goyal et al. (2013) provide stronger results for independent components analysis (ICA), where the tensor slices can be obtained in the Fourier domain.

There are other recent works which can learn overcomplete models, but under different settings than the ones considered in this paper. For instance, Arora et al. (2013); Agarwal et al. (2013) provide guarantees for the sparse coding problem. Anandkumar et al. (2013b) learn overcomplete sparse topic models, and provide guarantees for Tucker tensor decomposition under sparsity constraints. Specifically, the model is identifiable using $(2n)^{{\mbox{\tiny th}}}$ order moments when the latent dimension $k=O(d^{n})$ and the sparsity level of the factor matrix is $O(d^{1/n})$, where $d$ is the observed dimension. The Tucker decomposition is different from the CP decomposition considered here (it has weaker assumptions and guarantees), and the techniques in (Anandkumar et al., 2013b) differ significantly from the ones considered here.

The algorithm employed here falls under the general framework of alternating minimization. There are many recent works which provide guarantees on local/global convergence for alternating minimization, e.g., for matrix completion (Jain et al., 2013; Hardt, 2013), phase retrieval (Netrapalli et al., 2013) and sparse coding (Agarwal et al., 2013). However, the techniques in this paper are significantly different, since they involve tensors, while the previous works only required matrix analysis.

3 Notations and tensor preliminaries

Let $[n]$ denote the set $\{1,2,\dotsc,n\}$.

Tensor Decomposition Algorithm

In this section, we introduce the alternating tensor decomposition algorithm, and the guarantees are provided in Section 3. The goal of tensor decomposition algorithm is to recover the rank-1 components of tensor; see (4) for the notion of tensor rank. Figure 1 depicts the overview of our tensor decomposition method where the corresponding algorithms and procedures are also specified. Our algorithm includes two main steps as 1) alternating tensor power iteration, and 2) coordinate descent iteration for removing the residual error. The former one is performed in Algorithm 1 (see equation (5), and the latter one is done in Algorithm 4 (see equation (9)). We now describe these steps of the algorithm in more details as well as providing the auxiliary procedures required to complete the algorithm.

The main step of the algorithm is tensor power iteration which basically performs alternating asymmetric power updates This is exactly the generalization of asymmetric matrix power update to $3$rd order tensors. on different modes of the tensor as

Consider the problem of best rank-$1$ approximation of tensor $T$ as

where ${\cal S}^{d-1}$ denotes the unit $d$-dimensional sphere. This optimization program is non-convex, and has multiple local optima. It can be shown that the updates in (5) are the alternating optimization for this program where in each update, optimization over one vector is performed while the other two vectors are assumed fixed. This alternating minimization approach does not converge to the true components of tensor $T$ in general, and in this paper we provide sufficient conditions for the convergence guarantees.

We now provide an intuitive argument on the functionality of power updates in (5). Consider a rank-$k$ tensor $T$ as in (4), and suppose we start at the correct vectors $\hat{a}=a_{j}$ and $\hat{b}=b_{j}$, for some $j\in[k]$. Then for the numerator of update formula (5), we have

where the first term is along $c_{j}$ and the second term is an error term due to non-orthogonality. For orthogonal decomposition, the second term is zero, and the true vectors $a_{j},b_{j}$ and $c_{j}$ are stationary points for the power update procedure. However, since we consider non-orthogonal tensors, this procedure cannot recover the decomposition exactly leading to a residual error after running this step. Under incoherence conditions which encourages soft-orthogonality constraints See Assumption (A2) in Appendix A for precise description. (and some other conditions), we show that the residual error is small (see Lemma 1 where the guarantees for the tensor power iteration step is provided), and thus, with the additional step we propose in Section 2.2, we can also remove this residual error.

We discussed that the tensor power updates in (5) are the alternating iterations for the problem of rank-1 approximation of the tensor; see (6). This is a non-convex problem and has many local optima. Thus, the power update requires careful initialization to ensure convergence to the true rank-1 tensor components.

Notice that the algorithm is run for $L$ different initialization vectors for which we do not know the good ones in prior. In order to identify which initializations are successful at the end, we also need a clustering step proposed in Procedure 3 to obtain the final estimates of the vectors. The detailed analysis of clustering procedure is provided in Appendix D.

2 Coordinate descent iteration in Algorithm 4

We discussed in the previous section that the tensor power iteration recovers the tensor rank-1 components up to some residual error. We now propose Algorithm 4 to remove this additional residual error. This algorithm mainly runs a coordinate descent iteration as

The analysis of this algorithm requires that the estimate matrices $\hat{A},\hat{B},\hat{C}$ satisfy some bound on the spectral norm and some column-wise error bounds; see Definition 2 in Appendix B.2 for the details. The optimization program in (10) (which is only run in the first iteration) and projection Procedure 5 ensure that these conditions are satisfied.

3 Discussions

We now provide some further discussions and comparisons about the algorithm.

In many applications, the input tensor $T$ is not available in advance, and it is computed from samples. It is discussed in (Anandkumar et al., 2014b) that the tensor is not needed to be computed and stored explicitly, where the multilinear tensor updates (5) and (9) in the algorithm can be efficiently computed through multilinear operations on the samples directly.

Algorithm 1 is similar to the symmetric tensor power method analyzed by Anandkumar et al. (2014a) with the following main differences, viz.,

Symmetric and non-symmetric tensors: Our algorithm can be applied to both symmetric and non-symmetric tensors, while tensor power method in Anandkumar et al. (2014a) is only for symmetric tensors.

Linearity: The updates in Algorithm 1 are linear in each variable, while the symmetric tensor power update is a quadratic operator given a third order tensor.

Guarantees: In Anandkumar et al. (2014a), guarantees for the symmetric tensor power update under orthogonality are obtained, while here we consider non-orthogonal tensors under the alternating updates.

The updates in Algorithm 1 can be viewed as a rank-$1$ form of the standard alternating least squares (ALS) procedure. This is because the unnormalized update for $c$ in (5) can be rewritten as

Analysis

We require natural deterministic conditions on the tensor components to argue the convergence guarantees; see Appendix A for the details. We show that all of these conditions are satisfied if the true rank-1 components of the tensor are uniformly i.i.d. drawn from the unit $d$-dimensional sphere ${\cal S}^{d-1}$. Thus, for simplicity we assume this random assumption in the main part, and state the deterministic assumptions in Appendix A. Notice that it is also reasonable to assume these deterministic assumptions hold for some non-random matrices. Among the deterministic assumptions, the most important one is the incoherence condition which imposes a soft-orthogonality constraint between different rank-1 components of the tensor.

The convergence guarantees are provided in terms of distance between the estimated and the true vectors, defined below.

Note that distance function $\operatorname{dist}(u,v)$ is invariant w.r.t. norm of input vectors $u$ and $v$. Distance also provides an upper bound on the error between unit vectors $u$ and $v$ as (see Lemma A.1 of Agarwal et al. (2013))

Incorporating distance notion resolves the sign ambiguity issue in recovering the components: note that a third order tensor is unchanged if the sign of a vector along one of the modes is fixed and the signs of the corresponding vectors in the other two modes are flipped.

In the local convergence guarantee, we analyze the convergence properties of the algorithm assuming we have good initialization vectors for the non-convex tensor decomposition algorithm.

Rank-$k$ true tensor with random components: Let

where $a_{i},b_{i},c_{i},i\in[k],$ are uniformly i.i.d. drawn from the unit $d$-dimensional sphere ${\cal S}^{d-1}$. We state the deterministic assumptions in Appendix A, and show that random matrices satisfy these assumptions.

Rank condition: $k=o\left(d^{1.5}\right).$

Perturbation tensor $\Psi$ satisfies the bound

Weight ratio: The maximum ratio of weights $\gamma:=\frac{w_{\max}}{w_{\min}}$ satisfies the bound

Initialization: Assume we have good initialization vectors $\hat{a}^{(0)}_{j},\hat{b}^{(0)}_{j},j\in[k]$ satisfying

where $\gamma:=\frac{w_{\max}}{w_{\min}}$. In addition, given $\hat{a}^{(0)}_{j}$ and $\hat{b}^{(0)}_{j}$, suppose $\hat{c}^{(0)}_{j}$ is also calculated by the update formula in (5).

Same error bounds hold for other factor matrices $B:=[b_{1}\dotsb b_{k}]$ and $C:=[c_{1}\dotsb c_{k}]$.

Thus, we can efficiently decompose the tensor in the highly overcomplete regime $k\leq o\left(d^{1.5}\right)$ under incoherent factors and some other assumptions mentioned above. The deterministic version of assumptions are stated in Appendix A. We show that these assumptions are true for random components which is assumed here for simplicity. If $k$ is significantly smaller than $d^{1.5}$ ($k\ll d^{1.25}$), then many of the assumptions can be derived from incoherence. See Appendix A for the details.

The above local convergence result can be also interpreted as a local identifiability result for tensor decomposition under incoherent factors.

The $\sqrt{k}$ factor in the above theorem error bound is from the fact that the final recovery guarantee is on the Frobenius norm of the whole factor matrix $A$. In the following, we provide stronger column-wise guarantees (where there is no $\sqrt{k}$ factor) with the expense of having an additional residual error term. Recall that our algorithm includes two main update steps including tensor power iteration in (5) and residual error removal in (9). The guarantee for the first step — tensor power iteration — is provided in the following lemma.

Consider the same settings as in Theorem 1. Then, the outputs of tensor power iteration steps (output of Algorithm 1) satisfy w.h.p.

Same error bounds hold for other factor matrices $B$ and $C$.

The result in the above lemma is actually stated in the non-asymptotic form, where the details of constants are explicitly provided in Appendix A.

The above local convergence result also holds for recovering the components of a rank-$k$ symmetric tensor. Consider symmetric tensor $T$ with CP decomposition $T=\sum_{i\in[k]}w_{i}a_{i}\otimes a_{i}\otimes a_{i}$. The proposed algorithm can be also applied to recover the components $a_{i},i\in[k],$ where the main updates are changed to adapt to the symmetric tensor. The tensor power iteration is changed to

and the coordinate descent update is changed to the form stated in (28). Then, the same local convergence result as in Theorem 1 holds for this algorithm. The proof is very similar to the proof of Theorem 1 with some slight modifications considering the symmetric structure.

We also provide the generalization of the tensor decomposition guarantees to higher order tensors. We state and prove the result for the tensor power iteration part in details, while the generalization of coordinate descent part (for removing the residual error) to higher order tensors, can be argued by the same techniques we introduce in this paper

and similarly the other updates are changed. Then, we have the following generalization of Lemma 1 to higher order tensors.

Consider the same conditions and settings as in Lemma 1, unless tensor $T$ is $p$-th order with CP decomposition in (15) where $p\geq 3$ is a constant. In addition, the bounds on $\gamma:=\frac{w_{\max}}{w_{\min}}$ and $k$ are modified as

Then, the outputs of tensor power iteration steps (output of Algorithm 1) satisfy w.h.p.

2 Global convergence guarantee when k=O​(d)𝑘𝑂𝑑k=O(d)

Theorem 1 provides local convergence guarantee given good initialization vectors. In this section, we exploit SVD-based initialization method in Procedure 2 to provide good initialization vectors when $k=O(d)$. This method proposes the top singular vectors of random slices of the moment tensor as the initialization. Combining the theoretical guarantees of this initialization method (provided in Appendix C) with the local convergence guarantee in Theorem 1, we provide the following global convergence result.

The initialization in each run of Algorithm 1 is performed by SVD-based technique proposed in Procedure 2, with the number of initializations as

Rank-$k$ decomposition and perturbation conditions as Note that the perturbation condition is stricter than the corresponding condition in the local convergence guarantee (Theorem 1).

where $a_{i},b_{i},c_{i},i\in[k],$ are uniformly i.i.d. drawn from the unit $d$-dimensional sphere ${\cal S}^{d-1}$, and $\alpha_{0}>1$ is a constant.

Rank condition: $k=O(d)$, i.e., $k\leq\beta d$ for arbitrary constant $\beta>1$.

Consider noisy rank-$k$ tensor $\hat{T}=T+\Psi$ as the input to the tensor decomposition algorithm, and assume the conditions and settings mentioned above hold. Then, the same guarantees as in Theorem 1 hold.

Thus, we can efficiently recover the tensor decomposition, when the tensor is undercomplete or mildly overcomplete (i.e., $k\leq\beta d$ for arbitrary constant $\beta>1$), by initializing the algorithm with a simple SVD-based technique. The number of initialization trials $L$ is polynomial when $\gamma$ is a constant, and $k=O(d)$.

Note that the argument in Lemma 1 can be similarly adapted leading to global convergence guarantee of the tensor power iteration step.

Two undercomplete, and one overcomplete component

Since in the SVD initialization Procedure 2, two components $\hat{a}^{(0)}$ and $\hat{b}^{(0)}$ are initialized through SVD, and the third component $\hat{c}^{(0)}$ is initialized through update formula (5), we can generalize the global convergence result in Theorem 2 to the setting where $A$, $B$ are undercomplete, and $C$ is overcomplete.

Then, if $k=O(d_{u})$ and $d_{o}\geq\operatorname{polylog}(k)$, the same convergence guarantee as in Theorem 2 holds.

We observe that given undercomplete modes $A$ and $B$, mode $C$ can be arbitrarily overcomplete, and we can still provide global recovery of $A,B$ and $C$ by employing SVD initialization procedure along modes $A$ and $B$.

3 Proof outline

The global convergence guarantee in Theorem 2 is established by combining the local convergence result in Theorem 1 and the SVD initialization result in Appendix C.

The local convergence result in Theorem 1 is derived by establishing error contraction in each iteration of the tensor power iteration and the coordinate descent for removing the residual error. Note that these convergence properties are broken down in Lemmata 1 and 12, respectively.

Since we assume generic factor matrices $A,B$ and $C$, we utilize many useful properties such as incoherence, bounded spectral norm of the matrices $A,B$ and $C$, bounded tensor spectral norm and so on. We list the precise set of deterministic conditions required to establish the local convergence result in Appendix A. Under these conditions, with a good initialization (i.e., small enough $\max\{\operatorname{dist}(\hat{a},a_{j}),\operatorname{dist}(\hat{b},b_{j})\}\leq\epsilon_{0}$), we show that the iterative update in (5) provides an estimate $\hat{c}$ with

for some contraction factor $q<1/2$. The incoherence condition is crucial for establishing this result. See Appendix B for the complete proof.

The initialization argument for SVD-based technique in Procedure 2 has two parts. The first part claims that by performing enough number of initializations (large enough $L$), a gap condition is satisfied, meaning that we obtain a vector $\theta$ which is relatively close to $c_{j}$ compared to any $c_{i},i\neq j$. This is a standard result for Gaussian vectors, e.g., see Lemma B.1 of Anandkumar et al. (2014a). In the second part of the argument, we analyze the dominant singular vectors of $T(I,I,\theta)$, for a vector $\theta$ with a good relative gap, to obtain an error bound on the initialization vectors. This is obtained through standard matrix perturbation results (Weyl and Wedin’s theorems). See Appendix C for the complete proof.

Experiments

For each initialization $\tau\in[L]$, an alternative option of running the algorithm with a fixed number of iterations $N$ is to stop the iterations based on some stopping criteria. In this experiment, we stop the iterations when the improvement in subsequent steps is small as

where $t_{\operatorname{S}}$ is the stopping threshold. According to the bound in Theorem 1, we set

Algorithm 1 is applied to random tensors with $d=1000$ and $k=\{10,50,100,200,500,1000,2000\}$. The number of initializations is $L=2000$. The parameter $t_{1}$ in (17) is fixed as $t_{1}=1e-08$. Figure 2 and Table 1 illustrate the outputs of running experiments which is the average of 10 random runs.

Figure 2 depicts the ratio of recovered columns versus the number of initializations. Both horizontal and vertical axes are plotted in $\log$-scale. We observe that it is much easier to recover the columns in the undercomplete settings ($k\leq d$), while it becomes harder when $k$ increases. Linear start in Figure 2 suggests that recovering the first bunch of columns only needs polynomial number of initializations. For highly undercomplete settings like $d=1000$ and $k=10$, almost all columns are recovered in this linear phase. After this start, the concave part means that it needs many more initializations for recovering the next bunch of columns. As we go ahead, it becomes harder to recover true columns, which is intuitive.

Table 1 has the results from the experiments. Parameters $k$, stopping threshold $t_{\operatorname{S}}$, and the average square error of the output, the average weight error and the average number of iterations are stated. The output averages are over several initializations and random runs. The square error is given by

for the corresponding recovered $j$. The error in estimating the weights is defined as $|\hat{w}-w_{j}|^{2}/w_{j}^{2}$ which is the square relative error of weight estimate. The number of iterations performed before stopping the algorithm is mentioned in the last column. We observe that by increasing $k$, all of these outputs are increased which means we get less accurate estimates with higher computation. This shows that recovering the overcomplete components is much harder. Note that by running the coordinate descent Algorithm 4, we can also remove this additional residual error left after the tensor power iteration step. Similar results and observations as above are seen when $k$ is fixed and $d$ is changed.

Running experiments with SVD initialization instead of random initialization yields nearly the same recovery rates, but with slightly smaller number of iterations. But, since the SVD computation is more expensive, in practice, it is desirable to initialize with random vectors. Our theoretical results for random initialization appear to be highly pessimistic compared to the efficient recovery results in our experiments. This suggests additional room for improving our theoretical guarantees under random initialization.

Acknowledgements

We acknowledge detailed discussions with Sham Kakade and Boaz Barak. We thank Praneeth Netrapalli for discussions on alternating minimization. We also thank Sham Kakade, Boaz Barak, Jonathan Kelner, Gregory Valiant and Daniel Hsu for earlier discussions on the $2\to p$ norm bound for random matrices, used in Lemma 3. We also thank Niranjan U.N. for discussions on running experiments. A. Anandkumar is supported in part by Microsoft Faculty Fellowship, NSF Career award CCF-$1254106$, NSF Award CCF-$1219234$, and ARO YIP Award W$911$NF-$13$-$1$-$0084$. M. Janzamin is supported by NSF Award CCF-1219234, ARO Award W911NF-12-1-0404 and ARO YIP Award W911NF-13-1-0084.

More Matrix Notations

Appendix A Deterministic Assumptions

In the main text, we assume matrices $A$, $B$, and $C$ are randomly generated. However, we are not using all the properties of randomness. In particular, we only need the following assumptions.

Rank-$k$ decomposition: The third order tensor $T$ has a CP rank of $k\geq 1$ with decomposition

where ${\cal S}^{d-1}$ denotes the unit $d$-dimensional sphere, i.e. all the vectors have unit This normalization is for convenience and the results hold for general case. $2$-norm as $\|a_{i}\|=\|b_{i}\|=\|c_{i}\|=1,i\in[k]$. Furthermore, define $w_{\min}:=\min_{i\in[k]}w_{i}$ and $w_{\max}:=\max_{i\in[k]}w_{i}$.

Incoherence: The components are incoherent, and let

for some $\alpha=\operatorname{polylog}(d)$. In other words, $A^{\top}A=I+J_{A}$, $B^{\top}B=I+J_{B}$, and $C^{\top}C=I+J_{C}$, where $J_{A}$, $J_{B}$, and $J_{C}$, are incoherence matrices with zero diagonal entries. We have $\max\left\{\|J_{A}\|_{\infty},\|J_{B}\|_{\infty},\|J_{C}\|_{\infty}\right\}\leq\rho$ as in (19).

Spectral norm conditions: The components satisfy spectral norm bound

Bounds on tensor norms: Tensor $T$ satisfies the bound

for some constant $\alpha_{0}$ and $\alpha=\operatorname{polylog}(d)$.

Rank constraint: The rank of the tensor is bounded by $k=o\left(d^{1.5}/\operatorname{polylog}d\right)$.

Bounded perturbation: Let $\psi$ denote the spectral norm of perturbation tensor as

Suppose $\psi$ is bounded as Note that for the local convergence guarantee, only the first condition $\psi\leq\frac{w_{\min}}{6}$ is required.

Weights ratio: The maximum ratio of weights $\gamma:=\frac{w_{\max}}{w_{\min}}$ satisfies the bound

Contraction factor: The contraction factor $q$ in Theorem 1 is defined as

for some constants $\alpha_{0},\beta^{\prime}>0$, and $\alpha=\operatorname{polylog}(d)$. In particular, we need $\alpha\alpha_{0}\sqrt{k}/d+\beta^{\prime}<w_{\max}/10w_{\min}$ which ensures $q<1/2$. This is satisfied when $\sqrt{k}/d<w_{\max}/w_{\min}\operatorname{poly}\log d$ and $\beta^{\prime}<w_{\max}/20w_{\min}$. The parameter $\beta^{\prime}$ is determined by the following assumption (initialization).

denote the initialization error w.r.t. to some $j\in[k]$. Suppose it is bounded as

for some constants $\alpha_{0},\beta^{\prime}>0$, $\alpha=\operatorname{polylog}(d)$, and $0<q<1/2$ which is defined in (21).

$2\to p$ norm: For some fixed constant $p<3$, $\max\{\|A^{\top}\|_{2\to p},\|B^{\top}\|_{2\to p},\|C^{\top}\|_{2\to p}\}\leq 1+o(1)$.

Many of the assumptions are actually parameter choices. The only properties of random matrices required are (A2), (A3), (A4) and (A10),. See Appendix A.1 for detailed discussion.

Let us provide a brief discussion about the above assumptions. Condition (A1) requires the presence of a rank-$k$ decomposition for tensor $T$. We normalize the component vectors for convenience, and this removes the scaling indeterminacy issues which can lead to problems in convergence. Additionally, we impose incoherence constraint in (A2), which allows us to provide convergence guarantee in the overcomplete setting. Assumptions (A3) and (A4) impose bounds on the spectral norm of tensor $T$ and its decomposition components. Note that assumptions (A2)-(A4) and (A10) are satisfied w.h.p. when the columns of $A$, $B$, and $C$ are generically drawn from unit sphere ${\cal S}^{d-1}$ (see Lemma 2 and Guédon and Rudelson (2007)), all others are parameter choices. Assumption (A5) limits the overcompleteness of problem which is required for providing convergence guarantees. The first bound on perturbation in (A6) as $\psi\leq\frac{w_{\min}}{6}$ is required for local convergence guarantee and the second bound $\psi\leq\frac{w_{\min}\sqrt{\log k}}{\alpha_{0}\sqrt{d}}$ is needed for arguing initialization provided by Procedure 2. Assumption (A7) is required to ensure contraction happens in each iteration. Assumption (A8) defines contraction ratio $q$ in each iteration, and Assumption (A9) is the initialization condition required for local convergence guarantee.

The tensor-spectral norm and $2\to p$ norm assumptions (A4) and (A10) may seem strong as we cannot even verify them given the matrix. However, when $k<d^{1.25-\epsilon}$ for arbitrary constant $\epsilon>0$, both conditions are implied by incoherence. See Lemma 4. We only need these assumptions to go to the very overcomplete setting.

Here, we provide arguments that random matrices satisfy conditions (A2), (A3), (A4), and (A10). It is well known that random matrices are incoherent, and have small spectral norm (bound on spectral norm dates back to Wigner (1955)). See the following lemma.

for some $\alpha=O(\sqrt{\log k})$ and $\alpha_{0}=O(1)$.

The spectral norm of the tensor is less well-understood. However, it can be bounded by the $2\to 3$ norm of matrices. Using tools from Guédon and Rudelson (2007); Adamczak et al. (2011), we have the following result.

This directly implies Assumption (A10). In particular, since we only apply Assumption (A10) to unsupervised setting ($k\leq O(d)$) in Appendix D, for randomly generated tensor, Assumption (A10) holds for all $p>2$ (notice that we only need it to hold for some $p<3$).

We also give an alternative proof of $2\to p$ norm which does not assume randomness and only relies on incoherence.

Proof: Let $L=\sqrt{d}/\operatorname{poly}\log d$. By incoherence assumption we know every subset of $L$ columns in $A$ has singular values within $1\pm o(1)$ (by Gershgorin Disk Theorem).

Here the second inequality uses that every entry outside $S$ is small, and last inequality uses the fact that $p>3-2\epsilon$. $\Box$

The $2\to 3$ norm implies a bound on the tensor spectral norm by Hölder’s inequality.

When $1/p+1/q=1$, for two sequence of numbers $\{a_{i}\},\{b_{i}\}$, we have

Consequently, we have the following corollary.

For vectors $f,g,h$, and weights $w_{i}\geq 0$, we have

Proof: The proof applies Hölder’s inequality twice as

where in the first application, $p=3$ and $q=3/2$, and in the second application, $p=q=2$ (which is the special case known as Cauchy-Schwartz). $\Box$

In the following lemma, it is shown that the first bound in Assumption (A4) holds for random matrices w.h.p.

Proof: For any unit vectors $\hat{a},\hat{b},\hat{c}$, we have

where Corollary 3 is exploited in the first inequality, and Lemma 3 is used in the last inequality. $\Box$

For the case with two undercomplete and one overcomplete dimensions (see Corollary 2), we can prove the tensor spectral norm using basic properties of the matrices $A,B,C$.

The first inequality uses triangle inequality and the fact that $|\langle c_{i},w\rangle|\leq 1$. The Cauchy-Schwartz inequality is exploited in the second inequality. Therefore, the spectral norm of the tensor is bounded by $O(w_{\max})$. $\Box$

Finally, we show in the following lemma that the second bound in Assumption (A4) is satisfied for random matrices.

where $\delta_{j}:=w_{j}\langle A_{i},A_{j}\rangle\langle B_{i},B_{j}\rangle$ is independent of $C_{j}$. From Lemma 2, columns of $A$ and $B$ are incoherent, and therefore, for $j\neq i$, we have

The proof is completed by applying union bound. $\Box$

For the convergence guarantees of the second step of algorithm on removing residual error, we need the following additional bound on the spectral norm of Khatri-Rao product of random matrices.

Spectral Norm Condition on Khatri-Rao Products: The components satisfy the following spectral norm bound on the Khatri-Rao products as

for $\alpha_{0}\leq\operatorname{poly}\log d$.

We now prove that Assumption (A11) is satisfied with high probability, if the columns of $A$, $B$ and $C$ are uniformly i.i.d. drawn from unit $d$-dimensional sphere.

The key idea is to view $(A\odot B)^{\top}(A\odot B)$ as the sum of random matrices, and use the following Matrix Bernstein’s inequality to prove concentration results.

Although the lemma requires all $M_{i}$’s to have spectral norm at most $R$ almost surely, it suffices to have spectral norm bounded by $R$ with high probability and bounded by $R^{\infty}=\mbox{poly}(d,k)$ almost surely. This is because we can always condition on the fact that $\|M_{i}\|\leq R$ for all $i$. Such conditioning can only change the expectations by a negligible amount, and does not affect independence between $M_{i}$’s.

Random unit vectors are not easy to work with, as entries in the same column are not independent. Thus, we first prove the result for matrices $A$ and $B$ whose entries are independent Gaussian variables.

Note that when $d<k$, by standard random matrix theory we know $\|Q\|\leq O(k)$. Also, every row of $Q$ has norm smaller than the corresponding row of $B^{\top}B$, which is bounded by $\|B\|\|b_{(i)}\|\leq O(\sqrt{kd})$. When $d\geq k$, again by matrix concentration we know $\|Q\|\leq O(\sqrt{dk\log d})$. Every row of $Q$ has norm bounded by $O(\sqrt{kd})$ (because entries in a row are independently random, with variance equal to $d$).

First let us bound the spectral norm for each of the $M_{i}$’s. Notice that for any vector $v$, $v^{\top}[(a_{i}a_{i}^{\top})*Q]v=(v*a_{i})^{\top}Q(v*a_{i})$ by definition of Hadamard product. On the other hand, $\|v*a_{i}\|\leq\|v\|\|a_{i}\|_{\infty}$. With high probability $\|a_{i}\|_{\infty}\leq O(\sqrt{\log k})$, hence $\|M_{i}\|\leq\|a_{i}\|_{\infty}^{2}\|Q\|$. This is bounded by $O(k\log d)$ when $d<k$ and $O(\sqrt{kd}\log^{2}d)$ when $k\leq d$.

By Matrix Bernstein we know with high probability $\|M\|\leq O(d\sqrt{k\log d})$. $\Box$

Using this lemma, it is easy to get a bound when columns of $A$, $B$ are unit vectors. In this case, we just need to normalize the columns, the normalization factor is bounded between $d^{2}/2$ and $2d^{2}$ with high probability, and therefore, $\|(A^{\top}A)(B^{\top}B)-I\|\leq O(\sqrt{k\log d}/d)$.

Appendix B Proof of Convergence Results in Theorems 1 and 2

The main part of the proof is to show that error contraction happens in each iteration of Algorithms 1 and 4 as the two main parts of the algorithm. Then, the contraction result after $t$ iterations is directly argued.

In the following, we first provide a local contraction result for the tensor power iteration (5) in Algorithm 1 given noisy tensor $\hat{T}$. This leads to Lemma 1 which is the local convergence guarantee of the tensor power updates. Then, we provide a local contraction argument for the coordinate descent step (9) in Algorithm 4.

Combining the above convergence arguments for both updates conclude the overall local convergence guarantee in Theorem. 1. Then, combing this local convergence guarantee and the initialization result in Theorem 3 leads to the global convergence guarantee in Theorem 2. In addition, the result in Corollary 2 is similarly argued where the bound on the spectral norm of the tensor is argued in Lemma 6.

In this section, we prove Lemma 1 which is the local convergence guarantee of the tensor power updates in Algorithm 1.

where $\alpha=\operatorname{polylog}(d)$ and $\alpha_{0}=O(1)$. Notice that this function is a small constant when $k<d^{1.5}/\operatorname{poly}\log d$.

Consider $\hat{T}=T+\Psi$ as the input to Algorithm 1, where $T$ is a rank-$k$ tensor, and $\Psi$ is a perturbation tensor. Suppose Assumptions (A1)-(A5) hold, and estimates $\hat{a}$ and $\hat{b}$ satisfy distance bounds

for some $j\in[k]$, and $\epsilon_{a},\epsilon_{b}>0$. Let $\epsilon:=\max\{\epsilon_{a},\epsilon_{b}\}$, and suppose $\psi$ defined in (20) be small enough such that This is the denominator of bound provided in (23).

where $f(\epsilon;k,d)$ is defined in (22). Then, update $\hat{c}$ in (5) satisfies the following distance bound with high probability (w.h.p.)

Furthermore, if the bound in (23) is such that $\operatorname{dist}(\hat{c},c_{j})\leq\epsilon$, then the update $\hat{w}:=\hat{T}(\hat{a},\hat{b},\hat{c})$ in (7) also satisfies w.h.p.

In the asymptotic regime, $f(\epsilon;k,d)$ is

The local convergence result provided in Theorem 1 has a linear convergence rate. But, Algorithm 1 actually provides an almost-quadratic convergence rate in the beginning, and linear convergence rate later on. It can be seen by referring to one-step contraction argument provided in Lemma 10 where the quadratic term $\alpha_{0}\epsilon^{2}$ exists. In the beginning, this term is dominant over linear term involving $\epsilon$, and we have almost-quadratic convergence. Writing $\alpha_{0}\epsilon^{2}=\alpha_{0}\epsilon^{\zeta}\epsilon^{2-\zeta}$, we observe that we get rate of convergence equal to $2-\zeta$ as long as we have initialization error bounded as $\epsilon_{0}^{\zeta}=O(1)$. Therefore, we can get arbitrarily close to quadratic convergence with appropriate initialization error. Note that when the model is more overcomplete, the algorithm more rapidly reaches to the linear convergence phase. For the sake of clarity, in proposing Theorem 1, we approximated the almost-quadratic convergence rate in the beginning with linear convergence.

Lemma 10 is proposed in the general form. In Lemma 11, we provide explicit contraction result by imposing additional perturbation, contraction and initialization Assumptions (A6), (A8) and (A9). We observe that under reasonable rank, perturbation and initialization conditions, the denominator in (23) can be lower bounded by a constant, and the numerator is explicitly bounded by a term involving $\epsilon$, and a constant non-contracting term.

Consider $\hat{T}=T+\Psi$ as the input to Algorithm 1, where $T$ is a rank-$k$ tensor, and $\Psi$ is a perturbation tensor. Let Assumptions As mentioned in the assumptions, from perturbation bound in (A6), only the bound $\psi\leq\frac{w_{\min}}{6}$ is required here. (A1)-(A9) hold. Note that initialization bound in (A9) is satisfied for some $j\in[k]$. Then, update $\hat{c}$ in (5) satisfies the following distance bound with high probability (w.h.p.)

and contraction ratio $q<1/2$ is defined in (21). Note that $\alpha=\operatorname{polylog}(d)$. In addition, if the above bound be such that $\operatorname{dist}(\hat{c},c_{j})\leq\epsilon_{0}$, then the update $\hat{w}:=\hat{T}(\hat{a},\hat{b},\hat{c})$ in (7) also satisfies w.h.p.

Proof of Lemma 1: We incorporate condition (A7) to show that $q<1/2$ in assumption (A8) is satisfied. In addition, (A7) implies that the bound on $\epsilon_{0}$ in assumption (A9) holds where it can be shown that the bound in (A9) is bounded as $O(1/\gamma)$. Then, the result is directly proved by iteratively applying the result of Lemma 11. $\Box$

Proof of auxiliary lemmata: tensor power iteration in Algorithm 1

Before providing the proofs, we remind a few definitions and notations.

In Assumption (A2), matrices $J_{A}$, $J_{B}$, and $J_{C}$, are defined as incoherence matrices with zero diagonal entries such that $A^{\top}A=I+J_{A}$, $B^{\top}B=I+J_{B}$, and $C^{\top}C=I+J_{C}$. We have $\max\left\{\|J_{A}\|_{\infty},\|J_{B}\|_{\infty},\|J_{C}\|_{\infty}\right\}\leq\rho$ as in (19).

Proof of Lemma 10: Let $z_{a}^{*}\perp a_{j}$ and $z_{b}^{*}\perp b_{j}$ denote the vectors that achieve supremum value in (12) corresponding to $\operatorname{dist}(\hat{a},a_{j})$ and $\operatorname{dist}(\hat{b},b_{j})$, respectively. Furthermore, without loss of generality, assume $\|z_{a}^{*}\|=\|z_{b}^{*}\|=1$. Then, $\hat{a}$ and $\hat{b}$ are decomposed as

Substituting $\hat{a}$ and $\hat{b}$ from (25) and (26), we have

where equalities $A^{\top}A=I+J_{A}$ and $B^{\top}B=I+J_{B}$ are exploited in the second equality, and the assumption that $z_{c}\perp\overline{C}_{j}$ is used in the last equality. The last inequality is from Assumption (A4). For $S_{2}$, we have

for some $\alpha=\operatorname{polylog}(d)$ and $\alpha_{0}=O(1)$. Second inequality is concluded from $\|u*v\|\leq\|u\|_{\infty}\cdot\|v\|,$ and Assumptions (A2) and (A3) are exploited in the last inequality. Similarly, for $S_{3}$, we have

for some $\alpha_{0}=O(1)$. The bound on $\|T\|$ is from Assumption (A4). Note that for random components, we showed in Lemma 5 that this bound holds w.h.p. exploiting Assumption (A5) and results of Guédon and Rudelson (2007). For the error term $\Psi(\hat{a},\hat{b},z_{c})$, we have

which is concluded from the definition of spectral norm of a tensor. Note that all vectors $\hat{a}$, $\hat{b}$, $z_{c}$ have unit norm.

Let $\epsilon:=\max\{\epsilon_{a},\epsilon_{b}\}$. Then, combining all the above bounds, we have w.h.p.

Now, we provide the bound on $|w_{j}-\hat{w}|$. As assumed in the lemma, we have distance bounds

The estimate $\hat{w}=\hat{T}(\hat{a},\hat{b},\hat{c})$ proposed in (7) can be expanded as

Proof of Lemma 11: The result is proved by applying Lemma 10, and incorporating additional conditions (A6), (A8), and (A9). $f(\epsilon_{0};k,d)$ in (22) can be bounded as

where $\epsilon_{0}\leq\frac{\beta^{\prime}}{\alpha_{0}}$ from Assumption (A9) is exploited in the inequality. The last equality is concluded from definition of contracting factor $q$ in (21). On the other hand, the denominator in (23) can be lower bounded as

where Assumptions (A9) and (A6) are used in the inequality. Applying Lemma 10, the result on $\operatorname{dist}(\hat{c},c_{j})$ is proved.

where $\epsilon_{0}\leq\frac{w_{\min}q}{4w_{\max}}$ from Assumption (A9) is used in the last inequality. $\Box$

B.2 Convergence of removing residual error: Algorithm 4

In this section, we provide convergence of the coordinate descent of Algorithm 4 for removing the residual error. We first provide the following definition.

Given an approximate solution $\{\widehat{A},\widehat{B},\widehat{C},\widehat{w}\}$, we call it $(\eta_{0},\eta_{1})$-nice if matrix $\widehat{A}$ (similarly $\widehat{B}$ and $\widehat{C}$) satisfies

Given above conditions are satisfied, we prove the following guarantees for removing residual error, Algorithm 4.

Consider $T$ as the input to Algorithm 4, where $T$ is a rank-$k$ tensor. Suppose Assumptions (A1)-(A5) and (A11) hold (which are satisfied whp when the components are uniformly i.i.d. drawn from unit $d$-dimensional sphere). Given initial solution $\left\{\widehat{A}^{(0)},\widehat{B}^{(0)},\widehat{C}^{(0)},\widehat{w}^{(0)}\right\}$ which is $(\eta_{0},\eta_{1})$-nice, all the following iterations of Algorithm 4 are $(2\eta_{0},3\eta_{1})$-nice. Furthermore, given the exact tensor $T$, the Frobenius norm error $\max\{\|\Delta A\|_{F},\|\Delta B\|_{F},\|\Delta C\|_{F},\|\Delta w\|/w_{\min}\}$ shrinks by at least a factor of 2 in every iteration. In addition, if we have a noisy tensor $\hat{T}=T+\Psi$ such that $\|\Psi\|\leq\psi$, then

Proof: iteration for removing residual error in Algorithm 4

We now prove Lemma 12 as the local convergence guarantee of the iterations for removing residual error, Algorithm 4.

To prove this lemma, we first observe that the algorithm update formula in (9) is (before normalization) $w_{i}\langle a_{i},\widehat{a}_{i}\rangle\langle b_{i},\widehat{b}_{i}\rangle c_{i}+\epsilon_{i}$ where

In the following lemma, we show that the error terms $\epsilon_{i}$’s are small.

Proof: By the update formula in (9), we know

We expand it into several terms as follows.

The norm of three different types of terms mentioned above are bounded in Section B, which conclude the desired bound in the lemma. $\Box$

When we have noise, all the $\epsilon_{i}$’s have an additional term $\Psi(\widehat{a}_{i},\widehat{b}_{i},I)$ which is bounded by $\psi$, and thus, the second part of the lemma follows directly.

For symmetric tensors we should change the algorithm as computing the following:

The result of this will be a change in the term of type 1. Now the Q matrix will be $(A\odot A)^{T}(A\odot A)-(1-\frac{1}{d})I-\frac{1}{d}J$ which has desired spectral norm for random matrices.

Claims for proving Lemma 13

The first term deals with the difference between $C$ and $\widehat{C}$.

Proof: This sum is equal to the Frobenius norm of a matrix $M=QZ$. Here the matrix $Q$ is a matrix such that is equal to $Q=(A\odot B)^{\top}(A\odot B)-I$:

The matrix $Z$ has columns $Z_{i}=w_{i}c_{i}-\widehat{w}_{i}\widehat{c}_{i}$. By assumption we know $\|Q\|\leq o(1)$, and $\|Z\|_{F}\leq w_{\max}\|\Delta C\|_{F}+\|\widehat{w}-w\|$. Therefore we have

Of course, in the error $\epsilon_{i}$, we don’t have $\sum_{j\neq i}\langle a_{i},a_{j}\rangle\langle b_{i},b_{j}\rangle w_{i}c_{i}$, instead we have terms like $\sum_{j\neq i}\langle\widehat{a}_{i},a_{j}\rangle\langle\widehat{b}_{i},b_{j}\rangle w_{i}c_{i}$. The next two lemmas show that these two terms are actually very close.

Same is true if any $\widehat{\cdot}$ is replaced by the true value.

Proof: Similar as before, we treat the left hand side as the Frobenius norm of some matrix $M=QZ$. Here $Z_{i}=\widehat{w}_{i}\widehat{c}_{i}$, and $Q$ is the following matrix:

Notice that the proof works for both terms. $\Box$

The same is true if the inner-products are between $\langle\Delta A_{j},\widehat{a}_{i}\rangle$ or $\langle\Delta B_{j},\widehat{b}_{i}\rangle$, or if any $\widehat{\cdot}$ is replaced by the true value.

Proof: Similar as before, we treat the left hand side as the Frobenius norm of some matrix $M=QZ$. Here $Z_{i}=\widehat{w}_{i}\widehat{c}_{i}$, and $Q$ is the following matrix

Now using definition of $2\to 4$ norm and $2ab\leq a^{2}+b^{2}$ we first bound the Frobenius norm of the matrix $Q$:

Now we first bound the $2\to 4$ norm of the matrix $\widehat{A}^{\top}=A^{\top}+\Delta A^{\top}$. By assumption we already know $\|A^{\top}\|_{2\to 4}\leq O(1)$. On the other hand, for any unit vector $u$

On the other hand we know $\|Z\|\leq O(w_{\max}\sqrt{k/d})$, hence $\|M\|_{F}\leq\|Z\|\|Q\|_{F}\leq o(w_{\max})(\|\Delta A\|_{F}+\|\Delta B\|_{F})$.

Projection Procedure 5

By construction it is clear that the columns of the new solution is within $\eta_{0}\sqrt{k}/d$ to the columns of the initial solution, so they must be within $2\eta_{0}\sqrt{k}/d$ to the columns of the true solution. The only thing left to prove is that $\|\widehat{A}\|\leq 3\eta_{1}\sqrt{k/d}$.

First we observe that $\widehat{A}=\widehat{A}^{0}+Z$ where $Z$ is a matrix whose columns are multiples of $Q-\widehat{A}^{0}$, and the multiplier is never larger than 1. Therefore $\|\widehat{A}\|\leq\|h{A}^{0}\|+\|Z\|\leq\|\widehat{A}^{0}\|+\|Q-\widehat{A}^{0}\|\leq 2\|\widehat{A}^{0}\|+\|Q\|\leq 3\eta_{1}\sqrt{k/d}$. $\Box$

Appendix C SVD Initialization Result

In this section, we analyze the SVD-based initialization technique proposed in Procedure 2. The goal is to provide good initialization vectors close to the columns of true components $A$ and $B$ in the regime of $k=O(d)$.

Since $A$ and $B$ are not orthogonal matrices, the expansion in (29) is not the SVD Note that if $A$ and $B$ are orthogonal matrices, columns of $A$ and $B$ are directly recovered by computing SVD of $T(I,I,\theta)$. of $T(I,I,\theta)$. But, we show in the following theorem that if we draw enough number of random vectors $\theta$ in the regime of $k=O(d)$, we can eventually provide good initialization vectors through SVD of $T(I,I,\theta)$. Define

Consider tensor $\hat{T}=T+\Psi$ where $T$ is a rank-$k$ tensor, and $\Psi$ is a perturbation tensor. Let Assumptions (A1)-(A3) hold and $k=O(d)$. Draw $L$ i.i.d. random vectors $\theta^{(j)}\sim\mathcal{N}(0,I_{d}),j\in[L]$. Let $u_{1}^{(j)}$ and $v_{1}^{(j)}$ be the top left and right singular vectors of $\hat{T}(I,I,\theta^{(j)})$. This is $L$ random runs of Procedure 2. Suppose $L$ satisfies the bound

where $\psi:=\|\Psi\|$ is the spectral norm of perturbation tensor $\Psi$, and $\alpha_{0}>1$ is a constant.

From (30), with probability at least $1-2k^{-1}$, we have

From (31), with probability at least $1-k^{-7}$, we have

In the following Lemma, we show that the gap condition between the maximum and the second maximum of vector $\lambda$ required in Lemma 16 is satisfied under some number of random draws.

for some $0<\mu<\frac{w_{\min}}{w_{\max}\rho}-1$. Then, with probability at least $1-2k^{-1}-k^{-7}$, we have the following gap condition for at least one draw, say $j^{*}$,

From the given bound on $L$ in the lemma and inequalities (30) and (31), with probability at least $1-2k^{-1}-k^{-7}$, we have

where $\alpha=\operatorname{polylog}(d)$, and $\alpha_{0}>0$ is a constant.

Consider $\hat{T}=T+\Psi$, where $T$ is a rank-$k$ tensor, and $\Psi$ is a perturbation tensor. Let assumptions (A1)-(A3) hold for $T$. Let $u_{1}$ and $v_{1}$ be the top left and right singular vectors of $\hat{T}(I,I,\theta)$. Let

denote the vector that captures correlation of $\theta$ with different $c_{i},i\in[k]$, weighted by $w_{i},i\in[k]$. Without loss of generality, assume that $\lambda_{1}=\max_{i}|\lambda_{i}|$, and let $\lambda_{(2)}:=\max_{i\neq 1}|\lambda_{i}|$. Suppose the relative gap condition

is satisfied for some $\mu>\frac{\lambda_{1}}{\lambda_{1}-\|\Psi(I,I,\theta)\|}2\mu_{R}-1$, where $\mu_{R}$ and $\mu_{\min}$ are defined in (32). Then, with high probability (w.h.p.),

Proof: From Assumption (A1), $T(I,I,\theta)$ can be written as equation (29), Expanded as

From here, we prove the result in two cases. First when $\mu_{E}<\mu_{R}$ and therefore $\mu_{\min}=\mu_{E}$, and second when $\mu_{E}\geq\mu_{R}$ and therefore $\mu_{\min}=\mu_{R}$.

Case 1 ($\mu_{E}<\mu_{R}$): According to the subspaces spanned by $a_{1}$ and $b_{1}$, we decompose matrix $R$ to two components as $R={\cal P}_{\perp}(R)+{\cal P}_{\parallel}(R)$. First term ${\cal P}_{\perp}(R)$ is the component with column space orthogonal to $a_{1}$ and row space orthogonal to $b_{1}$, and ${\cal P}_{\parallel}(R)$ is the component with either the column space equal to $a_{1}$ or the row space equal to $b_{1}$. We have

Looking at $M$, it becomes more clear why we proposed the above decomposition for $R$. Since the column and row space of ${\cal P}_{\perp}(R)$ are orthogonal to $a_{1}$ and $b_{1}$, respectively, the SVD of $M$ has $a_{1}$ and $b_{1}$ as its left and right singular vectors, respectively. Hence, $M$ has the SVD form

where $u_{1}$ and $v_{1}$ are its top left and right singular vectors. We have

where the sub-multiplicative property of spectral norm is used in the second inequality, and the last inequality is from Assumption (A3). From Weyl’s theorem, we have

where (36) is used in the second inequality. Therefore, we have

Bounding the spectral norm of $E$: For any $i\neq j$, let $\rho_{ij}^{(a)}:=|\langle a_{i},a_{j}\rangle|$ and $\rho_{ij}^{(b)}:=|\langle b_{i},b_{j}\rangle|$. We have

Where the first equality is concluded from Lemma 19, and Assumptions (A2) and (A3) are exploited in the last inequality. Similarly, for $E_{2}$ and $E_{3}$, we have

Case 2 ($\mu_{R}\leq\mu_{E}$): The result can be similarly achieved when $\mu_{R}\leq\mu_{E}$. Here we directly apply Wedin’s theorem to $\hat{T}(I,I,\theta)=\lambda_{1}a_{1}b_{1}^{\top}+R+\Psi(I,I,\theta)$, treating $R+\Psi(I,I,\theta)$ as the error term. From Weyl’s theorem, we have

The above lemma concludes the proof for initialization procedure, except for a few auxiliary lemmata that we prove next.

First we use Gaussian tail bounds to prove that the largest entry of a Gaussian vector can be quite large with inverse polynomial probability:

Let $x\sim\mathcal{N}(0,\sigma)$ be a Gaussian random variable with mean zero and variance $\sigma^{2}$. Then, for any $t>0$, we have

where $f(t)=\frac{1}{\sqrt{2\pi}}e^{-t^{2}/2}$.

Proof: Let $z=\frac{x}{\sigma}$, where $z\sim\mathcal{N}(0,1)$ is a standard Gaussian random variable. Then, we have $\Pr[x\geq t]=\Pr[z\geq t/\sigma]$, and therefore, the result is proved by using standard tail bounds for Gaussian random variable. $\Box$

Proof: From Lemma 17, for any $i\in[k]$, we have

where the last inequality is concluded from the fact that $k\geq 2$. The result is then proved by taking a union bound. $\Box$

Next we prove a basic fact about spectral norm that is used in the proof of Lemma 16.

We have $Hx=\langle v,x\rangle h$, and therefore, $\|Hx\|=|\langle v,x\rangle|\|h\|$. This is maximized by $x=v/\|v\|$, and this finishes the proof. $\Box$

Finally, we show that noise matrix $\Psi(I,I,\theta)$ has bounded norm with high probability which is useful for initialization argument in Theorem 3.

where $\psi:=\|\Psi\|$ is the spectral norm of error tensor $\Psi$.

Proof: Let $\theta_{n}:=\frac{1}{\|\theta\|}\theta$ denote the normalized version of $\theta$. Then, we have

where the last inequality is from the definition of tensor spectral norm. Applying the bound on $\|\theta\|$ in Lemma 21 finishes the proof. $\Box$

The following lemma provides concentration bound for the norm of standard Gaussian vector which is basically a tail bound for the chi-squared random variable.

Let the random vector $\theta$ is distributed as $\mathcal{N}(0,I_{d})$. Then, for any $\alpha_{0}>1$, we have

Appendix D Clustering Process

In the last step of main algorithm, we need to cluster the generated 4-tuples into $k$ clusters. Theoretically, we only have convergence guarantees when the initialization vectors are good enough, while the other initializations can potentially generate arbitrary 4-tuples. In the worst case, these arbitrary 4-tuples can make the clustering process hard, and therefore, we provide specific Procedure 3 for which the output properties are provided in Lemma 24.

Note that the key observation for the algorithm is if $T(\hat{a},\hat{b},\hat{c})$ is large for some $(\hat{a},\hat{b},\hat{c})$, then these vectors are close to $(a_{i},b_{i},c_{i})$ for some $i\in[k]$.

For simplicity, we only prove this when the initialization procedure in Theorem 2 takes polynomial time, namely $k=O(d)$ and $w_{\max}/w_{\min}=O(1)$. Without loss of generality, we also assume $w_{\max}=w_{1}\geq w_{2}\geq\cdots\geq w_{k}=w_{\min}$. In this case, we choose the threshold $\epsilon$ in the following lemmata to be some small constant depending on $k/d$ and $w_{\max}/w_{\min}$. Also, we work in the case when noise $\Psi=0$, however the proof still works when the noise $\psi=\|\Psi\|=o(1)$.

for some $t\in[k]$. Let $\delta:=O\left(\frac{w_{\max}}{w_{\min}}\epsilon^{3-p}\right)$, and assume $|T(\hat{a},\hat{b},\hat{c})|\geq(1-\delta)w_{t}$. Then, there exists some $j$ such that

Proof: Partition tensor $T=\sum_{i\in[k]}w_{i}a_{i}\otimes b_{i}\otimes c_{i}$ to $T_{1}+T_{2}$, where $T_{1}$ contains all the terms indexed from $1$ to $t-1$, and $T_{2}$ contains the remaining terms. From Corollary 3, we have

where Assumption (A10) and the assumption in the lemma are exploited in the last step. Similar arguments hold for $b$ and $c$. Combining with the earliest inequality, we have

where the definition of $\delta$ is exploited in the last inequality. Applying assumption $|T(\hat{a},\hat{b},\hat{c})|\geq(1-\delta)w_{t}$ to the above bound, we have

Since all the 3-norms are bounded by $1+o(1)$, each of them must be at least $1-O(\delta)$ to let inequality (37) hold. Now we have

where the last inequality is from Assumption (A10). This implies $\max\{|\langle a_{j},\hat{a}\rangle|\}=1-O(\delta)$, which in turn implies there exists a $j$ such that

when $\epsilon$ and $\delta$ are small enough.

By symmetry we know there is also a $j^{\prime}$ such that $\operatorname{dist}(\hat{b},b_{j^{\prime}})<w_{\min}/10w_{\max}$. If $j\neq j^{\prime}$, then it is easy to check $T_{2}(\hat{a},\hat{b},\hat{c})$ cannot be large. Hence, $j=j^{\prime}$ and the Lemma is correct. $\Box$

On the other hand, we know if there is a good initialization, the largest $T(\hat{a},\hat{b},\hat{c})$ must be large.

Suppose there exists a good initialization (see initialization condition (13) in the local convergence theorem) for some column $t\in[k]$, and

Let $\delta:=O\left(\frac{w_{\max}}{w_{\min}}\epsilon^{3-p}\right)$. Then the corresponding output of iterations in Algorithm 1 denoted by $(\hat{a},\hat{b},\hat{c})$ satisfy

Furthermore, for any $i\neq t$, $\max\{|\langle\hat{a},a_{i}\rangle|,|\langle\hat{b},b_{i}\rangle|,|\langle\hat{c},c_{i}\rangle|\}\leq o(\epsilon)$.

Proof: Similar to the proof of Lemma 22, partition tensor $T=\sum_{i\in[k]}w_{i}a_{i}\otimes b_{i}\otimes c_{i}$ to $T_{2}=w_{t}a_{t}\otimes b_{t}\otimes c_{t}$ and $T_{1}=T-T_{2}$. Since the initialization is good, by the local convergence result in Theorem 1, we have

where the incoherence condition and $p>2$ are exploited in the last step. Therefore, $|T_{2}(\hat{a},\hat{b},\hat{c})|\geq(1-\delta/2)w_{t}$.

Similar to Lemma 22, by using Corollary 3, we have $|T_{1}(\hat{a},\hat{b},\hat{c})|\leq w_{t}\delta/2$. Applying these bounds, we have

The last part of the Lemma is trivial because $\operatorname{dist}(\hat{a},a_{t})$ is small and $\langle a_{i},a_{t}\rangle$ is small by incoherence. $\Box$

Finally we prove the clustering process succeeds.

Proof: We prove by induction to show that every step of the algorithm correctly computes one component.

By Lemma 23 we know there must be a 4-tuple with $|T(\hat{a},\hat{b},\hat{c})|>w_{t}(1-\delta)$. On the other hand, by Lemma 22 we know the 4-tuple we found must satisfy $\max\{\operatorname{dist}(\hat{a},a_{j}),\operatorname{dist}(\hat{b},b_{j})\}<w_{\min}/10w_{\max}$ for some $j$ (and this cannot be some $j$ that has already been found). This tuple then satisfies the conditions of the local convergence Theorem 1. Hence, after $N$ iterations it must have converged to $(a_{j},b_{j},c_{j})$. At this step the algorithm successfully found a new component of the tensor.