Successive Rank-One Approximations for Nearly Orthogonally Decomposable Symmetric Tensors

TL;DR

This paper demonstrates that successive rank-one approximations can robustly recover the symmetric canonical decomposition of nearly orthogonally decomposable tensors, even with small perturbations, without error accumulation over iterations.

Contribution

It extends the SROA method to nearly SOD tensors, proving robustness against perturbations and showing errors do not accumulate during iterations.

Findings

SROA effectively recovers tensor decomposition with small perturbations.

Approximation errors remain bounded and do not accumulate over iterations.

Numerical experiments support theoretical robustness results.

Abstract

Many idealized problems in signal processing, machine learning and statistics can be reduced to the problem of finding the symmetric canonical decomposition of an underlying symmetric and orthogonally decomposable (SOD) tensor. Drawing inspiration from the matrix case, the successive rank-one approximations (SROA) scheme has been proposed and shown to yield this tensor decomposition exactly, and a plethora of numerical methods have thus been developed for the tensor rank-one approximation problem. In practice, however, the inevitable errors (say) from estimation, computation, and modeling necessitate that the input tensor can only be assumed to be a nearly SOD tensor---i.e., a symmetric tensor slightly perturbed from the underlying SOD tensor. This article shows that even in the presence of perturbation, SROA can still robustly recover the symmetric canonical decomposition of the underlying tensor. It is shown that when the perturbation error is small enough, the approximation errors do not accumulate with the iteration number. Numerical results are presented to support the theoretical findings.

Figures

Equations

$$({\lambda }_1,{\mathbf{x}}_1)\in \underset{\lambda \in \mathbb{R},\Vert \mathbf{x}\Vert =1}{\mathrm{argmin}}{\Vert \mathbf{X}-\lambda \mathbf{x}{\mathbf{x}}^{\top }\Vert }_F;$$

$${\mathbf{v}}^{\otimes p}:=\underset{\text{p times}}{\underbrace{\mathbf{v}\otimes \mathbf{v}\otimes \cdots \otimes \mathbf{v}}},$$

$$({\lambda }_1,{\mathbf{v}}_1)\in \underset{\lambda \in \mathbb{R},\Vert \mathbf{v}\Vert =1}{\mathrm{argmin}}{\Vert \mathcal{T}-\lambda {\mathbf{v}}^{\otimes p}\Vert }_F,$$

$$\hat{\mathcal{T}}=\mathcal{T}+\mathcal{E},$$

$$\mathcal{A}=\left({\mathcal{A}}_{i_1,i_2,\dots ,i_p}\right),{\mathcal{A}}_{i_1,i_2,\dots ,i_p}\in \mathbb{R},1\le i_1,i_2,\dots ,i_p\le n,$$

$${\mathcal{A}}_{i_1,i_2,\dots ,i_p}={\mathcal{A}}_{i_{\pi (1)},i_{\pi (2)},\dots ,i_{\pi (p)}}$$

$${\left(\mathcal{A}({\mathbf{V}}_1,{\mathbf{V}}_2,\dots ,{\mathbf{V}}_p)\right)}_{i_1,i_2,\dots ,i_p}:=\sum _{j_1,j_2,\dots ,j_p\in [n]}{\mathcal{A}}_{j_1,j_2,\dots ,j_p}({\mathbf{V}}_1{)}_{j_1{i}_1}({\mathbf{V}}_2{)}_{j_2{i}_2}\cdots ({\mathbf{V}}_p{)}_{j_p{i}_p}.$$

$$\mathcal{A}({\mathbf{V}}_1,{\mathbf{V}}_2)={\mathbf{V}}_1^{\top }\mathcal{A}{\mathbf{V}}_2\in {\mathbb{R}}^{m_1\times m_2}.$$

$$\mathcal{A}({\mathbf{e}}_{i_1},{\mathbf{e}}_{i_2},\dots ,{\mathbf{e}}_{i_p})={\mathcal{A}}_{i_1,i_2,\dots ,i_p},$$

$$\mathcal{A}{\mathbf{x}}^{\otimes p}:=\mathcal{A}(\mathbf{x},\mathbf{x},\dots ,\mathbf{x})=\sum _{i_1,i_2,\dots ,i_p\in [n]}{\mathcal{A}}_{i_1,i_2,\dots ,i_p}{x}_{i_1}{x}_{i_2}\cdots x_{i_p},$$

$$\mathcal{A}{\mathbf{x}}^{\otimes p-1}$$

$${\left(\mathcal{A}{\mathbf{x}}^{\otimes p-1}\right)}_i$$

$$⟨\mathcal{A},\mathcal{B}⟩:=\sum _{i_1,i_2,\dots ,i_p\in [n]}{\mathcal{A}}_{i_1,i_2,\dots ,i_p}{\mathcal{B}}_{i_1,i_2,\dots ,i_p}.$$

$$\mathbb{P}\left[h=i\right]=w_i\text{for all i∈[n]}.$$

$$({\mathbf{M}}_2{)}_{i_1,i_2}$$

$$({\mathcal{M}}_p{)}_{i_1,i_2,\dots ,i_p}$$

$${\mathbf{M}}_2=\sum _{i\in [n]}{w}_i{\mathbf{\mu }}_i{\mathbf{\mu }}_i^{\top }\text{and}{\mathcal{M}}_p=\sum _{i\in [n]}{w}_i{\mathbf{\mu }}_i^{\otimes p}.$$

$$\mathbf{W}:=\mathbf{U}{\mathbf{D}}^{-1/2},{\lambda }_i:=w_i^{1-p/2},\text{and}{\mathbf{v}}_i:=\sqrt{w_i}{\mathbf{W}}^{\top }{\mathbf{\mu }}_i\in {\mathbb{R}}^n\text{for each i∈[n]}.$$

$${\mathbf{W}}^{\top }{\mathbf{M}}_2\mathbf{W}=\mathbf{I}={\mathbf{M}}_2(\mathbf{W},\mathbf{W})=\sum _{i\in [n]}{w}_i({\mathbf{W}}^{\top }{\mathbf{\mu }}_i)({\mathbf{W}}^{\top }{\mathbf{\mu }}_i{)}^{\top }=\sum _{i\in [n]}{\mathbf{v}}_i{\mathbf{v}}_i^{\top },$$

$$\mathcal{T}:={\mathcal{M}}_p(\mathbf{W},\mathbf{W},\dots ,\mathbf{W})=\sum _{i\in [n]}{w}_i({\mathbf{W}}^{\top }{\mathbf{\mu }}_i{)}^{\otimes p}=\sum _{i\in [n]}{\lambda }_i{\mathbf{v}}_i^{\otimes p}.$$

$${\hat{\mathbf{v}}}_i\in \underset{\Vert \mathbf{v}\Vert =1}{\mathrm{argmax}}{\hat{\mathcal{T}}}_{i-1}{\mathbf{v}}^{\otimes p},{\hat{\lambda }}_i={\hat{\mathcal{T}}}_{i-1}{\hat{\mathbf{v}}}_i^{\otimes p}.$$

$$(\hat{\lambda },\hat{\mathbf{x}})\in \underset{\lambda \in \mathbb{R},\Vert \mathbf{x}\Vert =1}{\mathrm{argmin}}{\Vert \hat{\mathbf{M}}-\lambda \mathbf{x}{\mathbf{x}}^{\top }\Vert }_F.$$

$$|\hat{\lambda }-{\lambda }_j|\le \epsilon ,{\Vert \hat{\mathbf{x}}-{\mathbf{v}}_j\Vert }_2\le 10\left(\frac{\epsilon }{{\lambda }_j}+{\left(\frac{\epsilon }{{\lambda }_j}\right)}^2\right).$$

