PaperlayerView exact revision on arXiv

Source-linked AI summary

On Polynomial Time Methods for Exact Low Rank Tensor Completion

Dong Xia, Ming Yuan

arXiv:1702.06980v1stat.MLcs.ITcs.LGmath.IT

TL;DR

Exact completion of low-rank higher-order tensors from subsets of entries is delicate because spectral-norm computation is difficult. The paper proposes a polynomial-time gradient-descent approach initialized spectrally, with a stated sample-size bound and a dimensionality limitation beyond third order.

  • Problem

    Exact completion of low-rank higher-order tensors is difficult because spectral-norm computation becomes delicate.

    Source: §1 · p. 2

  • Method

    The paper uses gradient descent with an initial value from a spectral method and proposes an alternative approach to estimating singular spaces.

    Source: §1 · pp. 4–5

  • Results

    Polynomial-time algorithms reconstruct higher-order tensors using a sample-size requirement whose stated bound includes r1r2r3(rd1d2d3)1/2 log7/2 d and (r1r2r3)2rd log6 d terms.

    Source: §1 · p. 3

  • Takeaways & Limitations

    The approach provides a polynomial-time route to tensor reconstruction, while its dimensionality dependence is worse than nuclear norm minimization for order k > 3.

    Source: §1 · p. 3, §7 · p. 25

    For order k > 3, the dimensionality dependence is worse than the d3/2 dependence of nuclear norm minimization.

    Source: §7 · p. 25

Abstract

from arXiv · show

In this paper, we investigate the sample size requirement for exact recovery of a high order tensor of low rank from a subset of its entries. We show that a gradient descent algorithm with initial value obtained from a spectral method can, in particular, reconstruct a ${d\times d\times d}$ tensor of multilinear ranks $(r,r,r)$ with high probability from as few as $O(r^{7/2}d^{3/2}\log^{7/2}d+r^7d\log^6d)$ entries. In the case when the ranks $r=O(1)$, our sample size requirement matches those for nuclear norm minimization (Yuan and Zhang, 2016a), or alternating least squares assuming orthogonal decomposability (Jain and Oh, 2014). Unlike these earlier approaches, however, our method is efficient to compute, easy to implement, and does not impose extra structures on the tensor. Numerical results are presented to further demonstrate the merits of the proposed approach.

1 Introduction

This section introduces the challenge of higher-order tensor completion and presents an efficient spectral-initialization and Grassmannian gradient-descent approach to polynomial-time exact recovery. It also outlines the paper’s organization, including initialization, optimization, numerical experiments, discussions, and proofs.

  • Higher-order tensor completion is delicate because tensor spectral norm computation is NP-hard in general, while several matrix-completion tools remain underdeveloped.These challenges are both computational and theoretical.Source: §1 · p. 2
  • Unfolding-based methods are easy to implement but may require unnecessarily many entries, whereas nuclear-norm minimization uses as few as O((r^(k−1)/2d^3/2 + r^(k−1)d)(log d)^2) entries for k ≥3.Norm-minimization approaches are based on convex optimization but are NP-hard to compute.Source: §1 · p. 3
  • The paper develops polynomial-time recovery for a d1 × d2 × d3 tensor with multilinear ranks (r1, r2, r3) using r1r2r3(rd1d2d3)^1/2 log^7/2 d + (r1r2r3)^2rd log^6 d observations.Here r = max{r1, r2, r3} and d = max{d1, d2, d3}.Source: §1 · p. 3
  • The sample complexity matches tensor nuclear-norm minimization in its dependence on dimensions but is worse in rank dependence; for r = O(1), it becomes O(d^3/2 log^7/2 d), nearly identical to nuclear-norm recovery.The method is therefore aimed especially at high-dimensional, low-rank instances.Source: §1 · pp. 3–4
  • Starting from this initialization, Grassmannian gradient descent converges locally because the objective behaves like a parabola near the truth, where the true tensor is the only stationary point.The analysis requires new tensor concentration bounds building on Yuan and Zhang (2016a,b).Source: §1 · pp. 4–5
  • Initialization estimates unfolding singular spaces through a carefully constructed second-moment U-statistic whose eigenspace concentrates more sharply than ordinary unfolding SVD.The resulting spectral computation is based on a d × d matrix and is described as fairly efficient.Source: §1 · p. 4
  • Sections 4 and 5 detail the initialization and local optimization steps.Source: §1 · p. 5
  • Section 6 presents numerical experiments, followed by discussions and remarks in Section 7 and proofs of the main results in Section 8.Source: §1 · p. 5

2 Preliminaries

