Provable Tensor Factorization with Missing Data

Introduction

Several real-world applications routinely encounter multi-way data with structure which can be modeled as low-rank tensors. Moreover, in several settings, many of the entries of the tensor are missing, which motivated us to study the problem of low-rank tensor factorization with missing entries. For example, when recording electrical activities of the brain, the electroencephalography (EEG) signal can be represented as a three-way array (temporal, spectral, and spatial axis). Oftentimes signals are lost due to mechanical failure or loose connection. Given numerous motivating applications, several methods have been proposed for this tensor completion problem. However, with the exception of 2-way tensors (i.e., matrices), the existing methods for higher-order tensors do not have theoretical guarantees and typically suffer from the curse of local minima.

In general, finding a factorization of a tensor is an NP-hard problem, even when all the entries are available. However, it was recently discovered that by restricting attention to a sub-class of tensors such as low-CP rank orthogonal tensors or low-CP rank incoherentThe notion of incoherence we assume in (2) can be thought of as incoherence between the fibers and the standard basis vectors. tensors , one can efficiently find a provably approximate factorization. In particular, exact recovery of the factorization is possible for a tensor with a low-rank orthogonal CP decomposition . We ask the question of recovering such a CP-decomposition when only a small number of entries are revealed, and show that exact reconstruction is possible even when we do not observe any entry in most of the fibers.

Problem formulation. We study tensors that have an orthonormal CANDECOMP/PARAFAC (CP) tensor decomposition with a small number of components. Moreover, for simplicity of notation and exposition, we only consider symmetric third order tensors. We would like to stress that our techniques generalizes easily to handle non-symmetric tensors as well as higher-order tensors. Formally, we assume that the true tensor $T$ has the the following form:

where $[n]=\{1,\ldots,n\}$ is the set of the first $n$ integers. Tensor completion becomes increasingly difficult for tensors with larger $\mu(T)$, because the ‘mass’ of the tensor can be concentrated on a few entries that might not be revealed. Out of $n^{3}$ entries of $T$, a subset $\Omega\subseteq[n]\times[n]\times[n]$ is revealed. We use ${\cal P}_{\Omega}(\cdot)$ to denote the projection of a matrix onto the revealed set such that

We want to recover $T$ exactly using the given entries ($P_{\Omega}(T)$). We assume that each $(i,j,k)$ for all $i\leq j\leq k$ is included in $\Omega$ with a fixed probability $p$ (since $T$ is symmetric, we include all permutations of $(i,j,k)$). This is equivalent to fixing the total number of samples $|\Omega|$ and selecting $\Omega$ uniformly at random over all ${n^{3}\choose|\Omega|}$ choices. The goal is to ensure exact recovery with high probability and for $|\Omega|$ that is sub-linear in the number of entries ($n^{3}$).

Ideally, one would like to minimize the rank of a tensor that explains all the sampled entries.

Recently, showed that an alternating minimization technique, can recover a matrix with missing entries exactly. We generalize and modify the algorithm for the case of higher order tensors and study it rigorously for tensor completion. However, due to special structure in higher-order tensors, our algorithm as well as analysis is significantly different than the matrix case (see Section 2.2 for more details).

The main novelty in our approach is that we refine all $r$ components iteratively as opposed to the sequential deflation technique used by the existing methods for tensor decomposition (for fully observed tensors). In sequential deflation methods, components $\{{\bf u}_{1},{\bf u}_{2},\dots,{\bf u}_{r}\}$ are estimated sequentially and estimate of say ${\bf u}_{2}$ is not used to refine ${\bf u}_{1}$. In contrast, our algorithm iterates over all $r$ estimates in the inner loop, so as to obtain refined estimates for all ${\bf u}_{i}$’s in the outer loop. We believe that such a technique could be applied to improve the error bounds of (fully observed) tensor decomposition methods as well.

As our method is directly solving a non-convex problem, it can easily get stuck in local minima. The key reason our approach can overcome the curse of local minima is that we start with a provably good initial point which is only a small distance away from the optima. To obtain such an initial estimate, we compute a low-rank approximation of the observed tensor using Robust Tensor Power Method (RTPM) . RTPM is a generalization of the widely used power method for computing leading singular vectors of a matrix and can approximate the largest singular vectors up to the spectral norm of the “error” tensor. Hence, the challenge is to show that the error tensor has small spectral norm (see Theorem 2.1). We perform a thresholding step similar to (see Lemma A.4) after the RTPM step to ensure that the estimates we get are incoherent, which is critical for our analysis.

Our analysis requires the sampled entries $\Omega$ to be independent of the current iterates ${\bf u}_{i},\forall i$, which in general is not possible as ${\bf u}_{i}$’s are computed using $\Omega$. To avoid this issue, we divide the given samples ($\Omega$) into equal $r\cdot\tau$ parts randomly where $\tau$ is the number of outer loops (see Algorithm 1).

2 Main Result

then the following holds with probability at least $1-n^{-5}\log_{2}(4\sqrt{r}\,\|T\|_{F}/\varepsilon)$:

the problem (4) has a unique optimal solution; and

$\log_{2}(\frac{4\sqrt{r}\,\|T\|_{F}}{\varepsilon})$ iterations of Algorithm 1 produces an estimate $\widehat{T}$ s.t. $\|T-\widehat{T}\|_{F}\leq\varepsilon\;.$

Note that the above result can be generalized to $k$-mode tensors in a straightforward manner, where exact recovery is guaranteed if, $p\geq C\,\frac{\mu^{6}\,r^{5}\,\sigma_{\rm max}^{4}\,(\log n)^{4}\,\log(r\|T\|_{F}/\varepsilon)}{\sigma_{\rm min}^{4}\,n^{k/2}}$. However, for simplicity of notation and to emphasize key points of our proof we present our proof for $3$-mode tensors only in Section 2.3.