$${\lambda }_1|x_1{|}^{p-2}\ge {\lambda }_2|x_2{|}^{p-2}\ge \cdots \ge {\lambda }_n|x_n{|}^{p-2}.$$

$${\lambda }_1\ge {\lambda }_{\max }-2\epsilon ,|x_1|\ge 1-2\epsilon /{\lambda }_1,x_1^2\ge x_1^{p-1}\ge 1-4\epsilon /{\lambda }_1,\text{and}|\hat{\lambda }-{\lambda }_1|\le \epsilon .$$

$$\hat{\lambda }=\hat{\mathcal{T}}{\hat{\mathbf{x}}}^{\otimes p}$$

$$=\sum _{i=1}^n{\lambda }_i{x}_i^p+\mathcal{E}{\hat{\mathbf{x}}}^{\otimes p}$$

$$\le \sum _{i=1}^n{\lambda }_i|x_i{|}^{p-2}{x}_i^2+\epsilon \le {\lambda }_1|x_1{|}^{p-2}+\epsilon ,$$

$$\hat{\lambda }\ge \underset{i\in [n]}{\max }\mathcal{T}{{\mathbf{v}}_i}^{\otimes p}-\Vert \mathcal{E}\Vert ={\lambda }_{\max }-\epsilon \ge {\lambda }_1-\epsilon .$$

$${\lambda }_1\ge {\lambda }_{\max }-2\epsilon ,|{\lambda }_1-\hat{\lambda }|\le \epsilon .$$

$$|x_1|\ge |x_1{|}^{p-2}\ge \frac{{\lambda }_1-2\epsilon }{{\lambda }_1}=1-\frac{2\epsilon }{{\lambda }_1}$$

$$\mathcal{L}(\mathbf{x},\lambda ):=\hat{\mathcal{T}}{\mathbf{x}}^{\otimes p}-\frac{\lambda }{2}\left({\Vert \mathbf{x}\Vert }^2-1\right),$$

Full text

\slugger

simaxxxxxxxxx–x

1 Introduction

The eigenvalue decomposition of symmetric matrices is one of the most important discoveries in mathematics with an abundance of applications across all disciplines of science and engineering. One way to explain such a decomposition is to express the symmetric matrix as a minimal sum of rank-one symmetric matrices. It is well known that the eigenvalue decomposition can be simply obtained via successive rank-one approximation (SROA). Specifically, for a symmetric matrix $\bm{X}\in\mathbb{R}^{n\times n}$ with rank $r$, one approximates $\bm{X}$ by a rank-one matrix to minimize the Frobenius norm error:

[TABLE]

then, one approximates the residual $\bm{X}-\lambda_{1}\bm{x}_{1}\bm{x}_{1}^{\scriptscriptstyle{\top}}$ by another rank-one matrix $\lambda_{2}\bm{x}_{2}\bm{x}_{2}^{\scriptscriptstyle{\top}}$, and so on. The above procedure continues until one has found $r$ rank-one matrices $(\lambda_{i}\bm{x}_{i}\bm{x}_{i}^{\scriptscriptstyle{\top}})_{i=1}^{r}$; their summation, $\sum_{i=1}^{r}\lambda_{i}\bm{x}_{i}\bm{x}_{i}^{\scriptscriptstyle{\top}}$, yields an eigenvalue decomposition of $\bm{X}$. Moreover, due to the optimal approximation property of the eigendecomposition, for any positive integer $k\leq r$, the best rank-$k$ approximation (in the sense of either the Frobenius norm or the operator norm) to $\bm{X}$ is simply given by $\sum_{i=1}^{k}\lambda_{i}\bm{x}_{i}\bm{x}_{i}^{\scriptscriptstyle{\top}}$ [(1)].

In this article, we study decompositions of higher-order symmetric tensors, a natural generalization of symmetric matrices. Many applications in signal processing, machine learning, and statistics, involve higher-order interactions in data; in these cases, higher-order tensors formed from the data are the primary objects of interest. A tensor $\bm{\mathcal{T}}\in\bigotimes_{i=1}^{p}\mathbb{R}^{n_{i}}:=\mathbb{R}^{n_{1}\times n_{2}\times\dotsb\times n_{p}}$ of order $p$ is called symmetric if $n_{1}=n_{2}=\dotsb=n_{p}=n$ and its entries are invariant under any permutation of their indices. Symmetric tensors of order two ($p=2$) are symmetric matrices. In the sequel, we reserve the term tensor (without any further specifiction) for tensors of order $p\geq 3$. A symmetric rank-one tensor can be naturally defined as a $p$-fold outer product

[TABLE]

where $\bm{v}\in\mathbb{R}^{n}$ and $\otimes$ denotes the outer product between vectors.111For any $1\leq i_{1},i_{2},\dotsc,i_{p}\leq n$ and any $\bm{v}\in\mathbb{R}^{n}$, the $(i_{1},i_{2},\dotsc,i_{p})$-th entry of $\bm{v}^{\otimes p}$ is $(\bm{v}^{\otimes p})_{i_{1},i_{2},\dotsc,i_{p}}=v_{i_{1}}v_{i_{2}}\dotsb v_{i_{p}}$.

The minimal number of rank-one symmetric tensors whose sum is $\bm{\mathcal{T}}$ is called the symmetric tensor rank in the literature, and any corresponding decomposition is called a symmetric canonical decomposition [(2)]. Such decompositions have applications in many scientific and engineering domains [(3), (4), (5), (6), (7), (8)].

By analogy to the matrix case, successive rank-one approximations schemes have been proposed for symmetric tensor decomposition [(9), (10), (11), (12)]. Just as in the matrix case, one first approximates $\bm{\mathcal{T}}$ by a symmetric rank-one tensor

[TABLE]

and then approximate the residual $\bm{\mathcal{T}}-\lambda_{1}\bm{v}_{1}^{\otimes p}$ again by another symmetric rank-one tensor, and so on. (The Frobenius norm $\left\|\cdot\right\|_{F}$ for tensors is defined later, but is completely analogous to the matrix case.) This process continues until a certain stopping criterion is met. However, unlike symmetric matrices, the above procedure for higher order tensors ($p\geq 3$) faces a number of computational and theoretical challenges.

Unlike problem (1.1)—which can be solved efficiently using simple techniques such as power iterations— solving the rank-one approximation to higher order tensors is much more difficult: it is NP-hard, even for symmetric third-order tensors [(13)]. Researchers in numerical linear algebra and numerical optimization have devoted a great amount of effort to solve problem (1.2). Broadly speaking, existing methods for problem (1.2) can be categorized into three types. First, as problem (1.2) is equivalent to finding the extreme value of a homogeneous polynomial over the unit sphere, general-purpose polynomial solvers based on the Sum-of-Squares (SOS) framework [(14), (15), (16), (17), (18)], such as GloptiPoly 3 [(19)] and SOSTOOLS [(20)], can be effectively applied to the rank-one approximation problem. The SOS approach can solve any polynomial problem to any given accuracy through a sequence of semidefinite programs; however, the size of these programs are very large for high-dimensional problems, and hence these techniques are generally limited to relatively small-sized problems. The second approach is to treat problem (1.2) as a nonlinear program [(21), (22)], and then to exploit and adapt the wealth of ideas from numerical optimization. The resulting methods—which include [(23), (9), (10), (12), (24), (25), (26), (27), (28)] to just name a few—are empirically efficient and scalable, but are only guaranteed to reach a local optimum or stationary point of the objective over the sphere. Therefore, to maximize their performance, these methods need to run with several starting points. The final approach is based on the recent trend of relaxing seemingly intractable optimzation problems such as problem (1.2) with more tractable convex optimization problems that can be efficiently solved [(29), (30), (31)]. The tensor structure in (1.2) has made it possible to design highly-tailored convex relaxations that appear to be very effective. For instance, the semidefinite relaxation approach in  [(29)] was able to globally solve almost all the randomly generated instances that they tested. Aside from the above solvers, a few algorithms have been specifically designed for the scenario where some side information regarding the solution of (1.2) is known. For example, when the signs of the optimizer $\bm{v}_{1}$ are revealed, polynomial time approximation schemes for solving (1.2) are available [(32)].

