Greedy Approaches to Symmetric Orthogonal Tensor Decomposition

TL;DR

This paper reviews and compares perturbation bounds for incremental rank-one approximation methods in symmetric orthogonal tensor decomposition, highlighting their theoretical properties and practical implications.

Contribution

It establishes and compares perturbation bounds for two natural incremental approaches to symmetric orthogonal tensor decomposition.

Findings

Perturbation bounds for two incremental methods are derived and analyzed.

Numerical experiments illustrate the theoretical results.

Open questions for future research are discussed.

Abstract

Finding the symmetric and orthogonal decomposition (SOD) of a tensor is a recurring problem in signal processing, machine learning and statistics. In this paper, we review, establish and compare the perturbation bounds for two natural types of incremental rank-one approximation approaches. Numerical experiments and open questions are also presented and discussed.

Equations

$$\mathcal{T}={\lambda }_1{\mathbf{v}}_1^{\otimes p}+{\lambda }_2{\mathbf{v}}_2^{\otimes p}+\cdots +{\lambda }_n{\mathbf{v}}_n^{\otimes p},$$

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

$$({\hat{\lambda }}_k,{\hat{\mathbf{v}}}_k)\in$$

$$\text{mbox}s.t.\Vert \mathbf{v}\Vert =1$$

$$-\theta \le ⟨\mathbf{v},{\hat{\mathbf{v}}}_i⟩\le \theta ,i=0,1,2,\dots ,k-1$$

$$\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}(V_1{)}_{j_1{i}_1}(V_2{)}_{j_2{i}_2}\cdots (V_p{)}_{j_p{i}_p}.$$

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

$$⟨\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}.$$

$$\min \left\{|{\lambda }_{\pi (j)}-{\hat{\lambda }}_j|,|{\lambda }_{\pi (j)}+{\hat{\lambda }}_j|\right\}=0,$$

$$\min \left\{\Vert {\mathbf{v}}_{\pi (j)}-{\hat{\mathbf{v}}}_j\Vert ,\Vert {\mathbf{v}}_{\pi (j)}+{\hat{\mathbf{v}}}_j\Vert \right\}=0.$$

$$\min \left\{|{\lambda }_{\pi (j)}-{\hat{\lambda }}_j|,|{\lambda }_{\pi (j)}+{\hat{\lambda }}_j|\right\}\le 2\epsilon ,$$

$$\min \left\{\Vert {\mathbf{v}}_{\pi (j)}-{\hat{\mathbf{v}}}_j\Vert ,\Vert {\mathbf{v}}_{\pi (j)}+{\hat{\mathbf{v}}}_j\Vert \right\}\le 20\epsilon /\left|{\lambda }_{\pi (j)}\right|.$$

$$\min \left\{|{\lambda }_{\pi (j)}-{\hat{\lambda }}_j|,|{\lambda }_{\pi (j)}+{\hat{\lambda }}_j|\right\}\le \epsilon ,$$

$$\min \left\{\Vert {\mathbf{v}}_{\pi (j)}-{\hat{\mathbf{v}}}_j\Vert ,\Vert {\mathbf{v}}_{\pi (j)}+{\hat{\mathbf{v}}}_j\Vert \right\}\le (6.2+4\kappa )\epsilon /|{\lambda }_{\pi (j)}|.$$

$$\min \left\{|{\lambda }_j-\hat{\lambda }|,|{\lambda }_j+\hat{\lambda }|\right\}\le \epsilon ,\text{mbox}and$$

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

$${\hat{\mathbf{v}}}_k\in \underset{\mathbf{v}\in {\mathbb{R}}^n}{\mathrm{argmax}}\hat{\mathcal{T}}{\mathbf{v}}^{\otimes p}\text{mbox}s.t.\Vert \mathbf{v}\Vert =1,\text{mbox}and|⟨\mathbf{v},{\hat{\mathbf{v}}}_i⟩|\le \theta \forall i\in [k-1],$$

$$\ast$$

$$|{\hat{\lambda }}_1-{\lambda }_j|\le \epsilon ,\text{mbox}and\Vert {\hat{\mathbf{v}}}_1-{\mathbf{v}}_j\Vert \le 10\frac{\epsilon }{{\lambda }_j}\left(1+\frac{\epsilon }{{\lambda }_j}\right)\le \frac{10.2\epsilon }{{\lambda }_j}\le (6.2+4\kappa )\frac{\epsilon }{{\lambda }_j},$$

$$|{\hat{\lambda }}_{k+1}-{\lambda }_l|\le \epsilon ,\text{mbox}and\Vert {\hat{\mathbf{v}}}_{k+1}-{\mathbf{v}}_l\Vert \le \frac{(6.2+4\kappa )\epsilon }{{\lambda }_{\pi (l)}}.$$

$$\hat{\mathbf{x}}\in \arg \underset{\mathbf{v}\in {\mathbb{R}}^n}{\max }\hat{\mathcal{T}}{\mathbf{v}}^{\otimes p}\text{mbox}s.t.\Vert \mathbf{v}\Vert =1,|⟨{\hat{\mathbf{v}}}_i,\mathbf{v}⟩|\le \theta \forall i\in [k],$$

$${\lambda }_1|x_1{|}^{p-2}\ge {\lambda }_2|x_2{|}^{p-2}\ge \dots \ge {\lambda }_k|x_k{|}^{p-2},\text{mbox}and$$

$${\lambda }_{k+1}|x_{k+1}{|}^{p-2}\ge {\lambda }_{k+2}|x_{k+2}{|}^{p-2}\ge \dots \ge {\lambda }_n|x_n{|}^{p-2}.$$

$$|⟨{\mathbf{v}}_i,{\hat{\mathbf{v}}}_j⟩|$$