We provide a proof of Theorem 1.1 in Section 2. For an incoherent, well-conditioned, and low-rank tensor with $\mu=O(1)$ and $\sigma_{\rm min}=\Theta(\sigma_{\rm max})$, alternating minimization requires $O(r^{5}n^{3/2}(\log n)^{4})$ samples to get within an arbitrarily small normalized error. This is a vanishing fraction of the total number of entries $n^{3}$. Each step in the alternating minimization requires $O(r|\Omega|)$ operations, hence the alternating minimization only requires $O(r|\Omega|\log(r\|T\|_{F}/\varepsilon))$ operations. The initialization step requires $O(r^{c}|\Omega|)$ operations for some positive constant numerical $c$. When $r\ll n$, the computational complexity scales linearly in the sample size up to a logarithmic factor.

A fiber in a third order tensor is an $n$-dimensional vector defined by fixing two of the axes and indexing over remaining one axis. The above theorem implies that among $n^{2}$ fibers of the form $\{T[{\mathds{I}},e_{j},e_{k}]\}_{j,k\in[n]}$, it is sufficient to have only $O(n^{3/2}(\log n)^{4})$ fibers with any samples. Most of the fibers are not sampled at all and, perhaps surprisingly, our approach can still recover the original low-rank tensor. This should be compared to the matrix completion setting where all fibers are required to have at least one sample.

However, unlike matrices, the fundamental limit of higher order tensor completion is not known. Building on the percolation of Erdös-Renýi graphs and the coupon-collectors problem, it is known that matrix completion has multiple rank-$r$ solutions when the sample size is less than $C\mu rn\log n$ , hence exact recovery is impossible. But, such arguments do not generalize directly to higher order; see Section 2.5 for more discussion. Interestingly, simulations in Section 1.3 suggests that for $r=O(\sqrt{n})$, the sample complexity scales as $(r^{1/2}n^{3/2}\log n)$. That is, assuming the sample complexity provided by simulations is correct, our result achieves optimal dependence on $n$ (up to $\log$ factors). However, the dependency on $r$ is sub-optimal (see Section 2.5 for a discussion).

3 Empirical Results

Theorem 1.1 guarantees exact recovery when $p\geq Cr^{5}(\log n)^{4}/n^{3/2}$. Numerical experiments show that the average recovery rate converges to a universal curve over $\alpha$, where $p^{*}=\alpha r^{1/2}\ln n/((1-\rho)n^{3/2})$ in Figure 1. Our bound is tight in its dependency $n$ up to a poly-logarithmic factor, but is loose in its dependency in the rank $r$. Further, it is able to recover the original matrix exactly even when the factors are not strictly orthogonal.

4 Related Work

Tensor decomposition and completion: The CP model proposed in is a multidimensional generalization of singular value decomposition of matrices. Computing the CP decomposition involves two steps: first apply a whitening operator to the tensor to get a lower dimensional tensor with orthogonal CP decomposition. Such a whitening operator only exists when $r\leq n$. Then, apply known power-method techniques for exact orthogonal CP decomposition . We use this algorithm as well as the analysis for the initial step of our algorithm. For motivation and examples of orthogonal CP models we refer to .

Recently, many heuristics for tensor completion have been developed such as the weighted least squares , Gauss-Newton , alternating least-squares , trace norm minimization . However, to the best of our knowledge, there is no tensor completion method with provable guarantees. In a different context, show that minimizing a weighted trace norm of flattened tensor provides exact recovery using $O(rn^{3/2})$ samples, but each observation needs to be a dense random projection of the tensor as opposed to observing just a single entry, which is the case in the tensor completion problem.

Relation to matrix completion: Matrix completion has been studied extensively in the last decade since the seminal paper by Candes and Recht . Since then, several provable approaches have been developed, such as, nuclear norm minimization , OptSpace , and Alternating Minimization . However, several aspects of tensor factorization makes it challenging to adopt matrix completion algorithms and analysis techniques directly.

First, there is no natural convex surrogate of the tensor rank function and developing such a function is in fact a topic of active research Next, even when all entries are revealed, tensor decomposition methods such as simultaneous power iteration are known to get stuck at local extrema, making it challenging to apply matrix decomposition methods directly. Third, for the initialization step, the best low-rank approximation of a matrix is unique and finding it is trivial. However, for tensors, finding the best low-rank approximation is notoriously difficult.

On the other hand, some aspects of tensor decomposition makes it possible to prove stronger results. Matrix completion aims to recover the underlying matrix only, since the factors are not uniquely defined due to invariance under rotations. However, for orthogonal CP models, we can hope to recover the individual singular vectors ${\bf u}_{i}$’s exactly. In fact, Theorem 1.1 shows that our method indeed recovers the individual singular vectors exactly.

Spectral analysis of tensors and hypergraphs: Theorem 2.1 and Lemma 2.2 should be compared to copious line of work on spectral analysis of matrices , with an important motivation of developing fast algorithms for low-rank matrix approximations. We prove an analogous guarantee for higher order tensors and provide a fast algorithm for low-rank tensor approximation. Theorem 2.1 is also a generalization of the celebrated result of Friedman-Kahn-Szemerédi and Feige-Ofek on the second eigenvalue of random graphs. We provide an upper bound the largest second eigenvalue of a random hypergraph, where each edge includes three nodes and each of the ${n\choose 3}$ edges is selected with probability $p$.

Analysis of the Alternating Minimization Algorithm

In this section, we provide a proof of Theorem 1.1 and the proof sketches of the required main technical theorems. Formal proofs of the technical theorems and lemmas are provided in the appendix. There are two key components: $a$) the analysis of the initialization step (Section 2.1); and $b$) the convergence of alternating minimization given a sufficiently accurate initialization (Section 2.2). We use these two analyses to prove Theorem 1.1 in Section 2.3.

We first show that $(1/p){\cal P}_{\Omega}(T)$ is close to $T$ in spectral norm, and use it bound the error of robust power method applied directly to ${\cal P}_{\Omega}(T)$. The normalization by $(1/p)$ compensates for the fact that many entries are missing. For a proof of this theorem, we refer to Appendix A.

For $p=\alpha/n^{3/2}$ satisfying $\alpha\geq\log n$, there exists a positive constant $C>0$ such that, with probability at least $1-n^{-5}$,

where $T_{\rm max}\equiv\max_{i,j,k}T_{ijk}$, and $\|T\|_{2}\equiv\max_{\|u\|=1}T[u,u,u]$ is the operator norm.