In contrast to the many active efforts and promising results on the computational side, the theoretical properties of successive rank-one approximations are far less developed. Although SROA is justified for matrix eigenvalue decomposition, it is known to fail for general tensors [(33)]. Indeed, much has been established about the failings of low-rank approximation concepts for tensors that are taken for granted in the matrix setting [(34), (35), (36), (37), (38)]. For instance, the best rank-$r$ approximation to a general tensor is not even guaranteed to exist (though several sufficient conditions for this existence have been recently proposed [(39), (40)]). Nevertheless, SROA can be still justified for certain classes of symmetric tensors that arise in applications.

Nearly SOD tensors

Indeed, in many applications (e.g., higher-order statistical estimation [(3)], independent component analysis [(4), (7)], and parameter estimation for latent variable models [(8)]), the input tensor $\widehat{\bm{\mathcal{T}}}$ may be fairly assumed to be a symmetric tensor slightly perturbed from a symmetric and orthogonally decomposable (SOD) tensor $\bm{\mathcal{T}}$. That is,

[TABLE]

where the underlying SOD tensor may be written as $\bm{\mathcal{T}}=\sum_{i=1}^{r}\lambda_{i}\bm{v}_{i}^{\otimes p}$ with $\left\langle\bm{v}_{i},\bm{v}_{j}\right\rangle=\mathds{1}\{i=j\}$ for $1\leq i,j\leq r$, and $\bm{\mathcal{E}}$ is a perturbation tensor. In these aforementioned applications, we are interested in obtaining the underlying pairs $\{(\lambda_{i},\bm{v}_{i})\}_{i=1}^{r}$. When $\bm{\mathcal{E}}$ vanishes, it is known that $\sum_{i=1}^{r}\lambda_{i}\bm{v}_{i}^{\otimes p}$ is the unique symmetric canonical decomposition [(2), (41)], and moreover, successive rank-one approximation exactly recovers $\{(\lambda_{i},\bm{v}_{i})\}_{i=1}^{r}$ [(9)]. However, because of the inevitable perturbation term $\bm{\mathcal{E}}$ arising from sampling errors, noisy measurements, model misspecification, numerical errors, and so on, it is crucial to understand the behavior of SROA when $\bm{\mathcal{E}}\neq\bm{0}$. In particular, one may ask if SROA provides an accurate approximation to $\{(\lambda_{i},\bm{v}_{i})\}_{i=1}^{r}$. If the answer is affirmative, then we can indeed take advantage of those sophisticated numerical approaches to solving (1.2) mentioned above for many practical problems. This is the focus of the present paper.

Algorithm-independent analysis

The recent work of [(8)] proposes a randomized algorithm for approximating SROA based on the power method of [(23)]. There, an error analysis specific to the proposed randomized algorithm (for the case $p=3$) shows that the decomposition $\{(\lambda_{i},\bm{v}_{i})\}_{i=1}^{r}$ of $\bm{\mathcal{T}}$ can be approximately recovered from $\widehat{\bm{\mathcal{T}}}$ in polynomial time with high probability—provided that the perturbation $\bm{\mathcal{E}}$ is sufficiently small (roughly on the order of $1/n$ under a natural measure). Our present aim is to provide a general analysis that is independent of the specific approach used to obtain rank-one approximations and it seems to be beneficial. Our analysis shows that the general SROA scheme in fact allows for perturbations to be of the order $1/\sqrt[p-1]{n}$, suggesting advantages of using more sophisticated optimization procedures and potentially more computational resources to solve each rank-one approximation step.

As motivation, we describe a simple and typical application from probabilistic modeling where perturbations of SOD tensors naturally arise. But first, we establish notations used throughout the paper, largely borrowed from [(42)].

Notation

A real $p$-th order $n$-dimensional tensor $\bm{\mathcal{A}}\in\bigotimes^{p}\mathbb{R}^{n}:=\mathbb{R}^{n\times n\times\dotsb\times n}$,

[TABLE]

is called symmetric if its entries are invariant under any permutation of their indices: for any $i_{1},i_{2},\dotsc,i_{p}\in[n]:=\left\{1,2,\dotsc,n\right\}$,

[TABLE]

for any permutation mapping $\pi$ on $[p]$.

In addition to being considered as a multi-way array, a tensor $\bm{\mathcal{A}}\in\bigotimes^{p}\mathbb{R}^{n}$ can also be interpreted as a multilinear map in the following sense: for any matrices $\bm{V}_{i}\in\mathbb{R}^{n\times m_{i}}$ for $i\in[p]$, we interpret $\bm{\mathcal{A}}(\bm{V}_{1},\bm{V}_{2},\dotsc,\bm{V}_{p})$ as a tensor in $\mathbb{R}^{m_{1}\times m_{2}\times\cdots\times m_{p}}$ whose $(i_{1},i_{2},\dotsc,i_{p})$-th entry is

[TABLE]

The following are special cases:

$p=2$ (i.e., $\bm{\mathcal{A}}$ is an $n\times n$ symmetric matrix):

[TABLE] 2. 2.

For any $i_{1},i_{2},\dotsc,i_{p}\in[n]$,

[TABLE]

where $\bm{e}_{i}$ denotes the $i$-th standard basis of $\mathbb{R}^{n}$ for any $i\in[n]$. 3. 3.

$\bm{V}_{i}=\bm{x}\in\mathbb{R}^{n}$ for all $i\in[p]$:

[TABLE]

which defines a homogeneous polynomial of degree $p$. 4. 4.

$\bm{V}_{i}=\bm{x}\in\mathbb{R}^{n}$ for all $i\in[p-1]$, and $\bm{V}_{p}=\bm{I}\in\mathbb{R}^{n\times n}$:

[TABLE]

For any tensors $\bm{\mathcal{A}}$, $\bm{\mathcal{B}}\in\bigotimes^{p}\mathbb{R}^{n}$, the inner product between them is defined as

[TABLE]

Hence, $\bm{\mathcal{A}}\bm{x}^{\otimes p}$ (defined above) can also be interpreted as the inner product between $\bm{\mathcal{A}}$ and $\bm{x}^{\otimes p}$.

Two tensor norms will be used in the paper. For a tensor $\bm{\mathcal{A}}\in\bigotimes^{p}\mathbb{R}^{n}$, its Frobenius norm is $\left\|\bm{\mathcal{A}}\right\|_{F}:=\sqrt{\left\langle\bm{\mathcal{A}},\bm{\mathcal{A}}\right\rangle}$, and its operator norm $\bm{\mathcal{A}}$, $\left\|\bm{\mathcal{A}}\right\|$, is defined as $\max_{\left\|\bm{x}_{i}\right\|=1}\bm{\mathcal{A}}(\bm{x}_{1},\bm{x}_{2},\dotsc,\bm{x}_{p})$. It is well-known that for symmetric tensors $\bm{\mathcal{A}}$, $\left\|\bm{\mathcal{A}}\right\|=\max_{\left\|\bm{x}\right\|=1}|\bm{\mathcal{A}}\bm{x}^{\otimes p}|$ (see, e.g., [(26), (27)]).

A motivating example.