$$\le \frac{(6.2+4\kappa )\epsilon }{{\lambda }_{\pi (j)}}\le \frac{(6.2+4\kappa ){\theta }^2\cancel{{\lambda }_{\min }}}{12.5\cancel{{\lambda }_{\min }}}=\frac{(6.2+4\kappa )\theta }{12.5}\cdot \theta =\frac{6.2+4\kappa }{25\kappa }\theta <\theta ,$$

$$\hat{\lambda }=\hat{\mathcal{T}}{\hat{\mathbf{x}}}^{\otimes p}\ge \underset{i\in [n]\backslash [k]}{\max }\hat{\mathcal{T}}{\mathbf{v}}_i^{\otimes p}\ge \underset{i\in [n]\backslash [k]}{\max }{\lambda }_i-\epsilon \ge {\lambda }_{k+1}-\epsilon .$$

$$\hat{\lambda }=\hat{\mathcal{T}}{\hat{\mathbf{x}}}^{\otimes p}=(\mathcal{T}+\mathcal{E}){\hat{\mathbf{x}}}^{\otimes p}$$

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

$$\le \sum _{i=1}^k{\lambda }_i|x_i{|}^{p-2}{x}_i^2+\sum _{i=k+1}^n{\lambda }_i|x_i{|}^{p-2}{x}_i^2+\mathcal{E}{\hat{\mathbf{x}}}^{\otimes p}$$

$$\le \max \{{\lambda }_1|x_1{|}^{p-2},{\lambda }_{k+1}|x_{k+1}{|}^{p-2}\}+\epsilon ,$$

$${\lambda }_{k+1}-\epsilon \le \underset{i\in [n]\backslash [k]}{\max }{\lambda }_i-\epsilon \le \hat{\lambda }\le \max \{{\lambda }_1|x_1{|}^{p-2},{\lambda }_{k+1}|x_{k+1}{|}^{p-2}\}+\epsilon .$$

$${\lambda }_1|x_1{|}^{p-2}+\epsilon$$

$$\le {\lambda }_1|x_1|+\epsilon ={\lambda }_1|⟨\hat{\mathbf{x}},{\mathbf{v}}_1⟩|+\epsilon ={\lambda }_1|⟨\hat{\mathbf{x}},{\hat{\mathbf{v}}}_{{\pi }^{-1}(1)}⟩+⟨\hat{\mathbf{x}},{\mathbf{v}}_1-{\hat{\mathbf{v}}}_{{\pi }^{-1}(1)}⟩|+\epsilon$$

$$\le {\lambda }_1|⟨\hat{\mathbf{x}},{\hat{\mathbf{v}}}_{{\pi }^{-1}(1)}⟩|+{\lambda }_1\Vert {\mathbf{v}}_1-{\hat{\mathbf{v}}}_{{\pi }^{-1}(1)}\Vert +\epsilon \le {\lambda }_1\theta +\cancel{{\lambda }_1}\frac{(6.2+4\kappa )\epsilon }{\cancel{{\lambda }_1}}+\epsilon$$

$$\le \frac{{\lambda }_1}{2\kappa }+6.2\epsilon +4\kappa \epsilon +\epsilon \le \frac{{\lambda }_{\min }}{2}+\frac{{\lambda }_{\min }}{12.5}+7.2\epsilon <{\lambda }_{\min }-\epsilon$$

Full text

\slugger

simaxxxxxxxxx–x

1 Introduction

A $p$-way $n$-dimensional tensor $\bm{\mathcal{T}}$, namely $\bm{\mathcal{T}}\in\bigotimes^{p}\mathbb{R}^{n}:=\mathbb{R}^{n\times n\times\dotsb\times n}$, is called symmetrically orthogonally decomposable (SOD) [(1), (2)] (a.k.a. odeco in [(3)]) if it can be expressed as a linear combination over the real field of symmetric $p$-th powers of $n$ vectors that generate an orthonormal basis of $\mathbb{R}^{n}$. Mathematically, $\bm{\mathcal{T}}$ is SOD if there exist $\bm{\lambda}=[\lambda_{1},\lambda_{2},\ldots,\lambda_{n}]\in\mathbb{R}^{n}$ and an orthogonal matrix $\bm{V}=[\bm{v}_{1},\bm{v}_{2},\ldots,\bm{v}_{n}]\in\mathbb{R}^{n\times n}$ such that

[TABLE]

where $\bm{v}^{\otimes p}$, the symmetric $p$-th power of the vector $\bm{v}$, denotes a $p$-way $n$-dimensional tensor with $(\bm{v}^{\otimes p})_{i_{1}i_{2}\cdots i_{n}}=v_{i_{1}}v_{i_{2}}\cdots v_{i_{n}}$. The decomposition $\left\{(\lambda_{i},\bm{v}_{i})\right\}_{i\in[n]}$ is called the symmetric orthogonal decomposition (also abbreviated as SOD) of $\bm{\mathcal{T}}$ with individual $\lambda_{i}$ and $\bm{v}_{i}$, respectively, called an eigenvalue and eigenvector of $\bm{\mathcal{T}}$.111For a more detailed discussion on eigenvalues and eigenvectors of SOD tensors, please see [(3)]. The gist of our paper is to find the SOD of $\bm{\mathcal{T}}$ (potentially with perturbations), a recurring problem arising in different contexts including higher-order statistical estimation [(4)], independent component analysis [(5), (6)], and parameter estimation for latent variable models [(7)], just to name a few.

From the expression (1.1), it is quite tempting to find $(\lambda_{i},\bm{v}_{i})$ one by one in a greedy manner using proper deflation procedures. Specifically, one first approximates $\bm{\mathcal{T}}$ by the best rank-one tensor,

[TABLE]