The preliminaries define multilinear rank through tensor unfoldings, represent low-rank tensors using factor matrices and a core tensor, and introduce coherence and tensor norms. They focus on multilinear ranks because these support stable numerical computation and can be related to CP rank.

Tensor representations

  • Multilinear ranks are the ranks of the tensor’s mode-specific matrix unfoldings, which rearrange tensor entries into matrices.The three unfolding maps are vector-space isomorphisms, and the ranks need not be equal across modes.Source: §2 · p. 6
  • A tensor with multilinear ranks (r1, r2, r3) can be represented by orthonormal factor matrices and a core tensor through multilinear multiplication.The factors lie on Stiefel manifolds, while the core has dimensions r1 × r2 × r3.Source: §2 · pp. 6–7

Rank notions

  • Multilinear rank is emphasized because it enables stable numerical computation, while the results can also be translated to CP-rank statements.The CP rank lies between max{r1, r2, r3} and min{r1r2, r1r3, r2r3}.Source: §2 · p. 7

Incoherence

  • Exact recovery from uniformly sampled entries requires low rank together with incoherence of the subspaces spanned by the tensor’s fibers.Tensor coherence is defined as the maximum coherence of its three mode-fiber subspaces.Source: §2 · pp. 7–8

Tensor norms

  • The section introduces inner-product, Frobenius, spectral, nuclear, and entrywise max norms for tensors.The nuclear norm is identified as the dual of the spectral norm.Source: §2 · pp. 8–9

3 Tensor Completion

Tensor completion is formulated as recovering a low-multilinear-rank, incoherent tensor from sampled entries by optimizing over factor subspaces and a core tensor. The proposed approach combines spectral initialization with local gradient-based optimization and yields a polynomial-time exact-recovery guarantee.

Sampling model

  • The framework focuses on the Bernoulli sampling model, while the analysis connects it to uniformly sampled subsets and motivates the choice through simpler technical arguments.The sampled entries are independently and uniformly drawn under the Bernoulli formulation.Source: §3 · pp. 9–10

Optimization formulation

  • The completion objective optimizes factor matrices and a core tensor so their reconstructed tensor matches the observed entries.The factor subspaces can be viewed as points on a direct product of Grassmannian manifolds.Source: §3 · pp. 10–11

Initialization

  • Gradient-based optimization needs a good initialization because the objective is highly nonconvex and iterative methods may converge only to local optima.Random initialization would require exponentially many seeds in d to ensure a good starting point.Source: §3 · p. 11

Initialization and convergence

  • The proposed spectral method efficiently produces an initial value sufficiently close to the optimum, overcoming the higher-order limitation that analogous spectral decompositions are NP hard to compute.The method then establishes local parabolic behavior of the objective around the truth, supporting local convergence of gradient descent.Source: §3 · pp. 11–12

Main result

  • Theorem 1 gives a polynomial-time exact-recovery algorithm for tensors with bounded multilinear rank, coherence, and condition number under independent uniform sampling.The guarantee holds with probability at least 1−d^-α when the sample size satisfies the stated rank-, dimension-, coherence-, condition-number-, and logarithmic-dependent bound.Source: §3 · p. 12

4 Second Order Method for Estimating Singular Spaces

The second-order method estimates tensor singular spaces by estimating an unfolding’s second moment rather than the unfolding itself. This improves behavior for unbalanced matrices and supplies consistent subspace estimates used for tensor initialization.

Tensor initialization

  • For tensor initialization, the leading singular vectors of the estimated unfolding second moments provide estimates of the factor spaces U, V, and W.The analysis measures the error between the estimated and target spaces using projection distance.Source: §4 · pp. 13–14

Unbalanced unfoldings

  • The second-order construction addresses the naive unfolding approach’s weakness when the unfolding is highly unbalanced and only the lower-dimensional singular space is needed.The naive estimate is driven largely by the larger dimension d2d3 even when the target space lies in dimension d1.Source: §4 · p. 13

Second-moment estimation

  • The method estimates MM^T from sampled matrix entries and extracts the leading left singular space from this second-moment estimate.Each rescaled single-entry observation is unbiased for M, and the resulting estimator is used to approximate MM^T.Source: §4 · pp. 13–14
  • Theorem 2 shows that the second-moment estimator can remain consistent even when directly estimating M in spectral norm is not.When m1 ≪ m2, the second-moment estimate remains consistent in a regime where the direct estimator is not consistent.Source: §4 · pp. 13–14

Tensor initialization

  • Corollary 1 establishes that the estimated factor space is consistent under bounded coherence and condition number with sufficiently many samples.The result holds with probability at least 1−d^-α, with analogous conclusions for the other modes.Source: §4 · p. 15

5 Exact Recovery by Optimizing Locally

The method optimizes a penalized objective over a product of Grassmannians, using spectral initialization and gradient descent to recover the tensor subspaces. Under stated conditions, the objective is locally well behaved and the algorithm converges to the true subspaces.

  • The algorithm operates on a product of Grassmannians using projection distance and geodesic updates for the three factor subspaces.The product space is J(d1, d2, d3, r1, r2, r3) = G(d1, r1) × G(d2, r2) × G(d3, r3).Source: §5 · pp. 15, 17
  • Theorem 3 establishes the local geometry under rank, conditioning, and sample-size conditions, with bounds involving r, r1, r2, r3, d1, d2, d3, and logarithmic factors.Source: §5 · p. 16
  • The objective behaves like a parabola near the true subspaces, where (U, V, W) is the unique stationary point.This local geometry enables gradient-based optimization within the neighborhood.Source: §5 · p. 16
  • A coherence penalty with sufficiently large ρ keeps iterates incoherent throughout gradient descent.Source: §5 · p. 17
  • With an appropriate initial value, Theorem 4 shows that Algorithm 1 converges to (U, V, W).The algorithm initializes from spectral estimates and repeatedly computes a negative gradient followed by geodesic updates.Source: §5 · pp. 17–18

6 Numerical Experiments

Numerical experiments evaluate recovery of cubic low-rank tensors and compare gradient descent with geometrical conjugate gradient descent. Recovery is strong with spectral initialization, while perturbing that initialization reduces reconstruction rates.

  • The experiments focus on cubic tensors with equal multilinear ranks and include orthogonally decomposable tensors whose core tensor is diagonal.Source: §6 · p. 19
  • Simulation results suggest recovery of cubic tensors from roughly r d^(3/2) observed entries, below the stated theoretical requirement n ≳ r^(7/2)d^(3/2).Source: §6 · p. 19
  • The proposed algorithm successfully recovers tensors when the relative Frobenius error is at most 10^-7, with experiments repeated over 50 runs.Source: §6 · p. 19
  • With spectral initialization, GoG and GeoCG perform similarly in the tested settings.The comparison uses d = 50, r = 5 and d = 100, r = 3.Source: §6 · pp. 22–24
  • Increasing perturbation noise in the initialization decreases reconstruction rates for GeoCG.The experiments use a random Gaussian perturbation of the spectral starting point.Source: §6 · pp. 22–23

7 Discussions

The discussion presents a polynomial-time spectral-initialized optimization framework for exact low-rank tensor completion, with sample complexity near d^3/2 for constant rank but weaker dimensional dependence for higher orders.

7 Discussions

  • The method combines a new spectral initialization based on second moments of tensor unfoldings with polynomial-time optimization over products of Grassmannians.The optimization is well behaved near the truth, so initialization largely drives the sample requirement.Source: §7 · p. 25
  • For r=O(1), the sample requirement is essentially d^3/2 up to polylogarithmic factors, matching nuclear norm minimization while remaining computationally efficient.Nuclear norm minimization has the same order but is NP hard to compute in general.Source: §7 · p. 25
  • The framework extends naturally to kth-order tensors through analogous Grassmannian optimization and unfolding-based initialization.The extension applies to tensors with equal multilinear rank and uses d × d^(k−1) unfoldings.Source: §7 · p. 25
  • For k>3, the dimensional dependence d^(k/2) is worse than the d^3/2 dependence of nuclear norm minimization, and its optimality remains unclear.The discussion leaves open whether this dependence is best attainable for polynomial-time algorithms.Source: §7 · p. 25

8 Proofs

This section proves the main theorems using spectral-initialization refinement, concentration and empirical-process bounds, and geometric control on Grassmannians. It establishes the local landscape and shows that gradient descent remains in the relevant neighborhood and converges to the unique stationary point.

Proof of Theorem 1

  • Theorem 1 follows by refining spectral estimates into an incoherent initialization in polynomial time.Lemma 2 gives an O(dr^2)-complexity procedure that increases incoherence by at most a factor of 3 and principal-angle distance by at most a factor of 4.Source: §8 · p. 26

Proof of Theorem 2

  • Theorem 2 is proved by decoupling the U-statistic sampling operator and separately controlling four error terms with Chernoff, Golden–Thompson, triangular, and matrix Bernstein inequalities.The proof combines the resulting probability bounds and adjusts the tail parameter to obtain the theorem’s final guarantee.Source: §8 · pp. 26, 27, 28, 31