To illustrate why we are particularly interested in nearly SOD tensors, we now consider the following simple probabilistic model for characterizing the topics of text documents. (We follow the description from [(8)].) Let $n$ be the number of distinct topics in the corpus, $d$ be the number of distinct words in the vocabulary, and $t\geq p$ be the number of words in each document. We identify the sets of distinct topics and words, respectively, by $[n]$ and $[d]$. The topic model posits the following generative process for a document. The document’s topic $h\in[n]$ is first randomly drawn according to a discrete probability distribution specified by $\bm{w}=(w_{1},w_{2},\dotsc,w_{n})$ (where we assume $w_{i}>0$ for each $i\in[n]$ and $\sum_{i\in[n]}w_{i}=1$):

[TABLE]

Given the topic $h$, each of the document’s $t$ words is then drawn independently from the vocabulary according to the discrete distribution specified by the probability vector $\bm{\mu}_{h}\in\mathbb{R}^{d}$; we assume that the probability vectors $\{\bm{\mu}_{i}\}_{i\in[n]}$ are linearly independent (and, in particular, $d\geq n$). The task here is to estimate these probability vectors $\bm{w}$ and $\{\bm{\mu}_{i}\}_{i\in[n]}$ based on a corpus of documents.

Denote by $\bm{M}_{2}\in\mathbb{R}^{d\times d}$ and $\bm{\mathcal{M}}_{p}\in\bigotimes^{p}\mathbb{R}^{d}$, respectively, the pairs probability matrix and $p$-tuples probability tensor, defined as follows: for all $i_{1},i_{2},\dotsc,i_{p}\in[d]$, ,

[TABLE]

It can be shown that $\bm{M}_{2}$ and $\bm{\mathcal{M}}_{p}$ can be precisely represented using $\bm{w}$ and $\{\bm{\mu}_{i}\}_{i\in[n]}$:

[TABLE]

Since $\bm{M}_{2}$ is positive semidefinite and $\mathop{\operator@font rank}\nolimits(\bm{M}_{2})=n$, $\bm{M}_{2}=\bm{U}\bm{D}\bm{U}^{\scriptscriptstyle{\top}}$ is its reduced eigenvalue decomposition. Here, $\bm{U}\in\mathbb{R}^{d\times n}$ satisfies $\bm{U}^{\scriptscriptstyle{\top}}\bm{U}=\bm{I}$, and $\bm{D}\in\mathbb{R}^{n\times n}$ is a diagonal matrix with $\mathop{\operator@font diag}\nolimits{\bm{D}}\succ\bm{0}$. Now define

[TABLE]

Then

[TABLE]

which implies that $\{\bm{v}_{1},\bm{v}_{2},\dotsc,\bm{v}_{n}\}$ are orthogonal. Moreover,

[TABLE]

Therefore, we can obtain $\{(\lambda_{i},\bm{v}_{i})\}_{i\in[n]}$ (and subsequently $\{(w_{i},\bm{\mu}_{i})\}_{i\in[n]}$ 222After obtaining $\{(\lambda_{i},\bm{v}_{i})\}_{i\in[n]}$, it is possible to obtain $\left\{(w_{i},\bm{\mu}_{i})\right\}_{i\in[n]}$ because for each $i\in[n]$, there exists $j\in[n]$ such that $w_{i}=\lambda_{j}^{2/(2-p)}$ and $\bm{\mu}_{i}=\lambda_{j}(\bm{W}^{\scriptscriptstyle{\top}})^{\dagger}\bm{v}_{j}$, where $(\bm{W}^{\scriptscriptstyle{\top}})^{\dagger}$ denotes the Moore-Penrose pseudoinverse of $\bm{W}^{\scriptscriptstyle{\top}}$.) by computing the (unique) symmetric canonical decomposition of tensor $\bm{\mathcal{T}}$, which can be perfectly achieved by SROA [(9)].

In order to obtain the tensor $\bm{\mathcal{T}}$, we need $\bm{M}_{2}$ and $\bm{\mathcal{M}}_{p}$, both of which can be estimated from a collection of documents. Due to their independence, all pairs (resp., $p$-tuples) of words in a document can be used in forming estimates of $\bm{M}_{2}$ (resp., $\bm{\mathcal{M}}_{p}$). However, the quantities $\bm{M}_{2}$ and $\bm{\mathcal{M}}_{p}$ are only known up to sampling errors, and hence, we are only able to construct a symmetric tensor $\widehat{\bm{\mathcal{T}}}$ that is, at best, only close to the one in (1.3). A critical question is whether we can still use SROA (Algorithm 1) to obtain an approximate decomposition and robustly estimate the model parameters.

Setting

Following the notation in the above example, in the sequel, we denote $\widehat{\bm{\mathcal{T}}}=\bm{\mathcal{T}}+\bm{\mathcal{E}}\in\bigotimes^{p}\mathbb{R}^{n}$. Here $\bm{\mathcal{T}}$ is a symmetric tensor that is orthogonally decomposable, i.e., $\bm{\mathcal{T}}=\sum_{i=1}^{n}\lambda_{i}\bm{v}_{i}^{\otimes p}$ with all $\lambda_{i}\neq 0$, $\left\{\bm{v}_{1},\bm{v}_{2},\dotsc,\bm{v}_{n}\right\}$ forming an orthonormal basis of $\mathbb{R}^{n}$, and $\bm{\mathcal{E}}$ is a symmetric perturbation tensor with operator norm $\varepsilon:=\left\|\bm{\mathcal{E}}\right\|$. Note that in some applications, we might instead have $\bm{\mathcal{T}}=\sum_{i=1}^{r}\lambda_{i}\bm{v}_{i}^{\otimes p}$ for some $r<n$. Our results nevertheless can be applied in that setting as well with little modification.

For simplicity, we also assume $p$ is odd and treat it as a constant in big-$O$ notations. (We discuss the even case in Section 3.4). Without loss of generality, we can assume $\lambda_{i}>0$ for all $i\in[n]$, as we can always change the sign of $\bm{v}_{i}$ to make it hold. Moreover, line 3 in Algorithm 1 simply becomes

[TABLE]

Organization

Section 2 analyzes the first iteration of Algorithm 1 and proves that $(\hat{\lambda}_{1},\hat{\bm{v}}_{1})$ is a robust estimate of a pair $(\lambda_{i},\bm{v}_{i})$ for some $i\in[n]$. A full decomposition analysis is provided in Section 3, in which we establish the following property of tensors: when $\left\|\bm{\mathcal{E}}\right\|$ is small enough, the approximation errors do not accumulate as the iteration number grows; in contrast, the use of deflation is generally not advised in the matrix setting for finding more than a handful of matrix eigenvectors due to potential instability. Numerical experiments are also reported to confirm our theoretical results.

2 Rank-One Approximation

In this section, we provide an analysis of the first step of SROA (Algorithm 1), which yield a rank-one approximation to $\widehat{\bm{\mathcal{T}}}$.

2.1 Review of Matrix Perturbation Analysis

We first state a well-known result about perturbations of the eigenvalue decomposition for symmetric matrices; this result serves as a point of comparison for our study of higher-order tensors. The result is stated just for rank-one approximations, and in a form analogous to what we are able to show for the tensor case (Theorem 2 below).

Theorem 1** ([(43), (44)]).**

Let $\bm{M}\in\mathbb{R}^{n\times n}$ be a symmetric matrix with eigenvalue decomposition $\sum_{i=1}^{n}\lambda_{i}\bm{v}_{i}\bm{v}_{i}^{\scriptscriptstyle{\top}}$, where $|\lambda_{1}|\geq|\lambda_{2}|\geq\dotsb\geq|\lambda_{n}|$ and $\{\bm{v}_{1},\bm{v}_{2},\dotsc,\bm{v}_{n}\}$ are orthonormal. Let $\widehat{\bm{M}}=\bm{M}+\bm{E}\in\mathbb{R}^{n\times n}$ for some symmetric matrix $\bm{E}$ with $\varepsilon:=\left\|\bm{E}\right\|$, and let

[TABLE]

The following holds.

• •

(Perturbation of leading eigenvalue.) $|\hat{\lambda}-\lambda_{1}|\leq\varepsilon$.

• •

(Perturbation of leading eigenvector.) Suppose $\gamma:=\min_{i\neq 1}|\lambda_{1}-\lambda_{i}|>0$. Then $\left\langle\hat{\bm{x}},\bm{v}_{1}\right\rangle^{2}\geq 1-(2\varepsilon/\gamma)^{2}$. This implies that if $2\varepsilon/\gamma\leq 1$, then $\min\{\left\|\bm{v}_{1}-\hat{\bm{x}}\right\|,\,\left\|\bm{v}_{1}+\hat{\bm{x}}\right\|\}\leq O(\varepsilon/\gamma)$.

For completeness, we give a proof of the eigenvector perturbation bound in Appendix A since it is not directly implied by results in [(44)] but essentially uses the same arguments.

2.2 Single Rank-One Approximation

The main result of this section concerns the first step of SROA (Algorithm 1) and establishes a perturbation result for nearly SOD tensors.

Theorem 2**.**

For any odd positive integer $p\geq 3$, let $\widehat{\bm{\mathcal{T}}}:=\bm{\mathcal{T}}+\bm{\mathcal{E}}\in\bigotimes^{p}\mathbb{R}^{n}$, where $\bm{\mathcal{T}}$ is a symmetric tensor with orthogonal decomposition $\bm{\mathcal{T}}=\sum_{i=1}^{n}\lambda_{i}\bm{v}_{i}^{\otimes p}$, $\left\{\bm{v}_{1},\bm{v}_{2},\dotsc,\bm{v}_{n}\right\}$ is an orthonormal basis of $\mathbb{R}^{n}$, $\lambda_{i}>0$ for all $i\in[n]$, and $\bm{\mathcal{E}}$ is a symmetric tensor with operator norm $\varepsilon:=\left\|\bm{\mathcal{E}}\right\|$. Let $\hat{\bm{x}}\in\arg\max_{\left\|\bm{x}\right\|_{2}=1}\widehat{\bm{\mathcal{T}}}\bm{x}^{\otimes p}$ and $\hat{\lambda}:=\widehat{\bm{\mathcal{T}}}\hat{\bm{x}}^{\otimes p}$. Then there exists $j\in[n]$ such that

[TABLE]

To prove Theorem 2, we first establish an intermediate result. Let $x_{i}:=\left\langle\bm{v}_{i},\hat{\bm{x}}\right\rangle$, so $\hat{\bm{x}}=\sum_{i=1}^{n}x_{i}\bm{v}_{i}$ and $\sum_{i=1}^{n}x_{i}^{2}=1$ since $\{\bm{v}_{1},\bm{v}_{2},\dotsc,\bm{v}_{n}\}$ is an orthonormal basis for $\mathbb{R}^{n}$ and $\left\|\hat{\bm{x}}\right\|=1$. We reorder the indicies $[n]$ so that

[TABLE]

Our intermediate result, derived by simply bounding $\hat{\lambda}$ from both above and below, is as follows.

Lemma 3**.**

In the notation from above,

[TABLE]

where $\lambda_{\max}=\max_{i\in[n]}\lambda_{i}$.

Proof.

To show (2.3), we will bound $\hat{\lambda}=\widehat{\bm{\mathcal{T}}}\hat{\bm{x}}^{\otimes p}$ from both above and below.

For the upper bound, we have

[TABLE]

where the last inequality follows since $\sum_{i=1}^{n}x_{i}^{2}=1$.

On the other hand,

[TABLE]

Combining (2.4) and (2.5), it can be easily verified that

[TABLE]

and moreover,

[TABLE]

which implies that $x_{1}^{p-1}=|x_{1}|^{p-2}\cdot|x_{1}|\geq 1-4\varepsilon/\lambda_{1}.$ ∎

Remark 1**.**

The higher-order requirement, $p\geq 3$, is crucial in the analysis. Specifically we can bound $|x_{1}|$ below by the lower bound of $|x_{1}|^{p-2}$, which can be done by bounding $\hat{\lambda}=\widehat{\bm{\mathcal{T}}}\hat{\bm{x}}^{\otimes p}$ from both above and below. This essentially explains why Lemma 3, different from the matrix case ($p=2$), does not rely on the spectral gap condition.

The bound $|\hat{\lambda}-\lambda_{1}|\leq\varepsilon$ proved in Lemma 3 is comparable to the matrix counterpart in Theorem 1, and is optimal in the worst case. Consider ${\bm{\mathcal{T}}}=\sum_{i=1}^{n}\lambda_{i}\bm{e}_{i}^{\otimes p}$ with $\lambda_{1}>\lambda_{2}\geq\ldots\geq\lambda_{n}\geq 0$ and $\bm{\mathcal{E}}=\varepsilon\bm{e}_{1}^{\otimes p}$, for some $\varepsilon>0$. Then clearly $\hat{\lambda}=\lambda_{1}+\varepsilon$ and $|\hat{\lambda}-\lambda_{1}|=\varepsilon$.

Moreover, when $\bm{\mathcal{E}}$ vanishes, Lemma 3 leads directly to the following result given in [(9)].

Corollary 4** ([(9)]).**

Suppose $\bm{\mathcal{E}}=\bm{0}$ (i.e., $\widehat{\bm{\mathcal{T}}}=\bm{\mathcal{T}}=\sum_{i=1}^{n}\lambda_{i}\bm{v}_{i}^{\otimes p}$ is orthogonally decomposable). Then Algorithm 1 computes $\{(\lambda_{i},\bm{v}_{i})\}_{i=1}^{n}$ exactly.

However, compared to Theorem 1, the bound for $\bm{x}$ is $|x_{1}|\geq 1-2\varepsilon/\lambda_{1}$ appears to be suboptimal; this is because the bound only implies $\left\|\hat{\bm{x}}-\bm{v}_{1}\right\|=O(\sqrt{\varepsilon/\lambda_{1}})$. In the following, we will proceed to improve this result to $\left\|\hat{\bm{x}}-\bm{v}_{1}\right\|=O({\varepsilon/\lambda_{1}})$ by using the first-order optimality condition [(22)]. See [(42)] for a discussion in the present context.

Consider the Lagrangian function corresponding to the optimization problem $\max_{\left\|\bm{x}\right\|=1}\widehat{\bm{\mathcal{T}}}\bm{x}^{\otimes p},$

[TABLE]

where $\lambda\in\mathbb{R}$ corresponds to the Lagrange multiplier for the equality constraint. As $\hat{\bm{x}}\in\arg\max_{\left\|\bm{x}\right\|_{2}=1}\widehat{\bm{\mathcal{T}}}\bm{x}^{\otimes p}$ (and the linear independent constraint qualification [(22)] can be easily verified), there exists $\bar{\lambda}\in\mathbb{R}$ such that

[TABLE]

Moreover, as $\left\|\hat{\bm{x}}\right\|=1$, $\bar{\lambda}=\bar{\lambda}\left\langle\bm{x},\bm{x}\right\rangle=\widehat{\bm{\mathcal{T}}}\bm{x}^{\otimes p}=\hat{\lambda}$. Thus we have $\hat{\lambda}\hat{\bm{x}}=\widehat{\bm{\mathcal{T}}}\hat{\bm{x}}^{\otimes p-1}$.

We are now ready to prove Theorem 2.

Proof of Theorem 2.