Notice that $T_{\rm max}$ is the maximum entry in the tensor $T$ and the factor $1/(T_{\rm max}n^{3/2}p)$ corresponds to normalization with the worst case operator norm of $p\,T$, since $\|pT\|_{2}\leq T_{\rm max}n^{3/2}p$ and the maximum is achieved by $T=T_{\rm max}({\mathds{1}}\otimes{\mathds{1}}\otimes{\mathds{1}})$. The following theorem guarantees that $O(n^{3/2}(\log n)^{2})$ samples are sufficient to ensure that we get arbitrarily small error. A formal proof is provided in Appendix.

Together with an analysis of robust tensor power method [1, Theorem 5.1], the next error bound follows from directly substituting (6) and using the fact that for incoherent tensors $T_{\rm max}\leq\sigma_{\rm max}\mu(T)^{3}r/n^{3/2}$. Notice that the estimates can be computed efficiently, requiring only $O(\log r+\log\log\alpha)$ iterations, each iteration requiring $O(\alpha n^{3/2})$ operations. This is close to the time required to read the $|\Omega|\simeq\alpha n^{3/2}$ samples. One caveat is that we need to run robust power method ${\rm poly}({r}\log n)$ times, each with fresh random initializations.

2 Alternating Minimization Analysis

We now provide convergence analysis for the alternating minimization part of Algorithm 1 to recover rank-$r$ tensor $T$. Our analysis assumes that $\|{\bf u}_{i}-{\bf u}^{*}_{i}\|_{2}\leq c\sigma_{\rm min}/r\sigma_{\rm max},\forall i$ where $c$ is a small constant (dependent on $r$ and the condition number of $T$). The above mentioned assumption can be satisfied using our initialization analysis and by assuming $\Omega$ is large-enough.

At a high-level, our analysis shows that each step of Algorithm 1 ensures geometric decay of a distance function (specified below) which is “similar” to $\max_{j}\|{\bf u}_{j}^{t}-{\bf u}^{*}_{j}\|_{2}$.

The next theorem shows that this distance function decreases geometrically with number of iterations of Algorithm 1. A proof of this theorem is provided in Appendix B.4.

If $d_{\infty}([U,\Sigma],[U^{*},\Sigma^{*}])\leq\frac{1}{1600r}\frac{\sigma^{*}_{min}}{\sigma^{*}_{max}}$ and ${\bf u}_{i}$ is $2\mu$-incoherent for all $1\leq i\leq r$, then there exists a positive constant $C$ such that for $p\geq\frac{Cr^{2}(\sigma^{*}_{\rm max})^{2}\mu^{3}\log^{2}n}{(\sigma_{\rm min}^{*})^{2}n^{3/2}}$ we have w.p. $\geq 1-\frac{1}{n^{7}}$,

Note that our number of samples depend on the number of iterations $\tau$. But due to linear convergence, our sample complexity increases only by a factor of $\log(1/\epsilon)$ where $\epsilon$ is the desired accuracy.

Difference from Matrix AltMin: Here, we would like to highlight differences between our analysis and analysis of the alternating minimization method for matrix completion (matrix AltMin) . In the matrix case, the singular vectors ${\bf u}^{*}_{i}$’s need not be unique. Hence, the analysis is required to guarantee a decay in the subspace distance $dist(U,U^{*})$; typically, principal angle based subspace distance is used for analysis. In contrast, orthonormal ${\bf u}^{*}_{i}$’s uniquely define the tensor and hence one can obtain distance bounds $\|{\bf u}_{i}-{\bf u}^{*}_{i}\|_{2}$ for each component ${\bf u}_{i}$ individually.

3 Proof of Theorem 1.1

Let $T=\sum_{q=1}^{r}\sigma^{*}_{q}(u^{*}_{q}\otimes u^{*}_{q}\otimes u^{*}_{q})$. Denote the initial estimates $U^{0}=[u^{0}_{1},\ldots,u^{0}_{r}]$ and $\sigma^{0}=[\sigma^{0}_{1},\ldots,\sigma^{0}_{r}]$ to be the output of robust tensor power method at step 5 of Algorithm 1. With a choice of $p\geq C(\sigma_{\rm max}^{*})^{4}\mu^{6}r^{4}(\log n)^{4}/(\sigma_{\rm min}^{*})^{4}n^{3/2}$ as per our assumption, Lemma 2.2 ensures that we have $\|u^{0}_{q}-u^{*}_{q}\|\leq\sigma^{*}_{\rm min}/(4800\,r\sigma_{\rm max})$ and $|\sigma_{q}^{0}-\sigma_{q}^{*}|\leq|\sigma_{q}^{*}|\sigma^{*}_{\rm min}/(4800\,r\sigma_{\rm max})$ with probability at least $1-n^{-5}$. This requires running robust tensor power method for $(r\log n)^{c}$ random initializations for some positive constant $c$, each requiring $O(|\Omega|)$ operations ignoring logarithmic factors.