After that, to find the next pair, one modifies the optimization problem (1.2) to exclude the found eigenpair $(\lambda^{\star},\bm{v}^{\star})$. We next review two natural deflation procedures—residual deflation [(8)] and constrained deflation [(9)]—which incorporate the information of $(\lambda^{\star},\bm{v}^{\star})$ into an optimization framework by altering, respectively, the objective and the feasible set of problem (1.2).

Residual deflation

In residual deflation, the rank-one approximation is subtracted from the original tensor, i.e., $\bm{\mathcal{T}}\leftarrow\bm{\mathcal{T}}-\lambda^{\star}\cdot(\bm{v}^{\star})^{\otimes p}$, and then finds the best rank-one approximation to the deflated tensor by solving (1.2) again. The complete scheme, referred to as Successive Rank-One Approximation with Residual Deflation (SROAwRD), is described in Algorithm 1.

Constrained deflation

In constrained deflation, one restricts $\bm{v}$ to be nearly orthogonal to $\pm\bm{v}^{\star}$ by solving problem (1.2) with the additional linear constraints $-\theta\leq\left\langle\bm{v}^{\star},\bm{v}\right\rangle\leq\theta$, where $\theta>0$ is a prescribed parameter. The complete scheme, referred to as Successive Rank-One Approximation with Constrained Deflation (SROAwCD), is described in Algorithm 2. At the $k$-th iteration, rather than deflating the original tensor $\widehat{\bm{\mathcal{T}}}$ by subtracting from it the sum of the $(k-1)$ rank-one tensors $\hat{\lambda}_{1}\hat{\bm{v}}_{1}^{\otimes p}$, $\hat{\lambda}_{2}\hat{\bm{v}}_{2}^{\otimes p}$, $\cdots$, $\hat{\lambda}_{k-1}\hat{\bm{v}}_{k-1}^{\otimes p}$ as the SROAwRD method does, the SROAwCD method imposes the near-orthogonality constraints $|\left\langle\bm{v},\hat{\bm{v}}_{i}\right\rangle|\leq\theta$ for $i=1,2,\ldots,k-1$.

It is not hard to prove that given the SOD tensor $\bm{\mathcal{T}}=\sum_{i\in[n]}\lambda_{i}\bm{v}_{i}^{\otimes p}$ as the input, both SROAwRD and SROAwCD methods are capable of finding the eigenpairs $\left\{(\lambda_{i},\bm{v}_{i})\right\}_{i\in[n]}$ exactly. In this paper, we focus on the more challenging case of tensors that are only close to being SOD. Specifically:

Problem 1**.**

Suppose the SOD tensor $\bm{\mathcal{T}}=\sum_{i\in[n]}\lambda_{i}\bm{v}_{i}^{\otimes p}$, and that the perturbed SOD tensor $\widehat{\bm{\mathcal{T}}}$ is provided as input to the SROAwRD and SROAwCD methods. Characterize the discrepancy between $\{(\lambda_{i},\bm{v}_{i})\}_{i\in[n]}$ and the components $\{(\hat{\lambda}_{i},\hat{\bm{v}}_{i})\}_{i\in[n]}$ found by these methods.

In this paper, we provide positive answers to Problem 1. The characterization for SROAwRD was done in our previous paper [(1)]; we review the results in Section 3. The charaterization for SROAwCD is the main contribution of the present paper. These results can be regarded as higher order generalizations of the Davis-Kahan perturbation result [(10)] for matrix eigen-eigenvalue decomposition, and is not only of mathematical interest but also crucial to applications in signal processing, machine learning and statistics [(4), (5), (6), (7)], where the common interest is to find the underlying eigenpairs $\left\{(\lambda_{i},\bm{v}_{i})\right\}$ but the tensor collected is subject to inevitable perturbations arising from sampling errors, noisy measurements, model specification, numerical errors and so on.

Organization

The rest of the paper is organized as follows. In Section 2, we introduce notation relevant to this paper. In Section 3, we review theoretically what is known about the SROAwRD method. In Section 4, we provide a perturbation analysis for the SROAwCD method.

2 Notation

In this section, we introduce some tensor notation needed in our paper, largely borrowed from [(11)].

Symmetric tensor

A real $p$-way $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, i.e. for any $i_{1},i_{2},\dotsc,i_{p}\in[n]:=\left\{1,2,\ldots n\right\}$,

[TABLE]

for every permutation mapping $\pi$ of $[p]$.

Multilinear map

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 define $\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 two special cases are quite frequently used in the paper:

$\rhd$

$\bm{V}_{i}=\bm{x}\in\mathbb{R}^{n}$ for all $i\in[p]$: $\bm{\mathcal{A}}\bm{x}^{\otimes p}:=\bm{\mathcal{A}}(\bm{x},\bm{x},\dotsc,\bm{x}),$ which defines a homogeneous polynomial of degree $p$.

$\rhd$

$\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 a symmetric tensor $\bm{\mathcal{A}}\in\bigotimes^{p}\mathbb{R}^{n}$, the differentiation result $\nabla_{\bm{x}}\left(\bm{\mathcal{A}}\bm{x}^{\otimes p}\right)=p\cdot\left(\bm{\mathcal{A}}\bm{x}^{\otimes p-1}\right)$ can be established.

Inner product

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

[TABLE]

Tensor norms

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 $\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 also well-known that for symmetric tensors $\bm{\mathcal{A}}$, $\left\|\bm{\mathcal{A}}\right\|$ can be equivalently defined as $\max_{\left\|\bm{x}\right\|=1}|\bm{\mathcal{A}}\bm{x}^{\otimes p}|$ (see, e.g., [(12), (13)]).

3 Review on SROAwRD