The first inequality in (2.1) has been proved in Lemma 3, so we are left to prove the second one.

From the first-order optimality condition above, we have

[TABLE]

Thus,

[TABLE]

by the triangle inequality. By Lemma 3, we have

[TABLE]

Moreover,

[TABLE]

where we have used Lemma 3 and the fact that $\lambda_{2}/\lambda_{1}\leq\lambda_{\max}/\lambda_{1}\leq 1+2\varepsilon/\lambda_{1}.$ Substituting (2.9) and (2.2) back into (2.8), we can easily obtain

[TABLE]

∎

Remark 2**.**

When $p>3$, we can slightly sharpen (2.11) to

[TABLE]

by replacing (2.2) with

[TABLE]

Theorem 2 indicates that the first step of SROA for a nearly SOD tensor approximately recovers $(\lambda_{j},\bm{v}_{j})$ for some $j\in[n]$. In particular, whenever $\varepsilon$ is small enough relative to $\lambda_{1}$ (e.g., $\varepsilon\leq\lambda_{1}/2$), there always exists $j\in[n]$ such that $|\hat{\lambda}-\lambda_{j}|\leq\varepsilon$ and $\left\|\hat{\bm{x}}-\bm{v}_{j}\right\|_{2}\leq 10\cdot(1+1/2)\varepsilon/\lambda_{j}=15\varepsilon/\lambda_{j}$. This is analogous to Theorem 1, except that the spectral gap condition required in Theorem 1 is not necessary at all for the perturbation bounds of SOD tensors.

2.3 Numerical Verifications for Theorem 2

We generate nearly symmetric orthogonally decomposable tensors $\widehat{\bm{\mathcal{T}}}=\bm{\mathcal{T}}+\bm{\mathcal{E}}\in\mathbb{R}^{10\times 10\times 10}$ in the following manner. We let the underlying symmetric orthogonally decomposable tensor $\bm{\mathcal{T}}$ be the diagonal tensor with all diagonal entries equal to one, i.e., $\bm{\mathcal{T}}=\sum_{i=1}^{10}\bm{e}_{i}^{\otimes 3}$ (where $\bm{e}_{i}$ is the $i$-th coordinate basis vector). The perturbation tensor $\bm{\mathcal{E}}$ is generated under the following three random models:

Binary:

independent entries $\mathcal{E}_{i,j,k}\in\{\pm\sigma\}$ uniformly at random;

Uniform:

indepedent entries $\mathcal{E}_{i,j,k}\in[-2\sigma,2\sigma]$ uniformly at random;

Gaussian:

independent entries $\mathcal{E}_{i,j,k}\sim\mathcal{N}(0,\sigma^{2})$;

where $\sigma$ is varied from $0.0001$ to $0.2$ with increment $0.0001$, and one instance is generated for each value of $\sigma$.

For every randomly generated instance, we solve the polynomial optimization problems

[TABLE]

using the general polynomial solver GloptiPoly 3 (henrion2009gloptipoly, ), and set $\hat{\lambda}_{1}:=\widehat{\bm{\mathcal{T}}}\hat{\bm{v}}_{1}^{\otimes 3}.$

In Figure 1, we plot the approximation error $|\hat{\lambda}_{1}-1|$ and $\min_{j\in[10]}\left\|\hat{\bm{v}}_{1}-\bm{e}_{j}\right\|$, respectively against the value of norm $\left\|\bm{\mathcal{E}}\right\|$. Each (red) point corresponds to one randomly generated instance, and the (blue) lines are the upper bounds given in Theorem 2. We observe no instance violating the theoretical bounds.

3 Full Decomposition Analysis

In the second iteration of Algorithm 1, we have

[TABLE]

where, for some $j\in[n]$,

[TABLE]

Theorem 2 can be directly applied again by bounding the error norm $\|\widehat{\bm{\mathcal{E}}}\|$. However, since

[TABLE]

it appears that the approximation error may increase dramatically with the iteration number.

Fortunately, a more careful analysis shows that approximation error does not in fact accmulate in this way. The high-level reason is that while the operator norm $\|\lambda_{j}\bm{v}_{j}^{\otimes p}-\hat{\lambda}_{1}\hat{\bm{v}}_{1}^{\otimes p}\|$ is of order $\varepsilon$, the relevant quantity is essentially $\left(\lambda_{j}\bm{v}_{j}^{\otimes p}-\hat{\lambda}_{1}\hat{\bm{v}}_{1}^{\otimes p}\right)$ operating on the direction of $\hat{\bm{v}}_{2}$, i.e. $|(\lambda_{j}\bm{v}_{j}^{\otimes p}-\hat{\lambda}_{1}\hat{\bm{v}}_{1}^{\otimes p})\hat{\bm{v}}_{2}^{\otimes p}|$, which only gives rise to a quantity of order $\varepsilon^{2}$ because $p\geq 3$. This enables us to keep the approximation errors under control.

The main result of this section is as follows.

Theorem 5**.**

Pick any odd positive integer $p\geq 3$. There exists a positive constant $c_{0}=c_{0}(p)>0$ such that the following holds. Let $\widehat{\bm{\mathcal{T}}}:=\bm{\mathcal{T}}+\bm{\mathcal{E}}\in\bigotimes^{p}\mathbb{R}^{n}$, where $\bm{\mathcal{T}}$ is a symmetric tensor with orthogonal decomposition $\bm{\mathcal{T}}=\sum_{i=1}^{n}\lambda_{i}\bm{v}_{i}^{\otimes p}$, $\left\{\bm{v}_{1},\bm{v}_{2},\dotsc,\bm{v}_{n}\right\}$ is an orthonormal basis of $\mathbb{R}^{n}$, $\lambda_{i}>0$ for all $i\in[n]$, and $\bm{\mathcal{E}}$ is a symmetric tensor with operator norm $\varepsilon:=\left\|\bm{\mathcal{E}}\right\|$. Assume $\varepsilon\leq c_{0}\lambda_{\min}/n^{1/(p-1)}$, where $\lambda_{\min}:=\min_{i\in[n]}\lambda_{i}$. Let $\{(\hat{\lambda}_{i},\hat{\bm{v}}_{i})\}_{i\in[n]}$ be the output of Algorithm 1 for input $\widehat{\bm{\mathcal{T}}}$. Then there exists a permutation $\pi$ on $[n]$ such that

[TABLE]

3.1 Deflation Analysis

The proof of Theorem 5 is based on the following lemma, which provides a careful analysis of the errors introduced in $\bm{\mathcal{T}}_{i}$ from steps $1,2,\dotsc,i$ in Algorithm 1. This lemma is a generalization of a result from JMLR:v15:anandkumar14b (which only dealt with the $p=3$ case) and also more transparently reveals the sources of errors that result from deflation.

Lemma 6**.**

Fix a subset $S\subseteq[n]$ and assume that $0\leq\hat{\varepsilon}\leq\lambda_{i}/2$ for each $i\in S$. Choose any $\{(\hat{\lambda}_{i},\hat{\bm{v}}_{i})\}_{i\in S}\subset\mathbb{R}\times\mathbb{R}^{n}$ such that

[TABLE]

and define $\bm{\Delta}_{i}:=\lambda_{i}\bm{v}_{i}^{\otimes p}-\hat{\lambda}_{i}\hat{\bm{v}}_{i}^{\otimes p}$ for $i\in S$. Pick any unit vector $\bm{x}=\sum_{i=1}^{n}x_{i}\bm{v}_{i}$. Let $S_{1}\subseteq S$ be the indices $i\in[n]$ such that $\lambda_{i}|x_{i}|\geq 4\hat{\varepsilon}$, and let $S_{2}:=S\setminus S_{1}$. Then

[TABLE]

