Fast and robust tensor decomposition with applications to dictionary learning

Introduction

Tensor decomposition is studied extensively across many disciplines including machine learning and signal processing. It is a powerful primitive for solving a wide range of other inverse / learning problems, for example: blind source separation / independent component analysis [LCC07], learning phylogenetic trees and hidden Markov models [MR05], mixtures of Gaussians [HK13], topic models [AFH+12], dictionary learning [BKS15, MSS16], and noisy-or Bayes nets [AGMR16].

A classical algorithm based on simultaneous diagonalization [Har70, LMV96] (often attributed to R. Jennrich) can decompose the input tensor (1.1) when the components are linearly independent, there is no error, and the order of the tensor is at least $3$. Current research on algorithms for tensor decomposition aims to improve over the guarantees of this classical algorithm in two important ways:

What conditions allow us to decompose tensors when components are linearly dependent?

What kind of errors can efficient decomposition algorithms tolerate? Can we tolerate errors $E$ with “magnitude” comparable to the low-rank part $\sum_{i=1}^{n}a_{i}^{\otimes k}$?

The focus of this work is on robustness. There are two ways in which errors arise in applications of tensor decomposition. The first is due to finite samples. For example, in some applications $T$ is the empirical $k$-th moment of some distribution and the error $E$ accounts for the difference between the empirical moment and actual moment (“population moment”). Errors of this kind can be made smaller at the expense of requiring a larger number of samples from the distribution. Therefore, robustness of decomposition algorithms helps with reducing sample complexity.

Another way in which errors arise is from modeling errors (“systematic errors”). These kinds of errors are more severe because they cannot be reduced by taking larger samples. Two important applications of tensor decomposition with such errors are learning Noisy-or Bayes networks [AGMR16] and sparsely-used dictionaries [BKS15]. For noisy-or networks, the errors arise due to non-linearities in the model. For sparsely-used dictionaries, the errors arise due to unknown correlations in the distribution of sparse feature representations. These examples show that robust tensor decomposition allows us to capture a wider range of models.

Robustness guarantees for tensor decomposition algorithms have been studied extensively (e.g., the work on tensor power iteration [AGH+14]). The polynomial-time algorithm with the best known robustness guarantees for tensor decomposition [BKS15, MSS16] are based on the sum-of-squares (SOS) method, a powerful meta-algorithm for polynomial optimization problems based on semidefinite programming relaxations. Unfortunately, these algorithms are far from practical and have polynomial running times with large exponents. The goal of this work is to develop practical tensor decomposition algorithms with robustness guarantees close to those of SOS-based algorithms.

In this work, we develop an easy-to-implement, randomized spectral algorithm to decompose 4-tensors with orthonormal components even when the error tensor $E$ has small but constant spectral norm as a $d^{2}$-by-$d^{2}$ matrix. This robustness guarantee is qualitatively optimal with respect to this norm in the sense that an error tensor $E$ with constant spectral norm (as a $d^{2}$-by-$d^{2}$ matrix) could change each component by a constant proportion of its norm. To the best of our knowledge, the only previous algorithms with this kind of robustness guarantee are based on SOS.We remark that the aforementioned analysis [MSS16] of Jennrich’s algorithm can tolerate errors $E$ if its spectral norm as a non-square $d^{3}$-by-$d$ matrix is constant. However, this norm of $E$ can be larger by a $\sqrt{d}$ factor than its spectral norm as a square matrix. Our algorithm runs in time $d^{2+\omega}\leqslant d^{4.373}$ using fast matrix multiplication. Even without fast matrix multiplication, our running time of $d^{5}$ is close to linear in the size of the input $d^{4}$ and significantly faster than the running time of SOS. As we will discuss later, an extension of this algorithm allows us to solve instances of dictionary learning that previously could provably be solved only by SOS.

A related previous work [HSSS16] also studied the question how to achieve similar guarantees for tensor decomposition as SOS using just spectral algorithms. Our algorithms follow the same general strategy as the algorithms in this prior work: In an algorithm based on the SOS semidefinite program, one solves a convex programming relaxation to obtain a “proxy” for a true, integral solution (such as a component of the tensor). Because the program is a relaxation, one must then process or “round” the relaxation or proxy into a true solution. In our algorithms, as in [HSSS16], instead of finding solutions to SOS semidefinite programs, the algorithms find “proxy objects” that behave in similar ways with respect to the rounding procedures used by SOS-based algorithms. Since these rounding procedures tend to be quite simple, there is hope that generating proxy objects that “fool” these procedures is computationally more efficient than solving general semidefinite programs.

However, our algorithmic techniques for finding these “proxy objects” differ significantly from those in prior work. The reason is that many of the techniques in [HSSS16], e.g., concentration inequalities for matrix-valued polynomials, are tailored to average-case problems and therefore do not apply in our setting because we do not make distributional assumption about the errors $E$.

Dictionary learning, also known as sparse coding, is studied extensively in neuroscience [OF97], machine learning [EP07, MRBL07], and computer vision [EA06, MLB+08, YWHM08]. Most algorithms for this problem used in practice do not come with strong provable guarantees. In recent years, several algorithms with provable guarantees have been developed for this problem [AAJ+14, ABGM14, AGMM15, BKS15, MSS16, HM16].