Algorithm 1 is intensively studied in the tensor community, though most papers [(14), (8), (15), (16), (17), (18), (12), (13), (19), (20), (21), (22), (7), (23), (24)] focus on the numerical aspects of how to solve the best tensor rank-one approximation (1.2). Regarding theoretical guarantees for the symmetric and orthogonal decomposition, Zhang and Golub [(8)] first prove that SROAwRD outputs the exact symmetric and orthogonal decomposition if the input tensor is symmetric and orthogonally decomposable:

Proposition 1**.**

[(8), Theorem 3.2]** Let $\bm{\mathcal{T}}\in\bigotimes^{p}\mathbb{R}^{n}$ be a symmetric tensor with orthogonal decomposition $\bm{\mathcal{T}}=\sum_{i\in[n]}\lambda_{i}\bm{v}_{i}^{\otimes p}$, where $\lambda_{i}\neq 0$ and $\left\{\bm{v}_{1},\bm{v}_{2},\ldots,\bm{v}_{n}\right\}$ forms an orthonormal basis of $\mathbb{R}^{n}$. Let $\{(\hat{\lambda}_{i},\hat{\bm{v}}_{i})\}_{i\in[n]}$ be the output of Algorithm 1 with input $\bm{\mathcal{T}}$. Then $\bm{\mathcal{T}}=\sum_{i\in[n]}\hat{\lambda}_{i}\hat{\bm{v}}_{i}^{\otimes p}$, and moreover there exists a permutation $\pi$ of $[n]$ such that for each $j\in[n]$,

[TABLE]

The perturbation analysis is recently addressed in [(1)]:

Theorem 2**.**

[(1), Theorem 3.1]** There exists a positive constant $c$ such that the following holds. Let $\widehat{\bm{\mathcal{T}}}:=\bm{\mathcal{T}}+\bm{\mathcal{E}}\in\bigotimes^{p}\mathbb{R}^{n}$, where the ground truth tensor $\bm{\mathcal{T}}$ is symmetric with orthogonal decomposition $\bm{\mathcal{T}}=\sum_{i\in[n]}\lambda_{i}\bm{v}_{i}^{\otimes p}$, $\left\{\bm{v}_{1},\bm{v}_{2},\ldots,\bm{v}_{n}\right\}$ forms an orthonormal basis of $\mathbb{R}^{n}$, $\lambda_{i}\neq 0$ and the perturbation tensor $\bm{\mathcal{E}}$ is symmetric with operator norm $\varepsilon:=\left\|\bm{\mathcal{E}}\right\|$. Assume $\varepsilon\leq c\cdot\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 with input $\widehat{\bm{\mathcal{T}}}$. Then there exists a permutation $\pi$ over $[n]$ such that for each $j\in[n]$:

[TABLE]

Theorem 2 generalizes Proposition 1, and provides perturbation bounds for the SROAwRD method. Specifically, when the operator norm of the perturbation tensor vanishes, i.e. $\varepsilon=0$, Theorem 2 is reduced to Proposition 1; when $\varepsilon$ is small enough (i.e. $\varepsilon=O(1/n^{1/(p-1)})$), the SROAwRD method is able to robustly recover the eigenpairs $\left\{(\lambda_{i},\bm{v}_{i})\right\}_{i\in[n]}$ of the underlying symmetric and orthogonal decomposable tensor $\bm{\mathcal{T}}$.

In Theorem 2, $\varepsilon$ is required to be at most on the order of $1/n^{1/(p-1)}$, which decreases with increasing tensor size. It is interesting to explore whether or not this dimensional dependency is essential:

Open Question 1**.**

Can we provide a better analysis for the SROAwRD method to remove the dimensional dependance on the noise level? Or can we design a concrete example to corroborate the necessity of this dimensional dependency?

The existence of the dimensional dependency, at least for the current proof in Mu et. al. [(1)], can be briefly explained as follows. At the end of the $k$-th iteration, we subtract the rank-one tensor $\hat{\lambda}_{k}\hat{\bm{v}}_{k}^{\otimes p}$ from $\widehat{\bm{\mathcal{T}}}_{k-1}$. Since $(\hat{\lambda}_{k},\hat{\bm{v}}_{k})$ only approximates the underlying truth, this deflation procedure introduces additional errors into $\widehat{\bm{\mathcal{T}}}_{k}$. Although [(1)] has made substantial efforts to reduce the accumulative effect from sequential deflation steps, the perturbation error $\varepsilon$ still needs to depend on the iteration number in order to control the perturbation bounds of the eigenvalue and eigenvector, and we tend to believe that the dimensional dependency in Theorem 2 is necessary.

In contrast, the SROAwCD method, instead of changing the objective, imposes additional constraints, which force the next eigenvector $\hat{\bm{v}}_{k}$ to be nearly orthogonal to $\left\{\hat{\bm{v}}_{1},\hat{\bm{v}}_{2},\ldots,\hat{\bm{v}}_{k-1}\right\}$. As the SROAwCD method alters the search space rather than the objective in the optimization, there is hope that the requirement on the noise level might be dimension-free. In the next section, we will confirm this intuition.

4 SROAwCD

In this section, we establish the first perturbation bounds that have been given for the SROAwCD method for tensor SOD. The main result can be stated as follows:

Theorem 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}\neq 0$ for all $i\in[n]$, and $\bm{\mathcal{E}}$ is a symmetric tensor with operator norm $\varepsilon:=\left\|\bm{\mathcal{E}}\right\|$. Assume $0<\theta\leq{1}/{(2\kappa)}$ and $\varepsilon\leq\theta^{2}\lambda_{\min}/12.5$, where $\kappa:=\lambda_{\max}/\lambda_{\min}$, $\lambda_{\min}:=\min_{i\in[n]}\left|\lambda_{i}\right|$ and $\lambda_{\max}:=\max_{i\in[n]}\left|\lambda_{i}\right|$. Let $\{(\hat{\lambda}_{i},\hat{\bm{v}}_{i})\}_{i\in[n]}$ be the output of Algorithm 2 for input $(\widehat{\bm{\mathcal{T}}},\theta)$. Then there exists a permutation $\pi$ of $[n]$ such that for all $j\in[n]$,