These imply that there exists positive constants $c_{1},c_{2}>0$, depending only on $p$, such that

[TABLE]

Remark 3**.**

Lemma 6 indicates that the accumulating error $\sum_{i\in S}\bm{\Delta}_{i}$ much less severely affects vectors that are incoherent with $\left\{\bm{v}_{i}:i\in S\right\}$. For instance, $\left\|\sum_{i\in S}\bm{\Delta}_{i}\bm{v}_{i}^{\otimes p-1}\right\|=O(\hat{\varepsilon}^{2})$ for $i\in[n]\setminus S$, while $\left\|\sum_{i\in S}\bm{\Delta}_{i}\bm{v}_{i}^{\otimes p-1}\right\|=O(\hat{\varepsilon})$ for $i\in S$.

3.2 Proof of Main Theorem

We now use Lemma 6 to prove the main theorem.

Proof of Theorem 5.

It suffices to prove that the following property holds for each $i\in[n]$:

[TABLE]

The proof is by induction. The base case of ($*$ ‣ 3.2) (where $i=1$) follows directly from by Theorem 2.

Assume the induction hypothesis  ($*$ ‣ 3.2) is true for some $i\in[n-1]$. We will prove that there exists an $l\in[n]\setminus\{\pi(j):j\in[i]\}$ that satisfies

[TABLE]

To simplify notation, we assume without loss of generality (by renumbering indices) that $\pi(j)=j$ for each $j\in[i]$. Let $\hat{\bm{x}}=\sum_{i\in[n]}x_{i}\bm{v}_{i}:=\hat{\bm{v}}_{i+1}$ and $\hat{\lambda}:=\hat{\lambda}_{i+1}$, and further assume without loss of generality (again by renumbering indices) that

[TABLE]

In the following, we will show that $l=i+1$ is an index satisfying (3.6). We use the assumption that

[TABLE]

(which holds with a suitable choice of $c_{0}$ in the theorem statement). Here, $c_{1}$ and $c_{2}$ are the constants from Lemma 6 when $\hat{\varepsilon}=10\varepsilon$. It can be verified that ($*$ ‣ 3.2) implies that the conditions for Lemma 6 are satisfied with this value of $\hat{\varepsilon}$.

Recall that $\hat{\lambda}=\widehat{\bm{\mathcal{T}}}_{i}\hat{\bm{x}}^{\otimes p}$, where

[TABLE]

We now bound $\hat{\lambda}$ from above and below. For the lower bound, we use (3.4) from Lemma 6 to obtain

[TABLE]

where $\lambda_{\max,i}:=\max_{j\in[n]\setminus[i]}\lambda_{j}$ and $\lambda_{\min}:=\min_{j\in[n]}\lambda_{j}$; the final inequality uses the conditions on $\varepsilon$ in (3.7). For the upper bound, we have

[TABLE]

The first inequality above follows from (3.5) in Lemma 6; the third inequality uses the fact that $\sum_{j=1}^{n}x_{j}^{2}=1$ as well as the conditions on $\varepsilon$ in (3.7). If the $\max$ is achieved by the second argument $10c_{1}\varepsilon$, then combining (3.8) and (3.9) gives

[TABLE]

a contradiction of (3.7). Therefore the $\max$ in (3.9) must be achieved by $\lambda_{i+1}|x_{i+1}|^{p-2}$, and hence combining (3.8) and (3.9) gives

[TABLE]

This in turn implies that

[TABLE]

Thus, we have shown that $\hat{\bm{x}}$ is indeed coherent with $\hat{\bm{v}}_{i+1}$. Next, we will sharpen the bound for $\left\|\hat{\bm{x}}-\hat{\bm{v}}_{i+1}\right\|$ by considering the first order optimality condition.

Since $\hat{\bm{x}}\in\arg\min_{\left\|\bm{x}\right\|_{2}=1}\widehat{\bm{\mathcal{T}}}_{i}\bm{x}^{\otimes p}$, a first-order optimality condition similar to (2.7) implies $\hat{\lambda}=\widehat{\bm{\mathcal{T}}}_{i}\hat{\bm{x}}^{\otimes p}$. Thus

[TABLE]

Therefore

[TABLE]

For the third term in (3.11), we use the fact that $|x_{i+2}|\leq\sqrt{1-x_{i+1}^{2}}$, the bounds from (3.10) and the conditions on $\varepsilon$ in (3.7) to obtain

[TABLE]

For the last term in (3.11), we use (3.5) from Lemma 6 and the conditions on $\varepsilon$ in (3.7) to get

[TABLE]

Therefore, substituting (3.10), (3.13) and $\left\|\bm{\mathcal{E}}\right\|\leq\varepsilon$ into (3.11) gives

[TABLE]

∎

3.3 Stability of Full Decomposition

Theorem 5 states a (perhaps unexpected) phenomenon that the approximation errors do not accumulate with iteration number, whenever the perturbation error is small enough. In this subsection, we numerically corroborate this fact.

We generate nearly symmetric orthogonally decomposable tensors $\widehat{\bm{\mathcal{T}}}=\bm{\mathcal{T}}+\bm{\mathcal{E}}\in\mathbb{R}^{10\times 10\times 10}$ as follows. We construct the underlying symmetric orthogonally decomposable tensor $\bm{\mathcal{T}}$ as the diagonal tensor with all diagonal entries equal to one, i.e., $\bm{\mathcal{T}}=\sum_{i=1}^{10}\bm{e}_{i}^{\otimes 3}$ (where $\bm{e}_{i}$ is the $i$-th coordinate basis vector). The perturbation tensor $\bm{\mathcal{E}}$ is generated under three random models with $\sigma=0.01$:

Binary:

independent entries $\mathcal{E}_{i,j,k}\in\{\pm\sigma\}$ uniformly at random;

Uniform:

indepedent entries $\mathcal{E}_{i,j,k}\in[-2\sigma,2\sigma]$ uniformly at random;

Gaussian:

independent entries $\mathcal{E}_{i,j,k}\sim\mathcal{N}(0,\sigma^{2})$.

For each random model, we generate $500$ random instances, and apply Algorithm 1 to each $\widehat{\bm{\mathcal{T}}}$ to obtain approximate pairs $\{(\hat{\lambda}_{i},\hat{\bm{v}}_{i})\}_{i\in[10]}$. Again, we use GloptiPoly 3 to solve the polynomial optimization problem in Algorithm 1.

In Figure 2, we plot the mean and the standard deviation of the approximation errors for $\hat{\lambda}_{i}$ and $\hat{\bm{v}}_{i}$ from the $500$ random instances (for each $i\in[10]$). These indeed do not appear to grow or accumulate as the iteration number increases. This is consistent with our results in Theorem 5.

3.4 When $p$ is Even

We now briefly discuss the case where the order of the tensor is even, i.e., $p\geq 4$ is an even integer.

Let $\widehat{\bm{\mathcal{T}}}:=\bm{\mathcal{T}}+\bm{\mathcal{E}}\in\bigotimes^{p}\mathbb{R}^{n}$, where $\bm{\mathcal{T}}$ is a symmetric tensor with orthogonal decomposition $\bm{\mathcal{T}}=\sum_{i=1}^{n}\lambda_{i}\bm{v}_{i}^{\otimes p}$, where $\left\{\bm{v}_{1},\bm{v}_{2},\dotsc,\bm{v}_{n}\right\}$ is an orthonormal basis of $\mathbb{R}^{n}$, $\lambda_{i}\neq 0$ for all $i\in[n]$, and $\bm{\mathcal{E}}$ is a symmetric tensor with operator norm $\varepsilon:=\left\|\bm{\mathcal{E}}\right\|$. Note that unlike the case when $p$ is odd, we cannot assume $\lambda_{i}>0$ for all $i\in[n]$, and correspondingly, line 3 in Algorithm 1 now becomes