Proof of Theorem 3

  • Theorem 3 analyzes tensor perturbations on Grassmannians using tangent spaces, geodesics, principal angles, and projection operators.The proof relates Grassmannian distance to principal angles, bounds multilinear perturbation terms, and controls empirical processes under incoherence and high-probability events.Source: §8 · pp. 32, 34, 37
  • With probability at least 1−e−t, the auxiliary empirical-process lemma supplies a uniform bound over γ1 and γ2 ranges.Specializing the parameters and setting t = α log d yields a probability of at least 1−d^−α over γ1 ∈ [(d1d2d3)^−1, 1] and γ2 ∈ [1, d].Source: §8 · p. 38
  • Combining the concentration estimates and projection bounds yields the stated high-probability control of the tensor perturbation and sampling operator.The final step combines equations (6), (8), and (13), then establishes the corresponding bound on event E1 ∩ E2 ∩ E3.Source: §8 · pp. 40–41

Proof of gradient bounds

  • The proof writes H = H1 + H2 and QTH1 = H1, then puts inequalities (19), (23), and (25) together on the event E1 ∩ E2 ∩ E3.Source: §8 · pp. 42, 46

Proof of Theorem 4

  • The additional penalty does not change the local behavior of F, so Theorem 3 and its first- and second-statement properties extend to the penalized objective ˜F.The proof reduces the remaining claim to bounding the penalty contribution, following the argument of Keshavan et al. (2009).Source: §8 · p. 46
  • For δ ≤ (Cακ0 log d)^−1, ˜F has a unique stationary point at (U, V, W), and Algorithm 1 converges to it by staying in the prescribed neighborhood.All iterates remain in N(δ/10, 4µ0), where the algorithm is gradient descent with exact line search.Source: §8 · p. 47

A Proof of Lemma 1

The proof reduces the tensor statements to its core tensor and bounds nuclear and Frobenius norms through successive mode slices.

A Proof of Lemma 1

  • Writing A=(U,V,W)·C preserves both the nuclear and Frobenius norms through the core tensor C.The core has dimensions determined by the multilinear ranks.Source: A Proof of Lemma 1 · p. 51
  • Convexity applied to the mode-1 slices of C provides the first norm bound, which is then repeated for modes 2 and 3.The same process on all three modes completes the argument.Source: A Proof of Lemma 1 · p. 51

B Proof of Corollary 1

The corollary proof applies the preceding theorem to a matrix unfolding and controls the resulting spectral bound with high probability.

B Proof of Corollary 1

  • Choosing m1=d1 and m2=d2d3 applies Theorem 2 to the unfolding dimensions of the tensor.The proof uses the condition n ≥ C1(α+1)(d1d2d3)^1/2.Source: B Proof of Corollary 1 · p. 51
  • The resulting bound holds with probability at least 1−d^−α after controlling the maximum norm under the incoherence condition.The proof explicitly recalls µ(T)≤µ0 before applying the bound.Source: B Proof of Corollary 1 · p. 52
  • The remaining steps use the definition of the relevant quantity and conclude the corollary.Source: B Proof of Corollary 1 · p. 52

C Proof of Lemma 3

The proof constructs a random coordinate-sampling tensor and uniformly bounds an empirical process over discretized parameter ranges using symmetrization and concentration inequalities.

  • The resulting bounds depend on logarithmic factors in d and hold under stated conditions involving the parameters and rank-related quantities.Source: C Proof of Lemma 3 · pp. 52–53
  • A random tensor E is defined with one randomly selected coordinate equal to 1, and independent copies E1, …, En support the empirical-process analysis.Source: C Proof of Lemma 3 · p. 52
  • The parameter ranges for γ1 and γ2 are partitioned into dyadic intervals, allowing bounds for each pair of intervals before a union argument makes them uniform.Source: C Proof of Lemma 3 · p. 52
  • Bousquet’s version of Talagrand’s concentration inequality supplies the high-probability control needed for the empirical process.Source: C Proof of Lemma 3 · p. 53
  • Symmetrization introduces i.i.d. Rademacher variables, while contraction uses the bound |⟨A, E⟩| ≤ γ1 to control the resulting process.Source: C Proof of Lemma 3 · p. 53
  • The argument follows related bounds such as Yuan and Zhang (2016a, Lemma 8).Source: C Proof of Lemma 3 · p. 53

D Proof of lower bound of ⟨QT(bT −T), H1⟩

The proof lower-bounds an inner product by decomposing it into a dominant squared-norm term and error terms, then controls the latter on a high-probability event.

Loading 1702.06980v1…