[TABLE]

Theorem 3 guarantees that for an appropriately chosen $\theta$, the SROAwCD method can approximately recover $\left\{(\lambda_{i},\bm{v}_{i})\right\}_{i\in[n]}$ whenever the perturbation error $\varepsilon$ is small. A few remarks immediately come to find. First, Theorem 3 specifies the choice of the parameter $\theta$, which depends on the ratio of the largest to smallest eigenvalues of $\bm{\mathcal{T}}$ in absolute value. In subsection 4.2, we will see this dependency is necessary through numerical studies. Second, in contrast to the SROAwRD method, Theorem 3 does not require the noise level to be dependent on the tensor size. This could be a potential advantage for the SROAwCD method.

The rest of this section is organized as follows. In subsection 4.1, we provide the proof for Theorem 3. In subsection 4.2, we present numerical experiments to corroborate Theorem 3. In subsection 4.3, we discuss issues related to determining the maximum spectral ratio $\kappa$ defined in Theorem 3.

4.1 Proof of Theorem 3

We will prove Theorem 3 by induction. For the base case, we need the perturbation result regarding the best rank-one tensor approximation, which is proven in [(1)] and can be regarded as a generalization of its matrix counterpart [(25), (10)]. In the following, we restate this result [(1), Theorem 2.2] with a minor variation:

Lemma 4**.**

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}\neq 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{\lambda},\hat{\bm{v}})\in\operatorname*{arg\,min}_{\lambda\in\mathbb{R},\left\|\bm{v}\right\|=1}\;\left\|\widehat{\bm{\mathcal{T}}}-\lambda\bm{v}^{\otimes p}\right\|_{F}$. Then there exist $j\in[n]$ such that

[TABLE]

Now we are ready to prove our main Theorem 3.

Proof.

Without loss of generality, we assume $p\geq 3$ is odd, and $\lambda_{i}>0$ for all $i\in[n]$ (as we can always flip the signs of the $\bm{v}_{i}^{\prime}s$ to ensure this). Then problem (1.3) can be equivalently written as

[TABLE]

and $\hat{\lambda}_{k}=\widehat{\bm{\mathcal{T}}}\hat{\bm{v}}_{k}^{\otimes p}$.

To prove the theorem, it suffices to prove that the following property holds for each $k\in[n]$: there is a permutation $\pi$ of $[n]$ such that for every $j\in[k]$,

[TABLE]

We will prove ($*$ ‣ 4.1) by induction.

For the base case $k=1$, Lemma 4 implies that there exists a $j\in[n]$ satisfying

[TABLE]

where we have used the fact that $\varepsilon/\lambda_{j}\leq\varepsilon/\lambda_{\min}\leq\theta^{2}/12.5\leq 1/50.$

Next we assume the induction hypothesis ($*$ ‣ 4.1) is true for $k\in[n-1]$, and prove that there exists an $l\in[n]\backslash\left\{\pi(j):j\in[k]\right\}$ that satisfies

[TABLE]

Denote $\hat{\bm{x}}:=\hat{\bm{v}}_{k+1}$ and $\hat{\lambda}:=\hat{\lambda}_{k+1}$. Then based on (4.3), one has

[TABLE]

and $\hat{\lambda}=\widehat{\bm{\mathcal{T}}}\hat{\bm{x}}^{\otimes p}$. Since $\left\{\bm{v}_{i}\right\}_{i\in[n]}$ forms an orthonormal basis, we may write $\hat{\bm{x}}=\sum_{i\in[n]}x_{i}\bm{v}_{i}$. Without loss of generality, we renumber $\left\{\left(\lambda_{\pi(i)},\bm{v}_{\pi(i)}\right)\right\}_{i\in[k]}$ to $\left\{\left(\lambda_{i},\bm{v}_{i}\right)\right\}_{i\in[k]}$ and renumber $\left\{\left(\lambda_{i},\bm{v}_{i}\right)\right\}_{i\in[n]\backslash\left\{\pi(i)|i\in[k]\right\}}$ to $\left\{\left(\lambda_{i},\bm{v}_{i}\right)\right\}_{i\in[n]\backslash[k]}$, respectively, to satisfy

[TABLE]

In the following, we will show that $l=k+1$ is indeed the index satisfying (4.4). The idea of the rest of the proof is as follows. We first provide lower and upper bounds for $\hat{\lambda}=\widehat{\bm{\mathcal{T}}}\hat{\bm{x}}^{\otimes p}$, based on which, we are able to show that $|\hat{\lambda}-\lambda_{l}|=O(\varepsilon)$ and $1-|\left\langle\hat{\bm{v}},\bm{v}_{l}\right\rangle|=O(\varepsilon)$. However, Theorem 3 requires $1-|\left\langle\hat{\bm{v}},\bm{v}_{l}\right\rangle|=O(\varepsilon^{2})$. To close this gap, we characterize the optimality condition of (4.5), use of which enables us to sharpen the upper bound of $1-|\left\langle\hat{\bm{v}},\bm{v}_{l}\right\rangle|$.

We first consider the lower bound for $\hat{\lambda}$ by finding a $\bm{v}$ that is feasible for (4.5). For each $(i,j)\in[n]\backslash[k]\times[k]$, one has

[TABLE]