[TABLE]

Nevertheless, the pair $(\hat{\lambda}_{i},\hat{\bm{v}}_{i})$ still satisfies the first-order optimality condition $\hat{\lambda}_{i}\hat{\bm{v}}_{i}=\widehat{\bm{\mathcal{T}}}\hat{\bm{v}}_{i}^{\otimes p-1}$.

Our proof for Theorem 5 can be easily modified and leads to the following result: there exists a positive constant $\hat{c}_{0}=\hat{c}_{0}(p)>0$ such that whenever $\varepsilon\leq\hat{c}_{0}\left(\min_{i\in[n]}|\lambda_{i}|\right)/n^{1/(p-1)}$, there exists a permutation $\pi$ on $[n]$ such that

[TABLE]

4 Conclusion

This paper sheds light on a problem at the intersection of numerical linear algebra and statistical estimation, and our results draw upon and enrich the literature in both areas.

From the perspective of numerical linear algebra, SROA was previously only known to exactly recover the symmetric canonical decomposition of an orthogonal decomposable tensor. Our results show that it can robustly recover (approximate) orthogonal decompositions even when applied to nearly SOD tensors; this substantially enlarges the applicability of SROA.

Previous work on statistical estimation via orthogonal tensor decompositions considered the specific randomized power iteration algorithm of JMLR:v15:anandkumar14b , which has been successfully applied in a number of contexts chaganty2013linear ; ZHPA13-contrast ; azar2013bandit ; huang2013hmm ; AGHK14-community ; doshi2014graph . Our results provide formal justification for using other rank-one approximation methods in these contexts, and it seems to be quite beneficial, in terms of sample complexity and statistical efficiency, to use more sophisticated methods. Specifically, the perturbation error $\left\|\bm{\mathcal{E}}\right\|$ that can be tolerated is relaxed from power iteration’s $O(1/n)$ to $O(1/\sqrt[p-1]{n})$. In future work, we plan to empirically investigate these potential benefits in a number of applications.

We also note that solvers for rank-one tensor approximation often lack rigorous runtime or error analyses, which is not surprising given the computational difficulty of the problem for general tensors (hillar2009most, ). However, tensors that arise in applications are often more structured, such as being nearly SOD. Thus, another promising future research direction is to sidestep computational hardness barriers by developing and analyzing methods for such specially structured tensors (see also JMLR:v15:anandkumar14b ; barak2014dictionary for ideas along this line).

Acknowledgements

Daniel Hsu acknowledges support from a Yahoo ACE award. Donald Goldfarb acknowledges support from NSF Grants DMS-1016571 and CCF-1527809.

Appendix A Proof of Theorem 1

Since $\widehat{\bm{M}}$ is symmetric, it has an eigenvalue decomposition $\sum_{i=1}^{n}\hat{\lambda}_{i}\hat{\bm{v}}_{i}\hat{\bm{v}}_{i}^{\scriptscriptstyle{\top}}$, where $|\hat{\lambda}_{1}|\geq|\hat{\lambda}_{2}|\geq\dotsb\geq|\hat{\lambda}_{n}|$ and $\{\hat{\bm{v}}_{1},\hat{\bm{v}}_{2},\dotsc,\hat{\bm{v}}_{n}\}$ are orthonormal. It is straightforward to obtain:

[TABLE]

By Weyl’s inequality weyl1912asymptotische ,

[TABLE]

To bound $\left\langle\hat{\bm{x}},\bm{v}_{1}\right\rangle^{2}$, we employ an argument very similar to one from davis1970rotation . Observe that

[TABLE]

Moreover,

[TABLE]

and therefore

[TABLE]

Combining the upper and lower bounds on $\left\|\bm{M}\hat{\bm{x}}-\lambda_{1}\hat{\bm{x}}\right\|^{2}$ gives $\left\langle\hat{\bm{x}},\bm{v}_{1}\right\rangle^{2}\geq 1-(2\varepsilon/\gamma)^{2}$ as claimed. ∎

Appendix B Proof of Lemma 6

Proof.

The lemma holds trivially if $\hat{\varepsilon}=0$. So we may assume $\hat{\varepsilon}>0$. Therefore, for every $i\in S_{1}$, we have $|x_{i}|\geq 4\hat{\varepsilon}/\lambda_{i}>0$. Let $c_{i}:=\left\langle\bm{v}_{i},\hat{\bm{v}}_{i}\right\rangle$, $\bm{w}_{i}:=(\hat{\bm{v}}_{i}-c_{i}\bm{v}_{i})/\left\|\hat{\bm{v}}_{i}-c_{i}\bm{v}_{i}\right\|_{2}$, and $y_{i}:=\left\langle\bm{w}_{i},\bm{x}\right\rangle$, so

[TABLE]

We first establish a few inequalities that will be frequently used later. Since $|\lambda_{i}-\hat{\lambda}_{i}|\leq\hat{\varepsilon}\leq\lambda_{i}/2$, one has $\hat{\varepsilon}/\lambda_{i}\leq 1/2$, and $1/2\leq\hat{\lambda}_{i}/\lambda_{i}\leq 3/2$. Also, since $c_{i}\geq 1-2(\hat{\varepsilon}/\lambda_{i})^{2}\geq 1/2$,

[TABLE]

For each $i\in S$,

[TABLE]

Therefore, due to the orthonormality of $\{\bm{v}_{i}\}_{i\in[n]}$ and the triangle inequality, for each $j\in\{1,2\}$,

[TABLE]

We now prove (3.1). For any $i\in S_{1}$, since $x_{i}\neq 0$, we may write (B.1) as

[TABLE]

Observe that

[TABLE]

because $|c_{i}|\geq 1/2$ and $\sqrt{1-c_{i}^{2}}\geq 2\hat{\varepsilon}/\lambda_{i}$. Moreover, since $1+(p-1)z\leq(1+z)^{p-1}\leq 1+(2^{p-1}-1)z$ for any $z\in[0,1]$,

[TABLE]

Therefore,

[TABLE]

The second inequality above is obtained using the inequality $|(1+a)(1+b)-1|\leq|a|+|b|+|ab|$ for any $a,b\in\mathbb{R}$, together with the inequality from (B.3) and the fact $|1-c_{i}^{p}|\leq 2p(\hat{\varepsilon}/\lambda_{i})^{2}\leq p(\hat{\varepsilon}/\lambda_{i})$. Using the resulting inequality in (B.4), the first summand in (B.2) can be bounded as

[TABLE]

To bound the second summand in (B.2), we have

[TABLE]

The second inequality uses the facts $c_{i}\geq 1/2$ and $\sqrt{1-c_{i}^{2}}\leq 2\hat{\varepsilon}/\lambda_{i}$; the last inequality uses the facts $\hat{\lambda}_{i}/\lambda_{i}\leq 3/2$ and $\lambda_{i}|x_{i}|\geq 4\hat{\varepsilon}$. Combining (B.5) and (B.6) gives the claimed inequality in (3.1) via (B.2).

It remains to prove (3.2). For each $i\in S_{2}$,

[TABLE]

The second inequality uses the facts $\sqrt{1-c_{i}^{2}}\leq 2\hat{\varepsilon}/\lambda_{i}$ and $\lambda_{i}|x_{i}|<4\hat{\varepsilon}$ for all $i\in S_{2}$; the third inequality uses the fact $\hat{\lambda}_{i}/\lambda_{i}\leq 3/2$. Therefore

[TABLE]

Moreover,

[TABLE]

Combining (B.7) and (B.8) establishes (3.2) via (B.1) (with $j=2$). ∎