For the case that the coordinates of the distribution $\{x\}$ are independent (and non-Gaussian) there is a well-known reductionThe reduction only requires $k$-wise independence where $k$ is the order to tensor the reduction produces. of this problem to tensor decomposition where the components are the columns of $A$ and the error $E$ can be made inverse polynomially small by taking sufficiently many samples. (In the case of independent coordinates, dictionary learning becomes a special case of independent component analysis / blind source separation, where this reduction originated.) Variants of Jennrich’s spectral tensor decomposition algorithm (e.g., [BCMV14, MSS16]) imply strong provable guarantees for dictionary learning in the case that $\{x\}$ has independent coordinates even in the “overcomplete” regime when $n\gg d$ (using a polynomial number of samples).

More challenging is the case that $\{x\}$ has non-independent coordinates, especially if those correlations are unknown. We consider a model proposed in [BKS15] (similar to a model in [ABGM14]): We say that the distribution $\{x\}$ is $\tau$-nice if

An extension of our aforementioned algorithm for orthogonal tensor decomposition with spectral norm error allows us to learn orthonormal dictionaries from $\tau$-nice distributions. To the best of our knowledge, the only previous algorithms to provably solve this problem use sum-of-squares relaxations [BKS15, MSS16], which have large polynomial running time. Our algorithm recovers a $0.99$ fraction of the columns of $A$ up to error $\tau$ from ${\widetilde{O}}(n^{3})$ samples, and runs in time $n^{3+O(\tau)}d^{4}$. By a standard reduction, our algorithm also works for non-singular dictionaries and the running time increases by a factor polynomial in the condition number of $A$.

1 Results

For tensor decomposition, we give an algorithm with close to linear running time that recovers the rank-1 components of a tensor with orthonormal components, so long as the spectral norm of the square unfoldings of the error tensor is small.

Furthermore, if $\varepsilon\leqslant O(\log\log n/\log^{3}n)$, the algorithm can recover all components up to error $O(\varepsilon)$ in time ${\widetilde{O}}(d^{2+\omega}+nd^{4})$ with high probability. Here, $\omega\leqslant 2.373$ is the matrix multiplication exponent.

The orthogonality condition may at first seem restrictive, but for most applications it is possible to take a tensor with linearly independent components and transform it to a tensor with orthogonal components, as we will do for our dictionary learning result below. Furthermore, our algorithm works as-is if the components are sufficiently close to orthonormal:

These robustness guarantees are comparable to those of the sum-of-squares-based algorithms in [BKS15, MSS16] for the undercomplete case, which are the best known. Meanwhile, the sum-of-squares based algorithms require solving large semidefinite programs, while the running time of our algorithms is close to linear in the size of the input, and our algorithms are composed of simple matrix-vector multiplications.

On the other hand, our algorithms fail to work in the overcomplete case, when the rank grows above $n$, and the components are no longer linearly independent. One interesting open question is whether the techniques used in this paper can be extended to the overcomplete case.

Dictionary learning

Using our tensor decomposition algorithm as a primitive, we give an algorithm for dictionary learning when the sample distribution is $\tau$-nice.

The total runtime is thus ${\widetilde{O}}(n^{3}d^{4})$—in the theorem statement, we write it in terms of the number of samples $m$ in order to separate the time spent processing the samples from the learning phase. We note that the sample complexity bound that we have, $m={\widetilde{O}}(n^{3})$, may very well be sub-optimal; we suspect that $m={\widetilde{O}}(n^{2})$ is closer to the truth, which would yield a better runtime.

We are also able to apply standard whitening operations (as in e.g. [AGHK13]) to extend our algorithm to dictionaries with linearly-independent, but non-orthonormal, columns, at the cost of polynomially many additional samples.

To our knowledge, our algorithms are the only remotely efficient dictionary learning algorithms with provable guarantees that permit $\tau$-nice distributions in which the coordinates of $x$ may be correlated by constant factors, the only other ones being the sum-of-squares semidefinite programming based algorithms of [BKS15, MSS16].

Preliminaries

We will also make frequent use of the following lemma, which states that the distance between two points cannot increase when both are projected onto a closed, convex set.

This lemma is well-known (see e.g. [Roc76]), but we will prove it for completeness in Appendix A.

Techniques

In this section we give a high-level overview of the algorithms in our paper, and of their analyses. We begin with the tensor decomposition algorithm, after which we’ll explain the (non-trivial) extension to the dictionary learning application. At the very end, we will discuss the relationship between our algorithms and sum-of-squares relaxations.

where $T_{ij}$ is the $i,j$th $d\times d$ matrix slice of the tensor $\mathbf{T}$. Since the $a_{i}$ are orthogonal, the coefficients $\langle g,a_{i}^{\otimes 2}\rangle$ are independent, and so we find ourselves in an ideal situation—$M_{g}$ is a sum of the orthogonal components we want to recover with independent Gaussian coefficients. A simple eigendecomposition will recover all of the $a_{i}$.

On the other hand, when we have a nonzero noise tensor $\mathbf{E}$, a random contraction along modes $\{1,2\}$ results in the matrix

where the $E_{ij}$ are $d\times d$ slices of the tensor $\mathbf{E}$. The last term, composed of the error, complicates things. Standard facts about Gaussian matrix series assert that the spectral norm of the error term behaves like $\|E_{\{1,2,3\}\{4\}}\|$, the spectral norm of a $d^{3}\times d$ reshaping of $E$, whereas we only have control over square reshapings such as $\|E_{\{1,2\}\{3,4\}}\|$.In fact, this observation was crucial in the analysis of [MSS16]—in that work, semidefinite programming constraints are used to control the spectral norm of the rectangular reshapings. These can be off by polynomial factors. If the Frobenius norm of $\mathbf{E}$ is $\|\mathbf{E}\|_{F}^{2}\approx\varepsilon^{2}d^{2}$, which is the magnitude one would expect from a tensor whose square reshapings are full-rank matrices with spectral norm $\varepsilon$, then we have that necessarily

