Square Deal: Lower Bounds and Improved Relaxations for Tensor Recovery
Cun Mu, Bo Huang, John Wright, Donald Goldfarb
Introduction
Tensors arise naturally in problems where the goal is to estimate a multi-dimensional object whose entries are indexed by several continuous or discrete variables. For example, a video is indexed by two spatial variables and one temporal variable; a hyperspectral datacube is indexed by two spatial variables and a frequency/wavelength variable. While tensors often reside in extremely high-dimensional data spaces, in many applications, the tensor of interest is low-rank, or approximately so [KB09], and hence has much lower-dimensional structure. The general problem of estimating a low-rank tensor has applications in many different areas, both theoretical and applied: e.g., estimating latent variable graphical models [AGH+12], classifying audio [MSS06], mining text [CC12], processing radar signals [DN10], to name a few.
In contrast, the correct generalization of these results to low-rank tensors is not obvious. The numerical algebra of tensors is fraught with hardness results [HL09]. For example, even computing a tensor’s (CP) rank,
is NP-hard in general. The nuclear norm of a tensor is also intractable, and so we cannot simply follow the formula that has worked for vectors and matrices.
With an eye towards numerical computation, many researchers have studied how to estimate or recover tensors of small Tucker rank [Tuc66]. The Tucker rank of a -way tensor is a -dimensional vector whose -th entry is the (matrix) rank of the mode- unfolding of :
The definition (1.2) suggests a very natural, tractable convex approach to recovering low-rank tensors: seek the that minimizes out of all satisfying . We will refer to this as the sum-of-nuclear-norms (SNN) model. Originally, proposed in [LMWY09], this approach has been widely studied [GRY11, SDS10, THK10, TSHK11, STDLS13] and applied to various datasets in imaging [SVdPDMS11, SHKM13, KS13, LL10, LYZY10].
Our theoretical results pertain to Gaussian operators . The motivation for studying Gaussian measurements is twofold. First, Gaussian measurements may be of interest for compressed sensing recovery [Don06], either directly as a measurement strategy, or indirectly due to universality phenomena [DT09, BLM12]. Second, the available theoretical tools for Gaussian measurements are very sharp, allowing us to rigorously investigate the efficacy of various regularization schemes, and prove both upper and lower bonds on the number of observations required. In simulation, our qualitative conclusions carry over to more realistic measurement models, such as random subsampling [LMWY09] (see Section 5). We expect our results to be of interest for a wide range of problems in tensor completion [LMWY09], robust tensor recovery / decomposition [LYZY10, GQ12] and sensing.
Our technical methodology draws on, and enriches, the literature on general structured model recovery. The surprisingly poor behavior of the SNN model is an example of a phenomenon first discovered by Oymak et. al. [OJF+12]: for recovering objects with multiple structures, a combination of structure-inducing norms is often not significantly more powerful than the best individual structure-inducing norm. Our lower bound for the SNN model follows from a general result of this nature, which we prove using the geometric framework of [ALMT13]. Compared to [OJF+12], our result pertains to a more general family of regularizers, and gives sharper constants. In addition, we demonstrate the possibility to reduce the number of generic measurements through a new convex regularizer that exploit several sparse structures jointly.
Bounds for Non-Convex Recovery
In this section, we introduce a non-convex model for tensor recovery, and show that it recovers low-rank tensors from near-minimal numbers of measurements. While our nonconvex formulation is computationally intractable, it gives a baseline for evaluating tractable (convex) approaches.
In vector optimization, a feasible point is called Pareto optimal if no other feasible point dominates it in every criterion. In a similar vein, we say that (2.1) recovers if there does not exist any other tensor that is consistent with the observations and has no larger rank along each mode:
We call recoverable by (2.1) if the set
This is equivalent to saying that is the unique optimal solution to the scalar optimization:
Whenever , with probability one, , and hence (2.1) recovers every .
The proof of Theorem 1 follows from a covering argument, which we establish in several steps. Let
The following lemma shows that the required number of measurements can be bounded in terms of the exponent of the covering number for , which can be considered as a proxy for dimensionality:
Suppose that the covering number for with respect to Frobenius norm, satisfies
It just remains to find the covering number of . We use the following lemma, which uses the triangle inequality to control the effect of perturbations in the factors of the Tucker decomposition
where the mode- (matrix) product of tensor with matrix of compatible size, denoted as , outputs a tensor such that .
Using this result, we construct an -net for by building -nets for each of the factors and . The total size of the resulting net is thus bounded by the following lemma:
With these observations in hand, Theorem 1 follows immediately.
Convexification: Sum of Nuclear Norms?
Since the nonconvex problem (2.1) is NP-hard for general , it is tempting to seek a convex surrogate. In matrix recovery problems, the nuclear norm is often an excellent convex surrogate for the rank [Faz02, RFP10, Gro11]. It seems natural, then, to replace the ranks in (2.1) with nuclear norms, and solve
Since is a convex function, the set is convex. For any pareto optimal point , there is a hyperplane supporting passing through , with normal vector . Therefore, is an optimal solution to the following scalar optimization:
Suppose that has Tucker rank , and . With high probability, is an optimal solution to (3.2), with each . Here, is numerical.
This result shows that there is a range in which (3.2) succeeds: loosely, when we undersample by at most a factor of . However, the number of observations is significantly larger than the number of degrees of freedom in , which is on the order of . Is it possible to prove a better bound for this model? Unfortunately, we show that in general measurements are also necessary for reliable recovery using (3.2):
Let be nonzero. Set . Then if the number of measurements , is not the unique solution to (3.2), with probability at least . Moreover, there exists for which .
This implies that Corollary 2 (and other results of [THK10]) is essentially tight. Unfortunately, it has negative implications for the efficacy of the sum of nuclear norms in (3.2): although a generic element of can be described using at most real numbers, we require observations to recover it using (3.2). Theorem 3 is a direct consequence of a much more general principle underlying multi-structured recovery, which is elaborated next.
The poor behavior of (3.2) is actually an instance of a much more general phenomenon, first discovered by Oymak et. al. [OJF+12]. Our target tensor has multiple low-dimensional structures simultaneously: it is low-rank along each of the modes. In practical applications, many other such simultaneously structured objects may be of interest – for example, matrices that are simultaneously sparse and low-rank [RSV12, OJF+12]. To recover such a simultaneously structured object, it is tempting to build a convex relaxation by combining the convex relaxations for each of the individual structures. In the tensor case, this yields (3.2). Surprisingly, this combination is often not significantly more powerful than the best single regularizer [OJF+12]. We obtain Theorem 3 as a consquence of a new, general result of this nature, using a geometric framework introduced in [ALMT13]. Compared to the proof strategy in [OJF+12], this approach has a clearer geometric intuition, covers a more general class of regularizers and yields sharper bounds.
where is a Gaussian measurement operator, and . Is the unique optimal solution to (3.3)? Recall that the descent cone of a function at a point is defined as
To control the size of , first consider a single norm , with dual norm . Suppose that is -Lipschitz: for all . Then for all as well. Noting that
for any , we have
A more geometric way of summarizing this is as follows: for , let
and denote the circular cone with axis and angle . Then if , and ,
For , notice that every element of is a conic combination of elements of the . Since each of the is contained in a circular cone with axis , is also contained in a circular cone:
Suppose that is -Lipschitz. For , set . Then
So, the subdifferential of our combined regularizer is contained in a circular cone whose angle is given by the largest of the .
How does this behavior affect the recoverability of via (3.3)? The informal reasoning above suggests that as becomes smaller, the descent cone becomes larger, and we require more measurements to recover . This can be made precise using an elegant framework introduced by Amelunxen et. al. [ALMT13]. They define the statistical dimension of the convex cone to be the expected norm of the projection of a standard Gaussian vector onto :
Using tools from spherical integral geometry, [ALMT13] shows that for linear inverse problems with Gaussian measurements, a sharp phase transition in recoverability occurs around . We will need only one side of their result; for more details see [ALMT13]. We state a slight variant here:
To apply this result to our problem, we lower bound the statistical dimension , of the descent cone of at . Using the Pythagorean theorem, monotonicity of , and Lemma 4, we calculate
Moreover, using the properties of statistical dimension, we are able to prove an upper bound for the statistical dimension of circular cone, which improves the constant in existing results [ALMT13, McC13].
Finally, by combining (3.11) and Lemma 5, we have . Using Corollary 4, we obtain:
Let . Suppose that for each , is -Lipschitz. Set
and . Then if ,
Thus, for reliable recovery, the number of measurements needs to be at least proportional to .E.g., if , the probability of success is at most . Notice that is determined by only the best of the structures. Per Table 1, is often on the order of the number of degrees of freedom in a generic object of the -th structure. For example, for a -sparse vector whose nonzeros are all of the same magnitude, .
Theorem 5 together with Table 1 leads us to the phenomenon that recently discovered by Oymak et. al. [OJF+12]: for recovering objects with multiple structures, a combination of structure-inducing norms tends to be not significantly more powerful than the best individual structure-inducing norm. As we demonstrate, this general behavior follows a clear geometric interpretation that the subdifferential of a norm at is contained in a relatively small circular cone with central axis .
We can specialize Theorem 5 to low-rank tensors as follows: if is a -mode tensor of Tucker rank , then for each , is -Lipschitz. Hence,
A Better Convexification: Square Norm
The number of measurements promised by Corollary 2 and Theorem 3 is actually the same (up to constants) as the number of measurements required to recover a tensor which is low-rank along just one mode. Since matrix nuclear norm minimization correctly recovers a matrix of rank when [CRPW12], solving
also exactly recovers with high probability when .
This suggests a more mundane explanation for the difficulty with (3.2): the term comes from the need to reconstruct the right singular vectors of the matrix . If we had some way of matricizing a tensor that produced a more balanced (square) matrix and also preserved the low-rank property, we could substantially reduce this effect, and reduce the overall sampling requirement. In fact, this is possible when the order of is four or larger.
We can view as a natural generalization of the standard tensor matricization. When , is nothing but . However, when some is selected, becomes a more balanced matrix. This reshaping also preserves some of the algebraic structures of . In particular, we will see that if is a low-rank tensor (in either the CP or Tucker sense), will be a low-rank matrix.
(1) If has CP decomposition , then
(2) If has Tucker decomposition , then
Thus, is not only more balanced but also maintains the low-rank property of tensor . In the following, we show how this new matricization can lead to better relaxations for tensor recovery. For ease of discussion, we assume has the same length (say ) along each mode and has Tucker rank . We write and call the square norm of tensor . Since is low-rank, we can attempt to recover by solving
Using Lemma 7 and Proposition 3.11 of [CRPW12], we can prove that this relaxation exactly recovers , when the number of measurements is sufficienly large:
(1) If has CP rank , using (4.4), is sufficient to recover with high probability. (2) If has Tucker rank , using (4.4), is sufficient to recover with high probability.
Compared with measurements required by the sum-of-nuclear-norms model, the sample complexity, , required by the square reshaping (4.4), is always within a constant of it, much better for small and – e.g., by a multiplicative factor of when is a constant. This is a significant improvement. However, there are also two clear limitations. First, no improvement is obtained for the case . Second, the improved sample complexity in Theorem 6 is still suboptimal compared to the nonconvex model (2.1).
It is also worth noting that for tensors with different lengths or ranks, Theorem 3 and Theorem 6 can be easily modified. It remains true that for a large class of tensors, our square reshaping is capable of reducing the number generic measurements required by SNN model. However, the comparison between sum-of-nuclear-norms and square norm becomes quite subtle then. Concrete instances can be definitely constructed so that square norm model does not have any advantage over the SNN model even for (e.g. a tensor of size with Tucker rank ). On the other hand, our square norm model can sometimes be blessed by unbalanced tensors. For example, consider a tensor of size with Tucker rank . Then our reshaping matrix is a square matrix with rank , which is a matrix with very good (perfect) conditions.
Simulation Results for Tensor Completion
Tensor completion attempts to reconstruct the low-rank tensor based on observations over a subset of its entries . By imposing appropriate incoherence conditions (and modifying slightly arguments in [Gro11]), it is possible to prove recovery guarantees for each of the following programs:
Unlike the recovery problem under Gaussian random measurements, due to the lack of sharp upper bounds, we have no proof that our square norm formulation outperforms the SNN model here. However, our simulation results below strongly suggest that (5.2) also performs much better than (5.1) for tensor completion case.
Conclusion
In this paper, we establish several theoretical bounds for the problem of low-rank tensor recovery using random Gaussian measurements. For the nonconvex model (2.1), we show that \big{(}(2r)^{K}+2nrK+1\big{)} measurements are sufficient to recover any almost surely. We highlight that though the nonconvex recovery program is NP-hard in general, it does serve a baseline for evaluating tractable (convex) approaches. For the conventional convex surrogate sum-of-nuclear-norms (SNN) model (3.2), we prove a necessary condition that Gaussian measurements are required for reliable recovery. This lower bound is derived from our study on multi-structured object recovery under a very general setting, which can be applied to many scenarios. To narrow the apparent gap between the non-convex model and the SNN model, we unfold the tensor into a more balanced matrix while preserving its low-rank property, leading to our square-norm model (4.4). We prove that measurements are sufficient to recover a tensor with high probability. Though the theoretical results only pertain to Gaussian measurements, our simulation result for tensor completion also suggests that square-norm model outperforms the SNN model.
Compared with measurements required by the sum-of-nuclear-norms model, the sample complexity, , required by the square reshaping (4.4), is always within a constant of it, much better for small and . Although this is a significant improvement, compared with the nonconvex model (2.1), the improved sample complexity achieved by square norm model is however still suboptimal. It remails an open problem to obtain near-optimal convex relaxations for all .
More broadly speaking, to recover objects with multiple structures, regularizing with a combination of individual structure-inducing norms is proven to be substantially suboptimal (Theorem 5 and also [OJF+12]). The resulting sample requirements tend to be much larger than the intrinsic degrees of freedom of the low-dimensional manifold that the structured signal lies in. Our square-norm model for the low-rank tensor recovery demonstrates the possibility that a better exploitation in those structures can significantly reduce this sample complexity. However, there are still no clear clues on how to intelligently utilize several simultaneous structures generally, and moreover how to design tractable method to recover multi-structured objects with near minimal number of measurements. These problems are definitely worth pursuing in future study and we hope that our work may also inspire researchers working in many other multi-structured recovery problems.
Acknowledgment
It is a great pleasure to acknowledge conversations with Michael McCoy (Caltech), Ju Sun (Columbia), Han-wen Kuo (Columbia), Martin Lotz (Manchester), Zhiwei Qin (WalmartLabs). JW was supported by Columbia University startup funding and Office of Naval Research award N00014-13-1-0492.
References
Appendix A Proofs for Section 2
The arguments we used below are primarily adapted from [ENP11], where their interest is to establish the number of Gaussian measurements required to recover a low rank matrix by rank minimization.
Notice that every , and every , is a standard Gaussian random variable, and so
Let be an -net for in terms of . Because the measurements are independent, for any fixed ,
Moreover, for any , we have
This follows from the basic fact that for any tensor and matrix of compatible size,
which can be established by direct calculation. Write
where the first inequality follows from triangle inequality and the second inequality follows from the fact that , , and . ∎
The idea of this proof is to construct a net for each component of the Tucker decomposition and then combine those nets to form a compound net with the desired cardinality.
Clearly . The rest is to show that is indeed an -net covering with respect to the Frobenius norm.
For any fixed where and , by our constructions above, there exist and such that and . Then is within -distance from , since by the triangle inequality derived in Lemma 2, we have
Appendix B Proofs for Section 3
Denote . Then following [ALMT13, Thm. 7.1], we have
Summing up the above inequalities, we have
Suppose is odd. Since the intersection of with any -dimensional linear subspace containing is an isometric image of , by [ALMT13, Prop. 4.1], we have
Thus, taking both cases ( is even and is odd) into consideration, we have
Notice that for any fixed , the function is decreasing for . Then due to Corollary 4 and the fact that , we have
Appendix C Proofs for Section 4
(1) By the definition of , it is sufficient to prove that the vectorization of the right hand side of (4.2) equals .
Since , we have
where the last equality follows from the fact that . Similarly, we can derive that the vectorization of the right hand side of (4.2),
Thus, equation (4.2) is valid. (2) The above argument can be easily adapted to prove the second claim. Since , we have
where the last equality follows from the fact that . Similarly, we can derive that the vectorization of the right hand side of (4.3),
Appendix D Algorithms for Section 5
In this section, we will discuss in detail our implementation of accelerated linearized Bregman algorithm for the following problem:
By introducing auxiliary variable and splitting into , it can be easily verified that problem (D.1) is equivalent to
whose objective function is now separable.
where we denote as the dual variable for the constraint and denote as the dual variable for the last constraint . Since the objective function in (D) is separable, each setup of the ALB algorithm is easy to solve as we can see from Algorithm 1 The Shrinkage operator in line 4 of Algorithm 1 performs the regular shrinkage on the singular values of the th unfolding matrix of , i.e. , and then folds the resulting matrix back into tensor..
For our numerical experiment (), we choose smoothing parameter and step size . Empirically, we observe that larger values of do not result in a better recovery performance. This is consistent with the theoretical results established in [LY12, ZCCZ11].