where we have used the Cauchy-Schwarz inequality, and the facts $\pi(j)\in[k]$, $\left\langle\bm{v}_{i},\bm{v}_{\pi(j)}\right\rangle=0$, $\varepsilon\leq{\theta^{2}\lambda_{\min}}/{12.5}$ and $\theta\leq{1}/{(2\kappa)}$. Hence, $\left\{\bm{v}_{i}\right\}_{i\in[n]\backslash[k]}$ are all feasible to problem (4.5) and then we can naturally achieve a lower bound for $\hat{\lambda}$, as

[TABLE]

Regarding the upper bound for $\hat{\lambda}$, one has

[TABLE]

where the last line is due to the assumptions made in (4.6), $\left\|\bm{x}\right\|=1$ and $\varepsilon=\left\|\bm{\mathcal{E}}\right\|$.

Combining (4.8) and (4.9), we have

[TABLE]

Also note that

[TABLE]

where we have used the facts that $\theta\leq 1/(2\kappa)$, $\varepsilon\leq\theta^{2}\lambda_{\min}/12.5\leq\lambda_{\min}/50$, and

[TABLE]

Therefore, in order to satisfy (4.10), we must have

[TABLE]

which simplifies (4.10) to

[TABLE]

Based on (4.13), we have that

[TABLE]

Thus, we have achieved the eigenvalue perturbation bound (4.1) promised in the theorem. Next, we will sharpen the eigenvector perturbation bound by exploiting the optimality conditions for problem (4.5).

The key observation is that, at the point $\hat{\bm{x}}$, the constraint $|\left\langle\hat{\bm{v}}_{i},\hat{\bm{x}}\right\rangle|\leq\theta$ is not active. To see this, for any $i\in[k]$,

[TABLE]

where the last line is due to (4.14) and the fact that $\kappa\theta\leq 1/2$ and $\varepsilon\leq\theta^{2}\lambda_{\min}/12.5.$ Therefore, only the equality constraint is active and will be involved in the optimality conditions at the point $\hat{\bm{x}}$. Consider the Lagrangian function at the point $\hat{\bm{x}}$,

[TABLE]

where $\lambda\in\mathbb{R}$ corresponds to the (scaled) Lagrange multiplier for the equality constraint on the norm of $\hat{\bm{x}}$, which we have squared. Since the linear independent constraint qualification [(26), Section 12.3] can be easily verified, by the first-order optimality conditions (a.k.a. KKT condition), there exists a $\bar{\lambda}\in\mathbb{R}$ such that

[TABLE]

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

[TABLE]

Consider the quantity

[TABLE]

Thanks to the intermediate result (4.14), we have

[TABLE]

Moreover, for the term $\left\|\sum_{i\neq k+1}\lambda_{i}x_{i}^{p-1}\bm{v}_{i}\right\|$, we can derive that

[TABLE]

The last line holds due to (4.6) and for $j\in\left\{1,k+2\right\}$,

[TABLE]

where we have used $\sum_{i\in[n]}x_{i}^{2}=1$ and (4.14).

Therefore, by substituting (4.17) and (4.18) into (4.16), one has

[TABLE]

which leads to the desired bound $\left\|\hat{\bm{x}}-\bm{v}_{k+1}\right\|\leq(6.2+4\kappa)\varepsilon/\lambda_{k+1}.$

By mathematical induction, we complete the proof. ∎

4.2 Numerical Experiments

In this subsection, we present three sets of numerical experiments to corroborate our theoretical findings in Theorem 4 regarding the SROAwCD method. We solve the main subproblem (1.3) via the general polynomial solver GloptiPoly 3 (henrion2009gloptipoly, ), which is a global solver based on the Sum-of-Squares (SOS) framework shor87poly ; nesterov00squared ; parrilo2000structured ; lasserre2001global ; parrilo2003semidefinite .

Experiment 1

In this experiment, we will synthetically verify the perturbation bounds stated in Theorem 4. We generate nearly symmetric orthogonally decomposable tensor $\widehat{\bm{\mathcal{T}}}=\bm{\mathcal{T}}+\bm{\mathcal{E}}\in\mathbb{R}^{5\times 5\times 5}$ in the following manner. The underlying symmetric orthogonally decomposable tensor $\bm{\mathcal{T}}$ is set as the diagonal tensor with all diagonal entries equal to 300, i.e., $\bm{\mathcal{T}}=\sum_{i=1}^{5}300\cdot\bm{e}_{i}^{\otimes 3}$, and the perturbation tensors $\bm{\mathcal{E}}$ are produced by symmetrizing a randomly generated $5\times 5\times 5$ tensor whose entries follow standard normal distribution independently. We set $\theta$ to be $1/(2\kappa)=1/2$ (as suggested in Theorem 3). 1000 random instances are tested. Figure 1 plots the histogram of perturbations in both eigenvalue and eigenvector. As depicted in Figure 1, both types of perturbations are well controlled by the bounds provided in Theorem 4.

Experiment 2

In Theorem 4, the parameter $\theta$ is suggested to be set to $1/(2\kappa)$. In this experiment, we compare the performance of SROAwCD with $\theta=1/(2\kappa)=1/2$ and $\theta=0$ based on the criterion

[TABLE]

The tensors are generated in the same way as in the first experiment. Among all the 1000 random cases, the SROAwCD method with $\theta=1/2$ consistently outperforms the one with $\theta=0$. This makes intuitive sense. As $(\hat{\lambda}_{k},\hat{\bm{v}}_{k})$ only approximate the underlying truth, setting $\theta=0$, which forces strict orthogonality, tends to introduce additional errors into the problem.

Experiment 3

In Theorem 4, the parameter $\theta$ is suggested to be set to $1/(2\kappa)$, which depends on $\kappa$. In this experiment, we will demonstrate the necessity of this dependency. We consider the SOD tensor $\bm{\mathcal{T}}=1000\cdot\bm{e}_{1}^{\otimes 3}+\sum_{i=2}^{5}100\cdot\bm{e}_{i}^{\otimes 3}$ with $\bm{\mathcal{E}}=\bm{0}$. We first apply the SROAwCD method with $\theta=1/2$. The output is as follows:

[TABLE]

which deviate greatly from the underlying eigenvalues and eigenvectors of $\bm{\mathcal{T}}$. Next, we apply the SROAwCD method again with $\theta={1}/{(2\kappa)}=1/20$ and the output is as follows:

[TABLE]

which exactly recovers (up to a permutation) the underlying eigenvalues and eigenvectors of $\bm{\mathcal{T}}$.

4.3 Determination of the maximum spectral ratio $\kappa$

As suggested by Theorem 3, to choose a proper $\theta$ for the SROAwCD method, we need a rough estimate for $\kappa$. This is not much of a problem, especially for applications in statistics and machine learning, due to several reasons. First, in many problems, we know the maximum spectral ratio $\kappa$ in advance. For example, if we apply independent component analysis comon1994independent ; comon2010handbook to the dictionary learning model considered in spielman2012exact ; sun2015complete , it is known that $\kappa=1$. Moreover, we can always pick the most favorable estimates for $\kappa$ using cross validation (friedman2001elements, , Section 7.10) based on the prediction errors. Furthermore, as a supplement, we can modify Algorithm 2 to allow the algorithm to determine the appropriate $\theta$ at each step through adaptive learning. We present the analysis of this modification in the appendix.

5 Conclusion

In this paper, we are concerned with finding the (approximate) symmetric and orthogonal decomposition of a nearly symmetric and decomposable (SOD) tensor. Two natural incremental rank-one approximation approaches, the SROAwRD and the SROAwCD methods, have been considered. We first reviewed the existing perturbation bounds for the SROAwRD method. Then we established the first perturbation results for the SROAwCD method that have been given, and discussed issues and potential advantages of this approach. Numerical results were also presented to corroborate our theoretical findings.

We hope our discussion can also shed light on the numerical side. In the SROAwRD method, the main computational bottleneck is the tensor best rank-one approximation problem (1.2), to which a large amount of attention from a numerical optimization point of view has been paid and for which many efficient numerical methods (e.g.  de2000best ; zhang2001rank ; kofidis2002best ; wang2007successive ; kolda2011shifted ; han2012unconstrained ; chen2012maximum ; zhang2012best ; l2014sequential ; jiang2012tensor ; nie2013semidefinite ; yang2014properties ) have been successfully proposed. In the SROAwCD method, the main computational concern is problem (1.3), which is similar to but slightly more complicated than problem (1.2) with additional linear inequalities. Though general-purpose polynomial solvers based on the sum-of-squares framework shor87poly ; nesterov00squared ; parrilo2000structured ; lasserre2001global ; parrilo2003semidefinite ; henrion2009gloptipoly ; sostools can be utilized (as we did in the subsection 4.2), more efficient and scalable methods (e.g. projected gradient method wright1999numerical , semidefinite programming relaxations jiang2012tensor ; nie2013semidefinite ; hu2016note ), specifically tailored to the structure of (1.3), may be anticipated. This is definitely a promising future research direction.

Acknowledgements

We are grateful to the associate editor Tamara G. Kolda and two anonymous reviewers for their helpful suggestions and comments that substantially improved the paper.

Appendix A Adaptive SROAwCD

In this section, we provide a modification to the SROAwCD method that adaptively learns an appropriate $\theta$ at each step based on the information collected so far. The complete algorithm is described in Algorithm 3. Note that we initially set $\theta$ as $1/2$ (which is the largest value allowed in Theorem 3), and gradually reduce it by checking certain conditions.

Theorem 5**.**

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}\neq 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\lambda_{\min}/70\kappa^{2}$, where $\kappa:=\lambda_{\max}/\lambda_{\min}$, $\lambda_{\min}:=\min_{i\in[n]}\left|\lambda_{i}\right|$ and $\lambda_{\max}:=\max_{i\in[n]}\left|\lambda_{i}\right|$. Then Algorithm 3 terminates in a finite number of steps and its output $\{(\theta_{i},\hat{\lambda}_{i},\hat{\bm{v}}_{i})\}_{i\in[n]}$ satisfies:

[TABLE]

and there exists a permutation $\pi$ of $[n]$ such that for each $j\in[n]$

[TABLE]

Remark 1**.**

The condition in line 5 of Algorithm 3 consists of two components, which are desired properties mainly inspired by the proof of Theorem 3. The first desired inequality $\left|\left\langle\hat{\bm{v}}_{k},\hat{\bm{v}}_{i}\right\rangle\right|<\hat{\lambda}_{k}/1.35\hat{\lambda}_{i}$ would allow us to establish properties similar to (4.12). The second desired inequality $\left|\left\langle\hat{\bm{v}}_{k},\hat{\bm{v}}_{i}\right\rangle\right|<\theta$ would help us make use of the optimality condition to sharpen the perturbation bounds for the eigenvectors. The constants chosen in Algorithm 3 and Theorem 5, mainly for illustrative purposes, might be better optimized.

Proof.

Similar to the proof for Theorem 3, without loss of generality, we assume $p\geq 3$ is odd, and $\lambda_{i}>0$ for all $i\in[n]$. Then problem (A.1) can be equivalently written as

[TABLE]

and $\hat{\lambda}_{k}=\widehat{\bm{\mathcal{T}}}\hat{\bm{v}}_{k}^{\otimes p}$.

Our proof is by induction.

The base case regarding $(\theta_{1},\hat{\lambda}_{1},\hat{\bm{v}}_{1})$ is the same as the base case in the proof of Theorem 3.

We now make the induction hypothesis that for some $k\in[n-1]$, $\{(\theta_{i},\hat{\lambda}_{i},\hat{\bm{v}}_{i})\}_{i\in[k]}$ satisfies

[TABLE]

and there exists a permutation $\pi$ of $[n]$ such that

[TABLE]

Then we are left to prove that

[TABLE]

and there exists an $l\in[n]\backslash\left\{\pi(j):j\in[k]\right\}$ that satisfies

[TABLE]

To prove that $\theta_{k+1}>\frac{0.96}{2\kappa}$, we show that whenever $\theta\in(\frac{0.96}{2\kappa},\frac{1}{2\kappa}]$, the condition for the while loop in line 5 of Algorithm 3 will always be satisfied, so $\theta$ can never be reduced to any value below $\frac{0.96}{2\kappa}$.

Consider $\theta\in(\frac{0.96}{2\kappa},\frac{1}{2\kappa}]$, and denote

[TABLE]

and $\hat{\lambda}=\widehat{\bm{\mathcal{T}}}\hat{\bm{x}}^{\otimes p}$. As $\theta\leq 1/2\kappa$ and $\varepsilon\leq\lambda_{\min}/70\kappa^{2}\leq\frac{0.96^{2}\lambda_{\min}}{4\cdot 12.5\cdot\kappa^{2}}\leq\theta^{2}\lambda_{\min}/12.5$, results from Theorem 3 and its proof can be directly borrowed. Based on (4.7), we know that for any $i\in\left\{\pi(j)\;|\;j\in[n]\backslash[k]\right\}$, $\bm{v}_{i}$ is feasible to problem (A.5). Then it can be easily verified that

[TABLE]

Moreover, for any $i\in[k]$,

[TABLE]

Using (A.6) and (A.7), we obtain that for each $i\in[k]$,

[TABLE]

Based on (4.15), we also know that for all $i\in[k]$,

[TABLE]

So, the combination of (A.8) and (A.9) leads to

[TABLE]

which implies that $\hat{\bm{x}}$ satisfies the condition in the while loop. Therefore, as we argued previously, we must have $\theta_{k+1}>\frac{0.96}{2\kappa}$. So $\theta_{k+1}$ is either in $(\frac{0.96}{2\kappa},\frac{1}{2\kappa}]$ or in $(\frac{1}{2\kappa},\frac{1}{2}]$.

For the first case, i.e. $\theta_{k+1}\in(\frac{0.96}{2\kappa},\frac{1}{2\kappa}]$, we can directly establish the result by using the argument for the induction hypothesis ($*$ ‣ 4.1) in the proof of Theorem 3.

Hence in the following, we only focus on the second case where $\theta_{k+1}\in(\frac{1}{2\kappa},\frac{1}{2}]$.

Denote $\hat{\bm{x}}=\sum_{i\in[n]}x_{i}\bm{v}_{i}:=\hat{\bm{v}}_{k+1}$ and $\hat{\lambda}:=\hat{\lambda}_{k+1}.$ Without loss of generality, we renumber $\left\{\left(\lambda_{\pi(i)},\bm{v}_{\pi(i}\right)\right\}_{i\in[k]}$ to $\left\{\left(\lambda_{i},\bm{v}_{i}\right)\right\}_{i\in[k]}$ and renumber $\left\{\left(\lambda_{i},\bm{v}_{i}\right)\right\}_{i\in[n]\backslash\left\{\pi(i)|i\in[k]\right\}}$ to $\left\{\left(\lambda_{i},\bm{v}_{i}\right)\right\}_{i\in[n]\backslash[k]}$, respectively, to satisfy

[TABLE]

In the following, we will show that $l=k+1$ is the index satisfying (A.4).

Based on (A.2),

[TABLE]

and $\hat{\lambda}=\widehat{\bm{\mathcal{T}}}\hat{\bm{x}}^{\otimes p}.$ We now bound $\hat{\lambda}$ from below and above.

We first consider the lower bound by finding a $\bm{v}$ that is feasible for (A.12). For any $(i,j)\in[n]\backslash[k]\times[k]$, one has

[TABLE]

Hence, $\left\{\bm{v}_{i}\right\}_{i\in[n]\backslash[k]}$ are all feasible to problem (A.12) and then we can easily achieve a lower bound for $\hat{\lambda}$, as

[TABLE]

Regarding the upper bound, one has

[TABLE]

as in (4.9).

Combining (A.13) and (A.14), we have

[TABLE]

Also note that

[TABLE]

Here, we have used $\varepsilon\leq\frac{\hat{\lambda}}{69}$, due to $\varepsilon\leq\frac{\lambda_{\min}}{70}$ and $\hat{\lambda}\geq\lambda_{\min}-\varepsilon.$ Therefore, in order to satisfy (A.15), we must have

[TABLE]

which simplifies (A.15) to

[TABLE]

Hence,

[TABLE]

The eigenvector perturbation bound can be sharpened by exploiting the optimality condition of problem (A.12) as what we did in the proof of Theorem 3. As explicitly required in the algorithm, the constraint $|\left\langle\hat{\bm{v}}_{i},\hat{\bm{x}}\right\rangle|\leq\theta_{k}$ is not active for any $i\in[k]$ at the point $\hat{\bm{x}}$. Then by the optimality condition, we again have

[TABLE]

By applying exactly the same argument as in the proof of Theorem 3, we can obtain

[TABLE]

which leads to $\left\|\hat{\bm{x}}-\bm{v}_{k+1}\right\|\leq(6.2+4\kappa)\varepsilon/\lambda_{k+1}.$

By mathematical induction, we have completed the proof. ∎