since there are at most $d$ nonzero singular values of rectangular reshapings of $\mathbf{E}$. In this case, unless $\varepsilon\ll 1/\sqrt{d}$, the components $a_{i}a_{i}^{\top}$ are completely drowned out by the contribution of the noise, and so the robustness guarantees leave something to be desired.

Basic idea

The above suggests that, as long as we allow the error $\mathbf{E}$ to have large Frobenius norm, an approach based on random contraction will not succeed. Our basic idea is to take $\mathbf{T}$, whose error has small spectral norm, and transform it into a tensor $\mathbf{T}^{\prime}$ whose error has small Frobenius norm.

Because we do not know the decomposition of $\mathbf{T}$, we cannot access the error $\mathbf{E}$ directly. However, we do know that for any $d^{2}\times d^{2}$ reshaping $T$ of $\mathbf{T}$,

where $S=\sum_{i=1}^{n}a_{i}^{\otimes 4}$, and $\|E\|\leqslant\varepsilon$. The rank of $S$ is $n\leqslant d$, and all eigenvalues of $S$ are $1$. Thus, if we perform the operation

where $(\cdot)_{+}$ denotes projection to the cone of positive semidefinite matrices, we expect that the signal term $S$ will survive, while the noise term $E$ will be dampened. More formally, we know that $T$ has $n$ eigenvalues of magnitude $1\pm\varepsilon$, and $d^{2}-n$ eigenvalues of magnitude at most $\varepsilon$, and therefore $\operatorname{rank}(T^{>\varepsilon})\leqslant n$. Also by definition, $\|T-T^{>\varepsilon}\|\leqslant\varepsilon$. Therefore, we have that

with $\|E^{\prime}\|\leqslant\varepsilon$, and thus

Since $S,T^{>\varepsilon}$ are both of rank at most $n$, and $E^{\prime},E$ have spectral norm bounded by $\varepsilon$, we have that

Finally, to eliminate these large singular values, we will project $\mathbf{T}^{>\varepsilon}$ into the set of matrices whose rectangular reshapings have singular values at most $1$—because $S$ is a member of this convex set, the projection can only decrease the Frobenius norm. After this, we will apply the random contraction algorithm, as originally suggested.

Variance of Gaussian matrix series

Recall that we wanted to sample a random $d^{2}$-dimensional Gaussian vector $g$, and the perform the contraction

The error term on the right is a matrix Gaussian series. The following lemma describes the behavior of the spectra of matrix Gaussian series:

It turns out that for us it suffices to perform two projections in sequence—we first reshape $\mathbf{T}^{>\varepsilon}$ to the matrix $\mathbf{T}^{>\varepsilon}_{\{123\}\{4\}}$,project it to the set of matrices with singular values at most $1$, and then reshape the result along modes $\{124\}\{3\}$, and project to the same set again. As mentioned before, because projection to a convex set containing $S$ cannot increase the distance from $S$, the Frobenius norm of the new error can only decrease. What is less obvious is that performing the second projection will not destroy the property that the reshaping along modes $\{123\}\{4\}$ has spectral norm at most $1$. By showing that each projection corresponds to either left- or right- multiplication of $\mathbf{T}^{>\varepsilon}_{\{123\}\{4\}}$ and $\mathbf{T}^{>\varepsilon}_{\{124\}\{3\}}$ by matrices of spectral norm at most $1$, we are able to show that the second projection does not create large singular values for the first flattening, and so two projections are indeed enough. Call the resulting tensor $(\mathbf{T}^{>\varepsilon})^{\leqslant 1}$.

Now, if we perform a random contraction in the modes $1,2$, we will have

where the spectral norm of $\|\hat{E}_{g}\|\leqslant\sqrt{\log n}$ with good probability. So, ignoring for the moment dependencies between $\langle a_{i}^{\otimes 2},g\rangle$ and $E_{g}$, $\max_{i}|\langle g,a_{i}^{\otimes 2}\rangle|>1.1\cdot\|E_{g}\|$ with probability at least $n^{-.2}$, which will give $a_{i}$ correlation $0.9$ with the top eigenvector of $M_{g}$ with good probability.

Improving accuracy of components

The algorithm described thus far will recover components $b_{i}$ that are $0.9$-correlated with the $a_{i}$, in the sense that $\langle b_{i},a_{i}\rangle^{2}\geqslant 0.9$. To boost the accuracy of the recovered components, we use a simple method which resembles a single step of tensor power iteration.

We’ll use the closeness of our original tensor $\mathbf{T}$ to $\sum_{i}a_{i}^{\otimes 4}$ in spectral norm. We let $T$ be a $d^{2}\times d^{2}$ flattening of $\mathbf{T}$, and compute the vector

Now, when the vector $v$ is reshaped to a $d\times d$ matrix $V$, the term $E(b_{i}\otimes b_{i})$ is a matrix of Frobenius norm (and thus spectral norm) at most $\varepsilon$. By the orthonormality of the $a_{j}$, the sum of the coefficients in the second term is at most 0.1, and so $a_{i}$ is $\varepsilon$-close to the top eigenvector of $V$.

Recovering every component

Because the Frobenius norm of the error in $(\mathbf{T}^{>\varepsilon})^{\leqslant 1}$ is $O(\varepsilon\sqrt{n})$, there may be a small fraction of the components $a_{i}^{\otimes 4}$ that are “canceled out” by the error—for instance, we can imagine that the error term is $\hat{\mathbf{E}}=-\sum_{i=1}^{\varepsilon^{2}n}a_{i}^{\otimes 4}$. So while only a constant fraction of the components $a_{i}^{\otimes 4}$ can be more than $\varepsilon$-correlated with the error, we may still be unable to recover some fixed $\varepsilon^{2}$-fraction of the $a_{i}$ via random contractions.

To recover all components, we must subtract the components that we have found already and run the algorithm iteratively—if we have found $m=0.99n$ components, then if we could perfectly subtract them from $\mathbf{T}$, we would end up with an even lower-rank signal tensor, and thus be able to make progress by truncating all but $0.1n$ eigenvalues in the first step.

The challenge is that we have recovered $b_{1},\ldots,b_{m}$ that are only $(1-\varepsilon)$-correlated with the $a_{i}$, and so naively subtracting $\mathbf{T}-\sum_{i=1}^{m}b_{i}^{\otimes 4}$ can result in a Frobenius and spectral norm error of magnitude $\varepsilon\sqrt{m}$—thus the total error is still proportional to $\sqrt{n}$ rather than $\sqrt{0.1n}$.

Dictionary learning

We can use our tensor decomposition algorithm to learn the dictionary $A$, as long as the columns of $A$ are linearly independent. For the sake of this overview, assume instead that the columns of $A$ are orthonormal. Then given samples $y^{(1)},\ldots,y^{(m)}$ for $m=\operatorname{poly}(n)$, we can compute the $4$th moment tensor to accuracy $\varepsilon$ in the spectral norm,

In the sum-of-squares relaxation, this issue is overcome easily: because of the symmetries required of the SDP solution matrix $X$, $\langle X,a_{i}a_{j}^{\top}\otimes a_{i}a_{j}^{\top}\rangle=\langle X,a_{i}a_{i}^{\top}\otimes a_{j}a_{j}^{\top}\rangle$, so by linearity this low-rank error term cannot influence the objective function any more than the $a_{i}a_{i}^{\top}\otimes a_{j}a_{j}^{\top}$ term.

This removes the spectrum in the direction of the “nice” error terms, corresponding to $a_{i}a_{i}^{\top}\otimes a_{j}a_{j}^{\top}$ and $a_{i}a_{j}^{\top}\otimes a_{j}a_{i}^{\top}$. The fact that the rest of the matrix is low-rank means that we can apply an analysis similar to the analysis in the first step of our tensor decomposition to argue that $\|M^{>\varepsilon}-\sum_{i}a_{i}^{\otimes 4}\|_{F}\leqslant\|M-\sum_{i}a_{i}^{\otimes 4}\|_{F}$.

Now, we re-shape $M^{>\varepsilon}$, so that if initially we had the flattening $\mathbf{M}\to M_{\{1,2\}\{3,4\}}$, we look at the flattening $M^{>\varepsilon}_{\{1,3\}\{2,4\}}$. In this flattening, the term $a_{i}a_{j}^{\top}\otimes a_{i}a_{j}^{\top}$ from $M^{>\varepsilon}$ is transformed to $a_{i}a_{i}^{\top}\otimes a_{j}a_{j}^{\top}$, and so the problematic error term from the original flattening has spectral norm $\varepsilon$ in this flattening! Applying the projection

eliminates the problematic term, and brings us again closer to the target matrix $\sum_{i}a_{i}^{\otimes 4}$. Thus, we end up with a tensor that is close to $\sum_{i}a_{i}^{\otimes 4}$ in Frobenius norm, and we can apply our tensor decomposition algorithm.

Connection to sum-of-squares algorithms

We take a moment to draw parallels between our algorithm and tensor decomposition algorithms in the sum-of-squares hierarchy. In noisy orthogonal tensor decomposition, we want to solve the non-convex program

The (somewhat simplified) corresponding sum-of-squares relaxation is the semidefinite program

To solve this semidefinite program, one should project $\mathbf{T}$ into the intersection of all of the convex feasible regions of the constraints. However, projecting to the intersection is an expensive operation in terms of runtime. Instead, we choose a subset of these constraints, and project $\mathbf{T}$ into the set of points satisfying each constraint sequentially, rather than simultaneously. These are not equivalent projection operations, but because we select our operations carefully, we are able to show that our $\mathbf{T}$ is close to the SDP optimum in a sense that is sufficient for successfully running Jennrich’s algorithm.

Decomposing orthogonal 444-tensors

Truncate to 0 all eigenvalues of that have magnitude less than $\varepsilon$:

where we have used $(M)_{+}$ to denote projection to the PSD cone.

To compute $T^{>\varepsilon}$, we can compute the top $n$ eigenvectors of $T$. Since $O(\varepsilon)\cdot\lambda_{n}\geqslant\lambda_{n+1}$, where $\lambda_{n},\lambda_{n+1}$ are the $n$th and $(n+1)$st eigenvalues, we can compute this in time ${\widetilde{O}}(\min\{nd^{4},d^{2+\omega}\})$ via subspace power iteration (see for example [HP14]). ∎

Project $\mathbf{T}$ to the set of tensors whose rectangular reshapings along modes $\{1,2,3\},\{4\}$ have spectral norm at most $1$ (obtaining a new tensor $\hat{\mathbf{T}}$ with $\|\hat{\mathbf{T}}_{\{1,2,3\}\{4\}}\|\leqslant 1$).

Project $\hat{\mathbf{T}}$ to the set of tensors whose rectangular reshapings along modes $\{1,2,4\},\{3\}$ have spectral norm at most $1$, obtaining a new tensor $\mathbf{T}^{\leqslant 1}=\mathbf{S}+\mathbf{E}^{\leqslant 1}$.

Output: $u_{L}$ and $u_{R}$, the top left- and right- unit singular vectors of $M_{g}$.

Under the appropriate conditions on $\mathbf{T}$, with probability ${\widetilde{O}}(n^{-\varepsilon})$, Algorithm 4.3 will output a vector that is $0.9$-correlated with $a_{j}^{\otimes 2}$ for some $j\in[n]$.

Recovering one component requires time ${\widetilde{O}}(d^{2+\omega}n^{O(\eta)})$, and recovering $m$ components requires time ${\widetilde{O}}(\max\{mn^{O(\eta)}d^{4},d^{2+\omega}\})$.

We can then post-process the vectors to obtain a vector that has correlation $1-\varepsilon$ with $a_{j}$—the details are given in Section 4.3 below.

We will now prove the correctness of Algorithm 4.3 step-by-step, tying details together at the end of this subsection. First, we argue that in step 1, the truncation of the large eigenvalues cannot increase the Frobenius norm of the error.

Suppose that we define $\mathbf{T}^{\leqslant 1}$ to be the result of projecting $\mathbf{T}=\mathbf{S}+\mathbf{E}$ to the set of tensors whose rectangular reshapings along modes $\{1,2,3\},\{4\}$ have spectral norm at most $1$, then projecting the result to the set of tensors whose rectangular reshapings along modes $\{1,2,3\},\{4\}$ have spectral norm at most $1$. Then $\|\mathbf{T}^{\leqslant 1}-\mathbf{S}\|_{F}\leqslant\|E\|_{F}$, and

This operation requires time ${\widetilde{O}}(d^{2+\omega})$.

To establish the first claim, we note that the tensor $\mathbf{T}^{\leqslant 1}$ was obtained by two projections of different rectangular reshapings of the matrix $S+E$ to the set of rectangular matrices with singular value at most $1$. This set is closed, convex, and contains $S$, and so the error can only decrease in Frobenius norm (see Lemma A.2 in Appendix A for a proof),

Finally, each reshaping step takes $O(d^{4})$ time. Since we are only interested in truncating large singular values, it suffices for us to compute the SVD corresponding to singular values between $\sqrt{n}$ and $1$. This can be done via subspace power iteration, which here involves the multiplication of a $d^{3}\times d$ matrix and a $d\times d$ matrix (with intermediate orthogonalization steps for the $d\times d$ matrix, see [HP14]), which requires time ${\widetilde{O}}(d^{2+\omega})$, where $\omega$ is the matrix multiplication constant. Going forward the representation of the matrix will be as the original matrix, with the subtracted SVD corresponding to large singular values. ∎

We will need to argue that if the Frobenius norm of the error matrix is small, this is a sufficient condition under which we succeed. For this, we will use the following two lemmas. The first tells us that with probability ${\widetilde{\Omega}}(n^{-O(\varepsilon)})$ over the choice of $g$, we will have for some $i\in[n]$ that

where $|c|\geqslant\|N\|$ and furthermore $\|Na_{i}\|$ and $\|N^{\top}a_{i}\|$ are small:

Then for a $1-3\delta$ fraction of $j\in[n]$,

In particular, if $\delta=\Omega(1)$, $\beta=O(\varepsilon)$, $\beta<1$, then this probability is ${\widetilde{\Omega}}(n^{-(1+O(\varepsilon))})$.

The proof consists primarily of the application of concentration inequalities, and we provide it below in Section 4.4.

The second lemma states that if $M_{g}$ indeed has the form above, the top singular vectors of $M_{g}$ must be close to the component $a_{i}$.

with $|c|\geqslant(1+\beta)\|N\|$ for $\beta>0$, and $\|Na_{1}\|,\|N^{\top}a_{1}\|\leqslant\varepsilon|c|$ so that the relationship $\frac{2\varepsilon(1+\beta)}{\beta}<0.01$ holds. Then letting $u$ be a top singular vector of $M_{g}$, it follows that

The proof requires some careful calculations, but is not complicated, and we will prove it below in Section 4.4.

Finally, we are ready to stitch these arguments together and prove that Algorithm 4.3 works.

After reshaping and truncating $\mathbf{T}$ in step $1$ of the algorithm, by Lemma 4.5 the matrix $T^{\leqslant 1}=\sum_{i}a_{i}^{\otimes 4}+E$ has the properties that

and also that still $\|E\|_{F}\leqslant\eta\sqrt{n}$.

We can now apply Lemma 4.6 with $\delta=\frac{1}{300}$ and $\beta=400\eta/\delta=O(\eta)$ to conclude that for at least a $0.99$-fraction of the $i\in[n]$, with probability at least ${\widetilde{O}}(n^{-1-O(\eta)})$, we will have

where $\|N\|\leqslant(1+\beta)c$ and $\|Na_{i}\|,\|N^{\top}a_{i}\|\leqslant 48\eta\cdot|c|$. Applying Lemma 4.7, we have that either the left- or right- top unit singular vector $u$ of $M_{g}$ has correlation at least

For runtime, by our arguments in Lemma 4.5 step 1 takes time ${\widetilde{O}}(d^{2+\omega})$. After this, with either representation of our matrix $T^{\leqslant 1}$ (whether we compute the full truncated SVD or have the original matrix minus the subtracted SVD), performing power iteration to find the top eigenvector with the flattening $M_{g}$ takes time ${\widetilde{O}}(d^{4})$, and finding a single component takes ${\widetilde{O}}(n^{-O(\eta)})$ samples of random contractions. Since we can reuse $T^{\leqslant 1}$ with new random contractions, the total runtime for recovering one component is ${\widetilde{O}}(d^{2+\omega}n^{-O(\eta)})$, and by the independence of the runs recovering $m$ components requires ${\widetilde{O}}(d^{2+\omega})+mn^{-O(\eta)}\cdot{\widetilde{O}}(d^{4})$ time. ∎

With the core of our algorithm in place, we now take care of the remaining technical issues: recovery precision, working with near-orthonormal vectors, and recovering the full set of component vectors.

Because the precision of recovery will be important in not amplifying the error, we begin with our precision-amplifying postprocessing algorithm.

Output: If for one of $v\in\{v_{L},v_{R}\}$, $(v^{\otimes 2})^{\top}Tv^{\otimes 2}\geqslant(1-3\varepsilon)^{2}-\varepsilon$, output $v$.

Suppose that $v$ is a unit vector with $\langle v,a_{i}\rangle^{2}\geqslant 0.99$, and $T=\sum_{i}a_{i}^{\otimes 4}+E$ for $\|E\|\leqslant\varepsilon$ and $a_{1},\ldots,a_{n}$ orthonormal. Then if we let $A$ be the reshaping of $T(v\otimes v)$ to a $d\times d$ matrix, and if we let $u_{L},u_{R}$ be the top left- and right- unit singular vectors of $M$, then

In other words, Algorithm 4.8 succeeds. Further, the time required is ${\widetilde{O}}(d^{4})$.

Therefore, defining $M$ to be the $n\times n$ reshaping of $T(v\otimes v)$ and defining $N$ to be the $n\times n$ reshaping of $E(v\otimes v)$,

where $\|N\|_{F}\leqslant\varepsilon$, since $\|E\|\leqslant\varepsilon$.

and that if we choose $\eta$ so that $1-\eta\geqslant 2\alpha\geqslant 1/50$,

Thus, $\|M-\eta a_{1}a_{1}^{\top}\|\leqslant\|M\|-\eta+2\varepsilon$, and we have by Fact A.3 (see Appendix A for a proof) that the top unit eigenvector $u$ of $M$ is such that

Choosing $\eta=49/50$, the result follows. ∎

2 Near-orthonormal components

We’ll now dispense with the discrepancy between the orthonormal and near-orthonormal cases.

3 Full Recovery

Now we give the full algorithm, which will remove the components we find in each step from the tensor without amplifying the spectral norm of the error too much.

Initialize the set of known components $K=\emptyset$, and the set of components under inspection $B=\emptyset$.

Initialize a working copy of $\mathbf{T}$, $\mathbf{T}^{(0)}_{work}$, and keep a clean copy of $\mathbf{T}$ called $\mathbf{T}_{clean}$.

Preprocess $\mathbf{T}^{(t)}_{work}$ with Algorithm 4.1, then run Algorithm 4.3 with $\mathbf{T}_{work}^{(t)}$ ${\widetilde{O}}(n)$ times, then postprocess with Algorithm 4.8 using $\mathbf{T}_{clean}$ and error parameter $\varepsilon$. If this produces an output vector $v$, add $v$ to $B$ (unless $K$ already contains a vector that is $1-\varepsilon$ correlated with $v$).

Given $\mathbf{T}=\sum_{i=1}^{n}a_{i}^{\otimes 4}+E$ where the $a_{i}$ are orthonormal and $\|E_{\{1,2\}\{3,4\}}\|\leqslant\varepsilon$, then if $\varepsilon<O(\eta^{2}/\log^{2}n)$, with probability $1-o(1)$, Algorithm 4.11 recovers orthonormal vectors $b_{1},\ldots,b_{n}$ so that there exists a permutation $\pi:[n]\to[n]$ such that for each $i\in[n]$,

Furthermore, this requires runtime ${\widetilde{O}}(n^{1+O(\eta)}d^{2+\omega})$.

First, we prove that if we have an orthonormal basis that approximates $a_{1},\ldots,a_{k}$, we can subtract it without introducing a large spectral norm error—this motivates and justifies steps 3(b)–3(e).

So it suffices for us to bound $\|U-V\|\cdot\|U+V\|$.

By the subadditivity of the norm, $\|U+V\|\leqslant\|U\|+\|V\|=2$. Meanwhile, the singular values of $U-V$ are the square roots of the eigenvalues of $\|(U-V)^{\top}(U-V)\|$, and so we bound

where the second line follows because $U$ and $V$ have orthonormal columns. Now, by assumption we know that

where for $i\neq j$, $E_{ij}=\langle b_{i},a_{j}\rangle^{2}$ and $E_{ii}=\langle b_{i},a_{i}\rangle^{2}-(1-\varepsilon)$, and by the orthonormality of the $a_{j}$,

So the $1$-norm of the rows of $E$ is at most $\varepsilon$. By the orthonormality of the $b_{j}$, the same holds for the $1$-norm of the columns, $\sum_{i}|E_{ij}|$. It follows that $\|E\|\leqslant\varepsilon$. Therefore,

where $\|\hat{E}\|\leqslant 2\varepsilon$. Returning to (4.1), we can conclude that $\|(U-V)\|\leqslant 2\sqrt{\varepsilon}$, and we have our result. ∎

Where we have used that the spectral norm and nuclear norm are dual. Therefore,

Finally, we are ready to prove that Algorithm 4.11 works.

We claim that in the $t$th iteration of step $3$, with high probability we have at most $0.45^{t}n$ components remaining to be found, and that $\mathbf{T}_{work}^{(t)}=\mathbf{T}+F^{(t)}$ where $\|F^{(t)}\|\leqslant 8t\sqrt{\varepsilon}$. For $t=0$, this is easily true.

Now assume this holds for $t$, and we will prove it for $t+1$. Since by assumption $t\sqrt{\varepsilon}\log n\leqslant 100\sqrt{\varepsilon}\log n\leqslant O(\eta)$, applying Lemma 4.2 and Theorem 4.4 to the running of preprocessing Algorithm 4.1 and the main step Algorithm 4.3 with $\mathbf{T}_{work}^{(t)}$ and Lemma 4.9 to the running of the postprocessing Algorithm 4.8 with $\mathbf{T}_{clean}$, in step $3(a)$ with high probability we will find $m\geqslant 0.9n_{t}$ vectors $b_{1},\ldots,b_{m}$ so that $\langle b_{i},a_{i}\rangle^{2}\geqslant 1-3\varepsilon$. Furthermore, this takes a total of ${\widetilde{O}}(mn^{O(\eta)}d^{2+\omega})$ time.

for a matrix $F$ with $\|F\|\leqslant 8\sqrt{\varepsilon}$. By induction, this implies that $\mathbf{T}_{work}^{(t+1)}=\mathbf{T}+F^{(t+1)}$ where $\|F^{(t+1)}\|\leqslant\|F\|+\|F^{(t)}\|\leqslant(t+1)8\sqrt{\varepsilon}$.

Taking a union bound over the high-probability success of Algorithm 4.3, we have that after $t=O(\log n)$ steps we have found all of the components. We have spent a total of ${\widetilde{O}}(n^{1+O(\eta)}d^{2+\omega})$ time in step 3(a). Finally, steps 3(b)-3(e) of Algorithm 4.11 require no more than ${\widetilde{O}}(d^{3})$ time, and since the entire loop runs ${\widetilde{O}}(1)$ times, we have our result. ∎

4 Supporting Lemmas

Now we circle back and prove the omitted supporting lemmas.

For convenience, fix $j:=1$. Let $g^{(1)}$ be the component of $g$ in the direction $a_{1}^{\otimes 2}$, and let $g^{(>1)}$ be the component of $g$ orthogonal to $a_{1}^{\otimes 2}$. Notice that $g^{(1)},g^{(>1)}$ are independent.

By the orthogonality of the $a_{i}$, our matrix $M_{g}$ can be written as

where $T_{j}$ is the $j$th matrix slice of $T$. For convenience, we can refer to the two sums on the right as

First, we get a lower bound on the probability that the coefficient of $a_{1}a_{1}^{\top}$ is large. Let $\mathcal{G}_{1}(\alpha)$ be the event that $|\langle g^{(1)},a_{1}^{\otimes 2}\rangle|=\|g^{(1)}\|\geqslant\sqrt{2\alpha\log n}$. By standard tail estimates on univariate Gaussians, we have that

Now, we bound $\|N\|$. Define the event $\mathcal{E}_{>1}(\rho)$ to be the even that

To bound $\|N\|$, it thus remains to understand the term

and then by an application of Markov’s inequality. ∎

Thus, combining the above we have a two-part upper bound on $\|N\|$.

Finally, define the event $\mathcal{E}_{a_{1},E}(\theta)$ to be the event that

We can apply the same arguments to $\left(\sum_{j}g_{j}^{(>1)}E_{j}\right)^{\top}a_{1}$, and we conclude that

Now, by the union bound $\mathcal{E}_{>1}(\rho)$ and $\mathcal{E}_{a_{1},E}$ both occur with probability at least $1-d^{-\rho}-2d^{-\theta}$. Also, we notice that $\mathcal{E}_{>1}\cup\mathcal{E}_{a_{1},E}$ and $\mathcal{G}_{1}$ are independent. Therefore, for $\rho,\theta\geqslant\log\log n/\log n$,

Conditioning on $\mathcal{E}_{>1}$ and $\mathcal{G}_{1}$,

where $|c|\geqslant\sqrt{2\alpha\log n}$, and $N$ is a matrix of norm at most $(\varepsilon/\delta)c+\sqrt{2(1+\rho)\log d}$, such that $\|Na_{1}\|,\|N^{\top}a_{1}\|\leqslant(\varepsilon/\delta)(c+\sqrt{2(1+\theta)})$.

We now set $\alpha$ so that $|c|\geqslant\beta\|N\|$. This occurs when

Choosing $\beta=1+\beta^{\prime}$, $\rho,\theta=\log\log n/\log n$, we have our conclusion for $a_{1}$, and by symmetry for all other $a_{i}$ in the $1-3\delta$ fraction of $i\in[n]$ for which the Frobenius norms of the contractions are small. ∎

Now expanding the $v^{\top}N^{\top}Nv$ term along the components of $v$,

It is easy to see that when $c/(1+\beta)<c-\kappa$, this quantity is maximized at $\alpha=1$, and so by our choice of $\kappa$ we have that

and thus $\|M_{g}-\kappa a_{1}a_{1}^{\top}\|\leqslant(1+\varepsilon)c-\kappa$.

Where we have applied the Cauchy-Schwarz inequality, and the assumption that $\|Na_{1}\|\leqslant\varepsilon c$. It follows that

Finally applying Fact A.3, we can conclude that for either the left- or right-singular unit vector $u$ of $M_{g}$,

Choosing $\delta=\frac{2\varepsilon(1+\beta)}{\beta}$ as small as possible, we have our result. ∎

Learning Orthonormal Dictionaries

Here, we show how to use our tensor decomposition algorithm to learn dictionaries with orthonormal basis vectors.

Below, in Section 5.3, we will prove that ${\widetilde{O}}(n^{3})$ samples suffice to estimate the $4$th moment tensor within $o(1)$ spectral norm error. Computing this matrix from ${\widetilde{O}}(n^{3})$ samples takes ${\widetilde{O}}(n^{3}d^{4})$ time. Thus we can equivalently formulate the problem as follows:

Finally, given access to a sufficiently large number of samples, we can reduce to the case where the columns of $A$ are orthogonal:

The proof is straightforward, involving a transformation by the empirical covariance matrix, and we give it below in Section 5.2.

Reshape $\mathbf{T}$ to $T_{\{1,2\}\{3,4\}}$, and perform the eigenvalue truncation

where $+$ denotes projection to the PSD cone.

Compute the truncation of a different reshaping of $T^{>\varepsilon}$,

It follows that when we perform the eigenvalue truncation in step 1 of Algorithm 5.4,

We now remark that $\sigma$ maps $N$ to one of the bounded-norm matrices from Claim 5.6,

Furthermore, because the positive semidefinite cone is a closed convex set, and because projection to closed convex sets can only decrease distances (see Lemma A.2),

Now we prove that the spectral norms of the symmetrizations of the tensor have small spectral norm.

The first matrix that we are interested in is PSD, and can dominated by a tensor power of the identity:

Without loss of generality, let $i:=1$ so that $\langle b,a_{i}\rangle^{2}=1-\eta\geqslant 0.99$ (henceforth, we use $i$ as an ordinary index). We have that

An identical proof, up to signs, gives us a lower bound of $2p\alpha$.

For a matrix $E$ with $\|E\|\leqslant 4\alpha$.

It remains to argue that the top eigenvector of $M_{u}$ is $a_{1}$. We have that

where the last line follows because $w$ is a unit vector, and we chose $\varepsilon$ so that $1-\eta-\varepsilon>\eta$. Therefore

Applying Fact A.3, we conclude that if $v$ is the top eigenvector of $M_{u}$, then $\langle v,a_{1}\rangle^{2}\geqslant\frac{\varepsilon-8\alpha}{\varepsilon}\geqslant 1-\frac{8\alpha}{\varepsilon}$. Since we can choose $\varepsilon=\frac{1}{2}$ and still have that $\varepsilon<0.98<1-2\eta$, the conclusion follows. ∎

2 From independent columns to orthonormal columns

We now use standard techniques to prove that one can reduce from a dictionary with independent columns to a dictionary with orthogonal columns, given sufficiently many samples.

Note that the expected covariance matrix of the samples is equal to a scaled version of the covariance matrix,

So given sufficiently many samples, we can compute a good spectral approximation of $\Sigma^{-1/2}$, $\hat{\Sigma}^{-1/2}$ with

Now, define $S_{diff}=\left((\hat{\Sigma}^{-1/2})^{\otimes 2}-(\Sigma^{-1/2})^{\otimes 2}\right)$ and $S_{sum}=\left((\hat{\Sigma}^{-1/2})^{\otimes 2}+(\Sigma^{-1/2})^{\otimes 2}\right)$. We can factor the difference

Since we can choose the number of samples so as to make this last quantity as small as we would like, as a function of the condition number, and then appealing to Fact 4.10, the reduction is complete. ∎

3 Sample complexity bounds

Below is our bound on the sample complexity of approximating the $4$th moment tensor, which we believe may be loose.

This is because of the orthonormality of the $a_{i}$, which guarantees that terms in the product $MM^{\top}$ in which we have an inner product between two non-identical vectors drop out to .

References

Appendix A Useful Tools

We compute the expectation and variance of our matrix,

The two variance terms correspond to $\|A_{\{\alpha\beta\}\{\gamma\}}\|^{2}$ and $\|A_{\{\alpha\gamma\}\{\beta\}}\|^{2}$ respectively. We can now apply concentration results for matrix Gaussian series to conclude the proof [Oli10]. ∎

The following lemma states that distances can only decrease under projections to a convex set, and is well-known (see e.g. [Roc76]).

If we let $D_{x}=x-\Pi(x)$, $D_{y}=y-\Pi(y)$,

Now the conclusion will follow from the fact that

where the second inequality is the triangle inequality. Rearranging, the conclusion follows. ∎