With this initialization, Theorem 2.3 tells us that $O(\log_{2}(4r^{1/2}\|T\|_{F}/\varepsilon)$ iterations (each iteration requires $O(r|\Omega|)$ operations) is sufficient to achieve:

4 Fundamental limit and random hypergraphs

For matrices, it is known that exact matrix completion is impossible if the underlying graph is disconnected. A refinement of this analysis for Erdös-Renýi graphs provides a lower bound on matrix completion: when sample size is less than $C\mu rn\log n$, no algorithm can recover the original matrix .

However, for tensor completion and random hyper graphs, such a simple connection does not exist. It is not known how the properties of the hyper graph is related to recovery. In this spirit, a rank-one third-order tensor completion has been studied in a specific context of MAX-3LIN problems. Consider a series of linear equations over $n$ binary variables $x=[x_{1}\,\ldots\,x_{n}]\in\{\pm 1\}^{n}$. An instance of a 3LIN problem consists of a set of linear equations on GF(2), where each equation involve exactly three variables, e.g.

We use $-1$ to denote true (or $1$ in GF(2)) and $+1$ to denote false (or in GF(2)). Then the exclusive-or operation denoted by $\oplus$ is the integer multiplication. the MAX-3LIN problem is to find a solution $x$ that satisfies as many number of equations as possible. This is an NP-hard problem in general, and hence random instances of the problem with a planted solution has been studied . Algorithm 1 provides a provable guarantee for MAX-3LIN with random assignments.

For random MAX-3LIN problem with a planted solution, under the hypotheses of Theorem 1.1, Algorithm 1 finds the correct solution with high probability.

Notice that the tensor has incoherence one and rank one. This implies exact reconstruction for $P\geq C(\log n)^{4}/n^{3/2}$. This significantly improves over a message-passing approach to MAX-3LIN in , which is guaranteed to find the planted solution for $p\geq C(\log\log n)^{2}/(n\log n)$. It was suggested that a new notion of connectivity called propagation connectivity is a sufficient condition for the solution of random MAX-3LIN problem with a planted solution to be unique [23, Proposition 2]. Precisely, it is claimed that if the hypergraph corresponding to an instance of MAX-3LIN is propagation connected, then the optimal solution for MAX-3LIN is unique and there is an efficient algorithm that finds it. However, the example in 7 is propagation connected but there is no unique solution: $,,,$ all satisfy the equations and corresponding tensor $T=x\otimes x\otimes x$ is not uniquely defined either. This proves that propagation connectivity is not a sufficient condition for uniqueness of the MAX-3LIN solution.

5 Open Problems and Future Directions

Tensor completion for non-orthogonal decomposition. Numerical simulations suggests that non-orthogonal CP models can be recovered exactly (without the usual whitening step). It would be interesting to analyze our algorithm under non-orthogonal CP model. However, we would like to point here that even with fully observed tensor, exact factorization is known only for orthonormal tensors. Now, given that our method guarantees not only completion but also factorization of the tensor (which is essential for large scale applications), it is natural that our method would require a similar condition.

Optimal dependence on $r$. The numerical results suggests threshold sample size scaling as $\sqrt{r}$. This is surprising since the degrees of freedom in describing CP model scales as linearly in $r$. This implies that the $\sqrt{r}$ scaling can only hold for $r=O(\sqrt{n})$. In comparison, for matrix completion we know that the threshold scales as $r$. It would be important to understand why this change in dependence in $r$ happens for higher order tensors, and identify how it depends on $k$ for $k$-th order tensor completion.

References

Appendix A Proof of Theorem 2.1 for Initialization Analysis

We prove the following bound on the spectrum of random tensors:

Here we prove the theorem for general case where $T$ is not symmetric and might even have different dimensions $n_{1}$, $n_{2}$ and $n_{3}$. Inspired by , our strategy is as follows:

Reduce to $x$,$y$, and $z$ which belongs to discretized sets $\widetilde{S}_{n_{1}}$, $\widetilde{S}_{n_{2}}$, and $\widetilde{S}_{n_{3}}$;

Bound the contribution of light triples using concentration of measure;

Bound the contribution of heavy triples using the discrepancy property of a random tripartite hypergraph.

Define a discretization of an $n$-dimensional ball as

It is therefore enough show that the bound holds for discretized vectors all discretized vectors $x$, $y$, and $z$. One caveat is that such a probabilistic bound must hold with probability sufficiently close to one such that we can apply the union bound over all discretized choices of $x$, $y$, and $z$. The following lemma bounds the number of such choices.

The size of the discretized set is bounded by |\widetilde{S}_{n}|\leq\big{(}\Delta/10\big{)}^{n}.

A naive approach to upper bound $({\cal P}_{\Omega}(T)-p\,T)[x,y,z]$ would be to consider it as a random variable and apply concentration inequalities directly. However, this naive approach fails since $x$, $y$ and $z$ can contain entries that are much larger than their typical value of $O(1/\sqrt{n})$. We thus separate the analysis into two contributions, and apply concentration inequalities to bound the contribution of the light triples and use graph topology of the random sampling to bound the contribution of the heavy triples. Define the light triples as

Heavy triples are defined as its complement $\overline{\cal L}=\{[n_{1}]\times[n_{2}]\times[n_{3}]\}\setminus{\cal L}$. Later we will set the appropriate value for $\epsilon=\Theta(p\sqrt{n_{1}n_{2}n_{3}})$. We can then write each contributions separately as

We will prove that both contributions are upper bounded by $CT_{\rm max}(\log n)^{2}\sqrt{(n_{1}n_{2}n_{3})^{1/2}p}$ with some positive constant $C$ for all $x\in\widetilde{S}_{n_{1}}$, $y\in\widetilde{S}_{n_{2}}$, and $z\in\widetilde{S}_{n_{3}}$. The bound on the light triples follows from Chernoff’s concentration inequalities. The bound on the heavy triples follows from the discrepancy property of random hyper graphs, which implies that there cannot be too many triples with large contributions. Theorem 2.1 then follows from Lemma A.1 with an appropriate choice of $\Delta=\Theta(1)$.

Let Z\equiv\sum_{(i,j,k)\in{\cal L}}\big{(}{\cal P}_{\Omega}(T)_{ijk}x_{i}y_{j}z_{k}\big{)}-p\,T[x,y,z] for some $x\in S_{n_{1}}$, $y\in S_{n_{2}}$, and $z\in S_{n_{3}}$. We claim that

We first show that the mean of $Z$ is bounded as

We next show concentration of $Z$ around item mean. Let $\lambda=\sqrt{n_{1}n_{2}n_{3}}/(2T_{\rm max}\sqrt{\epsilon})$ such that $|\lambda T_{ijk}x_{i}y_{j}z_{k}|\leq 1/2$ for all $(i,j,k)\in{\cal L}$. Then, $e^{\lambda T_{ijk}x_{i}y_{j}z_{k}}-1\leq\lambda T_{ijk}x_{i}y_{j}z_{k}+2\lambda(T_{ijk}x_{i}y_{j}z_{k})^{2}$.

Since $p=\epsilon/\sqrt{n_{1}n_{2}n_{3}}$, this proves that the contribution of light triples is bounded by $C\,T_{\rm max}\sqrt{\epsilon}$ with high probability.

Note that for the range of $p=\epsilon/n^{2}$, the contribution of light couples is bounded by $C\,T_{\rm max}\sqrt{\epsilon/n}$. However, even in this regime of $p$, the contribution of heavy triples is still $\Omega(1)$, which dominates the light triples by a factor of $\sqrt{n}$. This is the reason for the choice of $p=\Theta({\epsilon/n^{1.5}})$.

A.2 Bounding the contribution of heavy triples

The contribution of heavy triples is bounded by

In the following, we will show that the right-hand side of the above inequality is upper bounded by

for some positive numerical constant $C>0$ with probability larger than $1-n^{-5}$.

We consider a hypergraph $G=([n_{1}]\times[n_{2}]\times[n_{3}],E)$ with undirected hyper edges, where each edge connects three nodes, each one from each set $[n_{1}]$, $[n_{2}]$, and $[n_{3}]$. Given a sampling of entries in a tensor, we let the edges in $G$ denote the positions of the entries that is sampled. The proof is a generalization of similar proof for matrices in and is based on two properties of the hypergraph $G$. Define the degree of a node as the number of edges connected to that particular node such that ${\rm deg}_{1}(i)\equiv|\{(i,j,k)\in E\}|$, and similarly define ${\rm deg}_{2}(j)$ and ${\rm deg}_{3}(k)$. Define the degree of two nodes as the number of edges connected to both of the nodes such that ${\rm deg}_{12}(i,j)\equiv|\{(i,j,k)\in E\}|$, and similarly define ${\rm deg}_{13}(i,k)$ and ${\rm deg}_{23}(j,k)$.

Bounded degree property. A hyper graph $G$ satisfies the bounded degree property if the degree are upper bounded as follows:

for some positive numerical constant $\xi_{0}>0$ (independent of $n_{1},n_{2},n_{3}$ and $p$) where $p=|E|/(n_{1}n_{2}n_{3})$.

Discrepancy property. A hyper graph $G$ satisfies the discrepancy property if for any subset of nodes $A_{1}\in[n_{1}]$, $A_{2}\in[n_{2}]$, and $A_{3}\in[n_{3}]$, at least one of the following is true:

for some positive numerical constants $\xi_{1},\xi_{2}>0$ (independent of $n_{1},n_{2},n_{3}$ and $p$). Here, $e(A_{1},A_{2},A_{3})$ denotes the number of edges between the three subsets $A_{1},A_{2}$ and $A_{3}$, and ${\bar{e}}(A_{1},A_{2},A_{3})\equiv p\,|A_{1}|\,|A_{2}|\,|A_{3}|$ denotes the average number of edges between the three subsets.

We first prove that if the sampling pattern is defined by a graph $G$ which satisfies both the bounded degree and discrepancy properties, then the contribution of heavy triples is $O(\sqrt{\epsilon})$. Notice that this is a deterministic statement, that holds for all graphs with the above properties. We then finish the proof by showing that the random sampling satisfies both the bounded degree and discrepancy properties with probability at least $1-n^{-5}$.

We partition the indices according to the value of corresponding vectors:

for $u\in\{1,\ldots,\lceil\log_{2}(\sqrt{n_{1}}/\Delta)\rceil+1\}$, $v\in\{1,\ldots,\lceil\log_{2}(\sqrt{n_{2}}/\Delta)\rceil+1\}$, and $w\in\{1,\ldots,\lceil\log_{2}(\sqrt{n_{3}}/\Delta)\rceil+1\}$. We denote the size of each set by $a^{(u)}_{i}\equiv|A^{(u)}_{i}|$. We use $e_{uvw}$ to denote the number of edges between three subsets $A_{1}^{(u)}$, $A_{2}^{(v)}$, and $A_{3}^{(w)}$, and we use ${\bar{e}}_{uvw}\equiv p\,a_{1}^{(u)}a_{2}^{(v)}a_{3}^{(w)}$ to denote the average number of edges. Notice that the above definition of $A_{1}^{(u)}$’s cover all non-zero values of the entries of $x$, since, with discretization, the smallest possible positive value is $\Delta/\sqrt{n_{1}}$. The same applies to the entries of $y$ and $z$.

Note that since $\sum_{u}a_{1}^{(u)}2^{2(u-1)}\Delta^{2}/n_{1}\leq\|x\|^{2}\leq 1$, we get that

The contributions from various combinations of $(u,v,w)$ utilize various subsets of our assumptions. We prove that in each case the contribution is $O(\sqrt{\epsilon}(\log n)^{2})$ as follows.

Case1. For $(u,v,w)$ satisfying the first discrepancy property (13) : $e_{uvw}\leq\xi_{1}{\bar{e}}_{uvw}$.

In this case, using (15) and the fact that $p=\epsilon/\sqrt{n_{1}n_{2}n_{3}}$,

where $n\equiv\max\{n_{1},n_{2},n_{3}\}$ and in the last inequality we used the fact that we are summing over heavy triples satisfying $\Delta^{3}2^{u+v+w}>8\sqrt{\epsilon}$, and $\sum_{(u,v,w):2^{u+v+w}\leq 8\sqrt{\epsilon}/\Delta^{3}}2^{-(u+v+w)}\leq 2\log_{2}(\sqrt{n_{1}}/\Delta)\,\log_{2}(\sqrt{n_{2}}/\Delta)\,\Delta^{3}/(8\sqrt{\epsilon})$.

Case2. For $(u,v,w)$ satisfying the second discrepancy property in (14).

Case 2-1. For $(u,v,w)$ satisfying $\ln(e_{uvw}/{\bar{e}}_{uvw})\leq(1/2)\ln(en_{3}/a_{3}^{(w)})=(1/4)(\ln(en_{3}/(a_{3}^{(w)}2^{2w}))+\ln(2^{2w}))$.

Case 2-1-1. When $\ln(2^{2w})\leq\ln(en_{3}/(a_{3}^{(w)}2^{2w}))$, we have $\ln(e_{uvw}/{\bar{e}}_{uvw})\leq\ln(en_{2}/(a_{3}^{(w)}2^{2w}))$, which gives

It follows that $\sum\sigma_{uvw}\leq(16/\Delta)p\sqrt{n_{1}n_{2}n_{3}}2^{-u-v-w}\leq 2\Delta^{2}\sqrt{\epsilon}(\log n)^{2}$ using the fact that we are summing over heavy triples.

Case 2-1-2. When $\ln(2^{2w})>\ln(en_{3}/(a_{3}^{(w)}2^{2w}))$, we have $\ln(e_{uvw}/{\bar{e}}_{uvw})\leq\ln(2^{w})$.

Case 2-1-2-1. For $\sqrt{\epsilon}e_{uvw}>2^{u+v+w}{\bar{e}}_{uvw}$, it follows that $2^{u+v}\leq\sqrt{\epsilon}$. Since we are in the case where the first discrepancy does not hold, i.e. $e_{uvw}>\xi_{1}{\bar{e}}_{uvw}$, and the the second discrepancy property holds, we have $e_{uvw}\leq e_{uvw}\ln(e_{uvw}/{\bar{e}}_{uvw})\leq\xi_{2}a_{3}^{(w)}\ln(en_{3}/a_{3}^{(w)})\leq 2\xi_{2}a_{3}^{(w)}\ln(2^{2w})$ Then,

which is O\big{(}\sqrt{\epsilon}(\log n)^{2}\sqrt{(n_{3}/(n_{1}n_{2}))}\,\big{)}

Case 2-1-2-2. For $\sqrt{\epsilon}e_{uvw}\leq 2^{u+v+w}{\bar{e}}_{uvw}$,

Case 2-2. For $(u,v,w)$ satisfying $\ln(e_{uvw}/{\bar{e}}_{uvw})>(1/2)\ln(en_{3}/a_{3}^{(w)})$.

Case 2-2-1. For $2^{u+v}\leq\sqrt{n_{1}n_{2}\epsilon/n_{3}}2^{w}$, we know from the condition $\ln(e_{uvw}/{\bar{e}}_{uvw})>(1/2)\ln(en_{3}/a_{3}^{(w)})$, that $e_{uvw}\leq 2\xi_{2}a_{3}^{(w)}$. Then,

which is $O(\sqrt{\epsilon}(\log n)^{2})$.

Case 2-2-2. For $2^{u+v}>\sqrt{n_{1}n_{2}\epsilon/n_{3}}2^{w}$

Case 2-2-2-1. For $(u,v,w)$ satisfying bounded degree property with ${\rm deg}_{12}(i,j)\leq\xi_{0}pn_{3}$, we have $e_{uvw}\leq a_{1}a_{2}\xi_{0}pn_{3}$. Then,

which is $O(\sqrt{\epsilon}\sqrt{n_{3}/(n_{1}n_{2})}(\log_{2}n)^{2})$.

Case 2-2-2-2. For $(u,v,w)$ satisfying bounded degree property with ${\rm deg}_{12}(i,j)\leq\xi_{0}\log n_{3}$, we have $e_{uvw}\leq a_{1}a_{2}\xi_{0}\log n$.

which is $O((1/\sqrt{\epsilon})(\log n)^{3})$.

For $\epsilon\geq\log n$, this proves that the contribution of the heavy triples is $O(\sqrt{\epsilon}(\log n)^{2})$.

We are left to prove that the bounded degree and the bounded discrepancy properties hold for a random tripartite hypergraph $G=(V_{1}\cup V_{2}\cup V_{3},E)$ where each edge is selected with probability $p$. Precisely, let $n=\max\{|V_{1}|,|V_{2}|,|V_{3}|\}$, then the following lemma provides a bound on the degree and discrepancy, with high probability.

For any $\delta\in[0,1/e]$ and $p\geq(1/n^{2})\log n$, there exists numerical constants $C,C^{\prime}>0$ such that a random tripartite hyper graph satisfies the bounded degree property: for all $i\in V_{1}$, $j\in V_{2}$, and $k\in V_{3}$,

and the bounded discrepancy property: for all subsets $A_{1}\subseteq V_{1}$, $A_{2}\subseteq V_{2}$, and $A_{3}\subseteq V_{3}$, at least one of the following is true.

where $n_{1}=|V_{1}|$, $n_{2}=|V_{2}|$, $n_{3}=|V_{3}|$, $n=\max\{n_{1},n_{2},n_{2}\}$ and $\alpha\equiv\max n_{i}/n_{j}$.

Now, for the choice of $\delta=n^{-5}$, the bounded degree and discrepancy properties in (12), (13), and (14) hold for random tripartite hypergraphs. This finishes the proof of Theorem 2.1.

A.3 Proof of the bounded degree and discrepancy properties in Lemma A.3

We first prove the bounded degree properties of (12) hold with probability at least $1-\delta$. Applying standard concentration inequality, e.g. Bernstein inequality, we get that for some positive constant $\delta>0$,

for $n$ sufficiently large, and taking union bound over all choices of $i$, $j$ and $k$, ${\rm deg}_{1}(i)$, ${\rm deg}_{2}(j)$, and ${\rm deg}_{3}(k)$’s are uniformly bounded with probability at least $1-\delta/2$.

Similarly, we can apply concentration inequality to bound for some positive constant $\delta>0$

Applying the union bound over all choices of $(i,j)$, $(i,k)$ and $(j,k)$, we get that the bound holds uniformly with probability at least $1-\delta/2$.

Next, we prove that the random hyper graphs satisfy the discrepancy properties of (13) and (14). For any given subsets $A_{1}\subseteq[n_{1}]$, $A_{2}\subseteq[n_{2}]$, and $A_{2}\subseteq[n_{3}]$, let $a_{1},a_{2}$, and $a_{3}$ denote the cardinality of the subsets, and ${\bar{e}}(A_{1},A_{2},A_{3})=pa_{1}a_{2}a_{3}$.

Let’s assume, without loss of generality, that $a_{1}\leq a_{2}\leq a_{3}$. We divide the analysis into two cases depending on the size of the smallest subset. When at least two of the subsets are large, i.e. $a_{2}=\Omega(n)$ and $a_{3}=\Omega(n)$, then by bounded degree property, we can prove that (13) holds. However, when $a_{1}$ and $a_{2}$ are small, e.g. $O(1)$, then the first discrepancy no longer holds, and we need a different technique to show concentration.

From the bounded degree property, we know that ${\rm deg}_{1}(i)\leq 2pn_{2}n_{3}+(8/3)\ln(3n_{1}/\delta)$. Then,

We use the following bound on sum of indicator variables deviating from the mean :

where we denote ${\bar{e}}(A_{1},A_{2},A_{3})$ by ${\bar{e}}$, which holds for $t\geq 4$. For the bound holds with probability at least $1-\delta$, we require

where the term $1/(n_{1}n_{2}n_{3})$ is chosen to compensate for the union bound over all choices of $a_{1}$, $a_{2}$ and $a_{3}$. Simplifying the combinatorial terms, we get

We assumed that $a_{1}\leq n_{1}/e$, and since $x\ln(n_{1}/x)$ is monotone in $x\in[1,n_{1}/e]$, we know that $a_{1}\ln(en_{1}/a_{1})\geq\ln n_{1}$.

To lighten the notations, let’s suppose $a_{1}\ln(e\,n_{1}/a_{1})\leq a_{2}\ln(e\,n_{2}/a_{2})\leq a_{3}\ln(e\,n_{3}/a_{3})$. Let $t^{\prime}$ be the smallest number such that (3/{\bar{e}})\big{(}6a_{3}\ln(e\,n_{3}/a_{3})+\ln(\alpha^{2}/\delta)\,\big{)}=t^{\prime}\ln t^{\prime}.

For the regime of parameters such that $t^{\prime}\leq 4$, then $e(A_{1},A_{2},A_{3})\leq 4{\bar{e}}(A_{1},A_{2},A_{3})$ with probability at least $1-\delta$ the bounded discrepancy condition, in particular the first one, holds.

For the regime of parameters such that $t^{\prime}>4$, we can apply (16) to get that with probability at least $1-\delta$, the following holds uniformly for all choices of $A_{1}$, $A_{2}$, and $A_{3}$:

Since we defined $t^{\prime}$ to satisfy ${\bar{e}}t^{\prime}\ln t^{\prime}=18a_{3}\ln(e\,n_{3}/a_{3})+3\ln(\alpha^{2}/\delta)$, we have

As $t^{\prime}$ upper bounds $e(A_{1},A_{2},A_{3})/{\bar{e}}$, we have

A.4 Proof of Thresholding

Appendix B Alternating Minimization Analysis

If ${\bf u}_{1}$ and ${\bf u}_{2}$ are $2\mu$-incoherent, then there exists a positive constant $C$ such that for $p\geq C\frac{\mu^{3}\log^{2}n}{n^{1.5}}$ the following holds (w.p. $\geq 1-\log(1/\epsilon)/n^{8}$):

where $d_{\infty}([{\bf u}_{1},{\bf u}_{2}],[{\bf u}^{*}_{1},{\bf u}^{*}_{2}])=\max_{i,1\leq i\leq 2}\|{\bf u}_{i}-{\bf u}^{*}_{i}\|_{2}.$ Moreover, ${\bf u}_{1}^{t+1}$, ${\bf u}_{2}^{t+1}$ are both $2\mu$-incoherent.

We claim that with probability at least $1-1/n^{8}$,

for both $i\in\{1,2\}$. This proves the desired bound. Incoherence of $[{\bf u}^{t+1}_{1},{\bf u}^{t+1}_{2}]$ follows from Lemma B.2. Without loss of generality, we only prove the claim for $i=1$. Recall that $\widehat{{\bf u}}^{t+1}_{1}$ is the solution of the least squares problem in Step 11 of Algorithm 1, and can be written as

Note that the update that can be written in a vector form:

where $B$, $C$, $F,G$ are all diagonal matrices, s.t.,

Let $\widehat{{\bf u}}_{1}^{t+1}-\langle{\bf u}_{1},{\bf u}^{*}_{1}\rangle^{2}{\bf u}^{*}_{1}={\rm err}^{0}+{\rm err}^{1}+{\rm err}^{2}$, such that

We separate the analysis for each of the error terms. Using Lemma B.3, we have:

Setting $p\geq C\frac{\mu^{3}\log^{2}n}{\gamma^{2}n^{3/2}}$ for a $\gamma$ to be chosen appropriately later and using Lemma B.7 and Lemma B.5, we have (w.p. $\geq 1-2/n^{9}$):

Similarly, using Lemma B.4 and $p\geq C\frac{\mu^{3}\log^{2}n}{\gamma^{2}n^{3/2}}$, we have (w.p. $\geq 1-1/n^{9}$):

Since $1-\langle{\bf u}_{1},{\bf u}^{*}_{1}\rangle^{2}=(1/2)\|{\bf u}_{1}-{\bf u}^{*}_{1}\|_{2}^{2}$, we have from (21), (22), and (23) that (w.p. $\geq 1-10/n^{9}$):

Setting $\gamma\leq 1/200$ and for $d_{\infty}\left([{\bf u}_{1},{\bf u}_{2}],[{\bf u}^{*}_{1},{\bf u}^{*}_{2}]\right)\leq{1}/{200}$ as per our assumption, this proves the desired bound. $\Box$

B.2 Technical lemmas for rank-two analysis

The next lemma shows that all our estimates are $2\mu$-incoherent, which in turn allows us to bound the error in the above proof effectively. Note that the incoherence of the updates do not increase beyond a global constant ($2\mu$). Let $\widehat{{\bf u}}_{1}^{t+1}$ be obtained by update (17) and let ${\bf u}_{1}^{t+1}=\widehat{{\bf u}}_{1}^{t+1}/\|\widehat{{\bf u}}_{1}^{t+1}\|_{2}$.

Under the hypotheses of Theorem B.1, ${\bf u}_{1}^{t+1}$ is $2\mu$-incoherent with probability at least $1-1/n^{9}$.

Using (17) and the definitions of $B$, $C$, $F$, $G$ given in (19), we have:

where the second inequality follows by bounds on $B_{ii},C_{ii},F_{ii},G_{ii}$ obtained using Lemma B.5 and the distance bound $d_{\infty}\left([{\bf u}_{1},{\bf u}_{2}],[{\bf u}^{*}_{1},{\bf u}^{*}_{2}]\right)$. $\Box$

Next, we bound the first error term in (18).

Let ${\bf u}={\bf u}^{*}+{\bf d}^{u}$ and ${\bf v}={\bf v}^{*}+{\bf d}^{v}$, where ${\bf u},{\bf u}^{*},{\bf v},{\bf v}^{*}$ are all unit vectors and ${\bf u}^{*}\perp{\bf v}^{*}$. Also, let $\|{\bf d}^{u}\|_{2}\leq 1$ and $\|{\bf d}^{v}\|_{2}\leq 1$. Then, the following holds:

Now, $\langle{\bf u},{\bf v}^{*}\rangle=\langle{\bf d}^{u},{\bf v}^{*}\rangle\leq\|{\bf d}^{u}\|_{2}$. Also, $\langle{\bf u},{\bf v}\rangle^{2}\leq\langle{\bf u},{\bf v}\rangle=(\langle{\bf d}^{u},{\bf v}^{*}\rangle+\langle{\bf u}^{*},{\bf d}^{v}\rangle+\langle{\bf d}^{u},{\bf d}^{v}\rangle)\leq 2(\|{\bf d}^{u}\|+\|{\bf d}^{v}\|)$. Lemma now follows by combining the above observations with (26). $\Box$

Now, we bound the third error term in (18). Note that although the two individual terms ($(\langle{\bf u}_{1},{\bf u}^{*}_{2}\rangle^{2}B-F){\bf u}^{*}_{2}$ and $(\langle{\bf u}_{1},{\bf u}_{2}\rangle^{2}B-G){\bf v}$) are both small, still it is critical to bound the difference as the individual terms can be as large as a constant, even when ${\bf u}_{1}={\bf u}^{*}_{1}$ and ${\bf u}_{2}={\bf u}^{*}_{2}$. However, the difference goes down linearly with $\|{\bf u}_{2}-{\bf u}^{*}_{2}\|_{2}$.

Let $B,C,F,G$ be defined as in (19). Also, let the assumptions of Theorem B.1 hold. Also, let $p\geq C\frac{\mu^{3}\log^{2}n}{\gamma^{2}n^{3/2}}$, where $C>0$ is a global constant. Then, the following holds with probability $\geq 1-4/n^{9}$:

Let ${\bf u}_{2}={\bf u}^{*}_{2}+{\bf d}^{u}_{2}$ and ${\bf u}_{1}={\bf u}^{*}_{1}+{\bf d}^{u}_{1}$. Then,

Combining the above equation with (27), we get:

Lemma now follows using Lemma B.7, B.8, and the above equation. $\Box$

We now present a few technical lemmas that are critical to our proofs of the above given lemmas.

Then, the following holds with probability $\geq 1-1/n^{10}$:

where $\gamma\leq C/\log n$, where $C>0$ is a global constant.

Then, the following holds with probability $\geq 1-1/n^{10}$:

where $\gamma\leq C/\log n$, where $C>0$ is a global constant.

where $B,R$ are both diagonal matrices with $B(i,i)=\frac{1}{p}\sum_{jk}\delta_{ijk}({\bf u}(j))^{2}({\bf u}(k))^{2}$, and $R(i,i)=\frac{1}{p}\sum_{jk}\delta_{ijk}{\bf u}(j){\bf u}(k)a(j)b(k)$.

where $B,R$ are diagonal matrices s.t. $B(i,i)=\frac{1}{p}\sum_{jk}\delta_{ijk}({\bf u}(j))^{2}({\bf u}(k))^{2}$, $R(i,i)=\frac{1}{p}\sum_{jk}\delta_{ijk}{\bf u}(j){\bf u}(k)a(j)b(k)$.

B.3 Proofs of Technical Lemmas

Let $X_{jk}=\frac{1}{p}\delta_{jk}{\bf u}(j){\bf u}^{*}(j){\bf u}(k){\bf u}^{*}(k).$ Note that, $|X_{jk}|\leq\frac{\mu^{4}}{pn^{2}}$. Also,

Hence, using Bernstein’s inequality, we have:

Lemma now follows by selecting $t=C/\log n$. $\Box$

Let $X_{jk}=\frac{1}{p}\delta_{jk}{\bf u}(j)a(j){\bf u}(k)b(k).$ Note that, $|X_{jk}|\leq\frac{\mu^{3}\|b\|_{2}}{pn^{1.5}}$. Also,

Hence, using Bernstein’s inequality, we have:

Lemma now follows by selecting $t=\gamma\|b\|_{2}$. $\Box$

where $Z_{ijk}=\frac{1}{p}\delta_{ijk}c_{i}(\langle{\bf u},a\rangle\langle{\bf u},b\rangle(u(j))^{2}(u(k))^{2}-u(j)u(k)a(j)b(k)){\bf e}_{i}$. Note that,

as $p\geq\frac{C\mu^{3}(\log^{2}n)}{\gamma\cdot n^{3/2}}$. Also,

Hence, for $p$ and $\gamma$ mentioned above, we have:

Lemma now follows by using Bernstein’s inequality and the fact that $\sum_{ijk}Z_{ijk}=0$. $\Box$

Consider the $i$-th element of the diagonal matrix $(\langle{\bf u},a\rangle\langle{\bf u},b\rangle B-R)=\langle{\bf u},a\rangle\langle{\bf u},b\rangle B(i,i)-R(i,i)$. Now, using Lemma B.5, $|B(i,i)|\leq 1+\gamma$ w.p. $\geq 1-1/n^{10}$. Similarly, using Lemma B.6, $|R(i,i)-\langle{\bf u},a\rangle\langle{\bf u},b\rangle|\leq\gamma\|b\|_{2}$. Hence, w.p. $\geq 1-1/n^{10}$, we have:

Lemma now follows by observing that $\|(\langle{\bf u},a\rangle\langle{\bf u},b\rangle B-R)\|_{2}=\max_{i}|\langle{\bf u},a\rangle\langle{\bf u},b\rangle B(i,i)-R(i,i)|$ and using the above mentioned bound with union bound. $\Box$

B.4 Proof of Theorem 2.3 and general rank-r𝑟r analysis of alternating minimization

We prove the theorem by showing the following for all $q$:

The update for $\widehat{{\bf u}_{q}}^{t+1}$ is given by:

Using Lemma B.7, we have for all $p$ satisfying $p\geq(C\mu^{3}(\log n)^{2})/(\gamma^{2}\,n^{3/2})$, with probability at least $1-2/n^{10}$: