Low-Rank Principal Eigenmatrix Analysis

TL;DR

This paper introduces low-rank principal eigenmatrix analysis, a novel approach that leverages low-rank structures in eigenvectors for high-dimensional data, offering an alternative to sparse PCA with efficient algorithms.

Contribution

It proposes a new low-rank eigenmatrix analysis method, including a rank-truncated power algorithm, with theoretical guarantees and practical efficiency.

Findings

The method performs competitively on synthetic datasets.

The proposed algorithm is computationally efficient.

It effectively captures low-rank structures in eigenvectors.

Abstract

Sparse PCA is a widely used technique for high-dimensional data analysis. In this paper, we propose a new method called low-rank principal eigenmatrix analysis. Different from sparse PCA, the dominant eigenvectors are allowed to be dense but are assumed to have a low-rank structure when matricized appropriately. Such a structure arises naturally in several practical cases: Indeed the top eigenvector of a circulant matrix, when matricized appropriately is a rank-1 matrix. We propose a matricized rank-truncated power method that could be efficiently implemented and establish its computational and statistical properties. Extensive experiments on several synthetic data sets demonstrate the competitive empirical performance of our method.

Equations

$${\mbox{\boldmath$A$}}={}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmath$A$}}}+{\mbox{\boldmath$E$}},$$

$$\displaystyle\underset{{\mbox{\boldmath$x$}}\in\mathbb{R}^{d}}{\text{maximize}}$$

$$\displaystyle\|{\mbox{\boldmath$x$}}\|_{\ell_{2}}=1,\;\|{\mbox{\boldmath$x$}}\|_{\ell_{0}}\leq s.$$

$$\displaystyle\underset{{\mbox{\boldmath$X$}}\in\mathbb{R}^{p\times p}}{\text{maximize}}$$

$$\displaystyle\|{\mbox{\boldmath$X$}}\|_{F}=1,$$

$$\displaystyle{\rm rank}\left({\mbox{\boldmath$X$}}\right)\leq k\ll p.$$

$$RankTrunc({\mbox{\boldmath$X$}},k)=P_{U}{\mbox{\boldmath$X$}}P_{V}$$

$$RankTrunc({\mbox{\boldmath$x$}},k)={\textsc{vec}}\left(RankTrunc({\textsc{mat}}\left({\mbox{\boldmath$x$}}\right),k)\right)$$

$$\displaystyle\underset{{\mbox{\boldmath$x$}}}{\text{maximize}}$$

$$\displaystyle\|{\mbox{\boldmath$x$}}\|_{\ell_{2}}=1,$$

$$\displaystyle{\rm rank}\left({\mbox{\boldmath$X$}}\right)\leq k<p.$$

$$\displaystyle({\mbox{\boldmath$x$}}_{t}-{\mbox{\boldmath$x$}}_{t-1})^{T}\,{\mbox{\boldmath$A$}}\,({\mbox{\boldmath$x$}}_{t}-{\mbox{\boldmath$x$}}_{t-1})$$

$$\displaystyle=({\mbox{\boldmath$x$}}_{t}-{\mbox{\boldmath$x$}}_{t-1})^{T}\,P_{V\otimes U}{\mbox{\boldmath$A$}}P_{V\otimes U}\,({\mbox{\boldmath$x$}}_{t}-{\mbox{\boldmath$x$}}_{t-1})$$

$$\ge 0.$$

$$\displaystyle Q({\textsc{vec}}\left({\mbox{\boldmath$X$}}_{t}\right))-Q({\textsc{vec}}\left({\mbox{\boldmath$X$}}_{t-1}\right))$$

$$\displaystyle=L({\mbox{\boldmath$x$}}_{t},{\mbox{\boldmath$x$}}_{t-1})+({\mbox{\boldmath$x$}}_{t}-{\mbox{\boldmath$x$}}_{t-1})^{T}\,{\mbox{\boldmath$A$}}\,({\mbox{\boldmath$x$}}_{t}-{\mbox{\boldmath$x$}}_{t-1})$$

$$\displaystyle\geq L({\mbox{\boldmath$x$}}_{t},{\mbox{\boldmath$x$}}_{t-1})\geq L({\mbox{\boldmath$x$}}_{t-1},{\mbox{\boldmath$x$}}_{t-1})$$

$$=0.$$

$$\displaystyle\Big{\lVert}P_{U}{\mbox{\boldmath$W$}}_{i}P_{V}^{\top}\Big{\rVert}_{F}\leq\sqrt{{\rm rank}\left(P_{U}{\mbox{\boldmath$W$}}_{i}P_{V}^{\top}\right)}\,\Big{\lVert}P_{U}{\mbox{\boldmath$W$}}_{i}P_{V}^{\top}\Big{\rVert}$$

$$\displaystyle\qquad\leq\sqrt{k}\Big{\lVert}{\mbox{\boldmath$W$}}_{i}\Big{\rVert}=O(\sqrt{k/p}).$$

$$\displaystyle\rho^{2}({\mbox{\boldmath$E$}}_{UV})=\lVert{\mbox{\boldmath$E$}}_{UV}\rVert^{2}=\lVert{\mbox{\boldmath$Q$}}{\mbox{\boldmath$\Lambda$}}{\mbox{\boldmath$Q$}}^{\top}{\mbox{\boldmath$Q$}}{\mbox{\boldmath$\Lambda$}}{\mbox{\boldmath$Q$}}^{\top}\rVert$$

$$\displaystyle=\lVert{\mbox{\boldmath$Q$}}{\mbox{\boldmath$\Lambda$}}{\mbox{\boldmath$Q$}}^{\top}{\mbox{\boldmath$Q$}}{\mbox{\boldmath$\Lambda$}}{\mbox{\boldmath$Q$}}^{\top}\rVert$$

$$\displaystyle\leq\lVert{\mbox{\boldmath$\Lambda$}}{\mbox{\boldmath$Q$}}^{\top}{\mbox{\boldmath$Q$}}{\mbox{\boldmath$\Lambda$}}\rVert\lVert{\mbox{\boldmath$Q$}}\rVert^{2}$$

$$\displaystyle=\lVert P_{V\otimes U}{\mbox{\boldmath$E$}}{\mbox{\boldmath$W$}}\rVert^{2}\lVert{\mbox{\boldmath$Q$}}\rVert^{2}$$

$$\displaystyle=O(k/p)\cdot\rho^{2}({\mbox{\boldmath$E$}}).$$

$$\displaystyle\rho\left({\mbox{\boldmath$E$}}_{UV}\right)=O(\sqrt{k/p})\cdot\rho\left({\mbox{\boldmath$E$}}\right).$$

$$\displaystyle\gamma\coloneqq\frac{\lambda-\triangle\lambda+\rho\left({\mbox{\boldmath$E$}}_{UV}\right)}{\lambda-\rho\left({\mbox{\boldmath$E$}}_{UV}\right)},$$

$$\delta :=\frac{\sqrt{2}}{\sqrt{1+{\left(\kappa -2\right)}^2}}.$$

$$\mu :=\sqrt{1+2((\bar{k}/k{)}^{1/2}+\bar{k}/k)}\sqrt{1-0.5\theta (1+\theta )(1-{\gamma }^2)}<1.$$

$$\sqrt{1-|\langle{\mbox{\boldmath$X$}}_{0},{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmath$X$}}}\rangle|}\leq\sqrt{10}\delta/(1-\mu)$$

$$\sqrt{1-|\langle{\mbox{\boldmath$X$}}_{t},{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmath$X$}}}\rangle|}\leq\mu^{t}\sqrt{1-|\langle{\mbox{\boldmath$X$}}_{0},{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmath$X$}}}\rangle|}+\sqrt{10}\delta/(1-\mu).$$

$$|\langle{\mbox{\boldmath$X$}}_{0},{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmath$X$}}}\rangle|\geq C\left(\rho({\mbox{\boldmath$E$}}_{UV})+\sqrt{\bar{k}/k}\right)$$

$$\|{\mbox{\boldmath$X$}}_{t}-{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmath$X$}}}\|_{F}=O\left(\rho({\mbox{\boldmath$E$}}_{UV})\right)=O\left(\sqrt{k/p}\,\rho({\mbox{\boldmath$E$}})\right).$$

$${}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmath$A$}}}=\lambda_{1}\bar{\mbox{\boldmath$x$}}\bar{\mbox{\boldmath$x$}}^{\top}+{\mbox{\boldmath$\Sigma$}}_{\epsilon},$$

$${\mbox{\boldmath$A$}}=\frac{1}{n}\sum_{i=1}^{n}{\mbox{\boldmath$y$}}_{i}{\mbox{\boldmath$y$}}_{i}^{\top}.$$

$$\displaystyle\gamma\coloneqq\frac{\lambda-\triangle\lambda+\rho\left({\mbox{\boldmath$E$}}_{UV}\right)}{\lambda-\rho\left({\mbox{\boldmath$E$}}_{UV}\right)}.$$

$$\displaystyle\lVert{\mbox{\boldmath$x$}}_{1}({\mbox{\boldmath$A$}}_{UV})-\bar{\mbox{\boldmath$x$}}\rVert_{\ell_{2}}\leq\delta\coloneqq\frac{\sqrt{2}}{\sqrt{1+\left(\kappa-2\right)^{2}}}.$$

$$\displaystyle\lambda_{1}({\mbox{\boldmath$A$}}_{UV})\geq\lambda_{1}({}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmath$A$}}}_{UV})+\lambda_{p}({\mbox{\boldmath$E$}}_{UV})\geq\lambda-\rho({\mbox{\boldmath$E$}}_{UV})$$

Full text

1 Introduction

Principal Component Analysis (PCA) is a popular technique for analyzing big data with a rich theoretical literature, along with several real-world applications in analyzing financial data, matrix completion, network analysis, gene expression analysis, to name a few. We refer the reader to Izenman (2008); Jolliffe (2011); Fan et al. (2018) for a recent survey. Given a $d\times d$ symmetric positive semidefinite matrix $A$, PCA finds the unit vector {\mbox{\boldmathx}}\in\mathbb{R}^{d}, also called as the Principal Component (PC), that maximizes the quadratic form {\mbox{\boldmathx}}^{\top}{\mbox{\boldmathA}}{\mbox{\boldmathx}}. In statistical machine learning, $A$ is often taken to be the sample covariance matrix of a $d$-dimensional random vector, estimated with $n$ sample points. Several theoretical results concerning the consistency and rates of convergence of the estimated principal components are well-known in the setting when $d$ is fixed and $n\to\infty$. See for example Anderson (1963); Muirhead (2009).

On the other hand, contemporary datasets often have the number of input variables ($d$) comparable with or even much larger than the number of samples ($n$). Asymptotic properties of the PCs under the setting that $n\to\infty$ such that $d/n\to c\in(0,\infty)$ have been studied by several authors; See for example Nadler (2008); Jung and Marron (2009); Johnstone and Lu (2012); Johnstone and Paul (2018). At the risk of sounding non-rigorous, the main conclusion of that line of work is that the estimated PC is not a consistent estimator of the true PC in this setting. Such a concern lead to the development of sparse PCA Zou et al. (2006), where it is assumed that only $s$ of the $d$ components in the unit vector $x$ are non-zero. Such a technique, improves the interpretation of the PC and alleviates some of the problems associated with the consistency of PCA in high-dimensional low-sample size scenario Birnbaum et al. (2013); Vu and Lei (2012, 2013); Cai et al. (2013). But the question, when are the PCs of a covariance matrix truly sparse?, is not well understood. So far, the main motivation for sparse PCA has been the fact that high-dimensional estimation is possible only under the structural sparsity assumptions, which is rather unconvincing. To the best of our knowledge, Lei and Vu (2015) is the only work that analyzes sparse PCA in an agnostic setting, i.e., when the true PC is not assumed to be sparse.

In this paper we propose a novel method called Low-Rank Principal Eigenmatrix Analysis, where the PC is assumed to be low-rank when matricized appropriately. A main motivation for proposing such a structure is from the case of circulant and Toeplitz covariance matrices that arises in the context of stochastic processes used to model financial applications and medical imaging dataset Christensen (2007); Snyder et al. (1989); Roberts and Ephraim (2000); Cai et al. (2013). Consider the case when $d=p^{2}$. It is known that circulant matrices have FFT matrices as eigenvectors Gray (2006); Davis (2012).

Example 1** (Circulant Matrix).**

Theorem 3.1 in Gray (2006) shows that every circulant matrix {\mbox{\boldmathC}}\in\mathbb{R}^{p^{2}\times p^{2}} has eigenvectors of the form {\mbox{\boldmathy}}=(1,\rho,\cdots,\rho^{p^{2}-1})^{\top}\in\mathbb{R}^{p^{2}}. If we reshape the vector $y$ according to Definition 1.2, then the matrix {\textsc{mat}}\left({\mbox{\boldmathy}}\right)\in\mathbb{R}^{p\times p} is an exactly rank one matrix, i.e. {\textsc{mat}}\left({\mbox{\boldmathy}}\right)={\mbox{\boldmathu}}{\mbox{\boldmathv}}^{\top} where {\mbox{\boldmathu}}=\begin{pmatrix}1&\rho&\cdots&\rho^{p}\end{pmatrix}^{\top} and {\mbox{\boldmathv}}=\begin{pmatrix}1&\rho^{p+1}&\cdots&\rho^{p(p-1)+1}\end{pmatrix}^{\top}.

Example 2** (Toeplitz, Diagonally Dominant and Kronecker Structred Covariance Matrices).**

In Figure 1, we plot the spectrum of the matricized version of top eigenvector of a Toeplitz, diagonally dominant, kronecker-product structured covariance matrices Werner et al. (2008), along that of a unstructured covariance matrix. The dimensions of the covariance matrix is $10000\times 10000$. The top eigenvector is hence of dimension $10000$. The reshaping according to Definition 1.2, leads to a $100\times 100$ matrix. We see that, for the case of Toeplitz matrices and Diagonally Dominant matrices, the top eigenmatrix could almost always be well-aproximated by an extremely low-rank matrix. Furthermore, for covariance matrices with Kronecker product structure, the rank of the top eigenmatrix is well approximated by a low-rank approximation: for a $10000$ dimensional principal component represented as a $100\times 100$ matrix a rank-40 approximation captures the bulk. Note that this obviously does not hold in the case of general covariance matrices.

A noteworthy observation in the above examples is that while the reshaped eigenvectors low-rank, they are never sparse. Thus in these cases, while the sparsity assumption on the eigenvector might appear rather superficial, the low-rank assumption appears naturally. This observation serves as a motivation for the low-rank model proposed in this work.

Organization: The remaining of this paper is organized as follows: §2 describes the low-rank principal eigenmatrix problem and proposes a iterative algorithm for computing the estimator. §2.2 analyzes the convergence rates of the proposed algorithm. §4 provides numerical results that confirm the relevance of our theoretical predictions. We end this section with the list of notations we follow in the rest of the paper.

Notation: We use bold-faced letter $y$ to denote vectors, and bold-faced capital letters such as $A$ to denote matrices. For a vector $y$, \lVert{\mbox{\boldmathy}}\rVert_{\ell_{0}} and \lVert{\mbox{\boldmathy}}\rVert_{\ell_{2}} represents the $L_{0}$ and $L_{2}$ norms respectively. For any matrix {\mbox{\boldmathX}}\in\mathbb{R}^{d_{1}\times d_{2}}, we denote its singular values by s_{1}({\mbox{\boldmathX}})\geq\cdots\geq s_{d}({\mbox{\boldmathX}})\geq 0 where $d=\min(d_{1},d_{2})$. The inner product in {\mbox{\boldmathX}}\in\mathbb{R}^{d_{1}\times d_{2}} can be written in matrix form as \langle{\mbox{\boldmathX}},{\mbox{\boldmathY}}\rangle=\mathbf{Tr}({\mbox{\boldmathX}}^{\top}{\mbox{\boldmathY}}). Let \lVert{\mbox{\boldmathX}}\rVert_{F}, \lVert{\mbox{\boldmathX}}\rVert and \lVert{\mbox{\boldmathX}}\rVert_{\star} denote the Frobenius, operator and nuclear norm of $X$, respectively. Let $\mathbb{S}^{d\times d}$ denote the set of symmetric matrices. For any matrix {\mbox{\boldmathA}}\in\mathbb{S}^{d\times d}, we denote its eigenvalues by \lambda_{\max}({\mbox{\boldmathA}})=\lambda_{1}({\mbox{\boldmathA}})\geq\cdots\geq\lambda_{d}({\mbox{\boldmathA}})=\lambda_{\min}({\mbox{\boldmathA}}). Furthermore, we call the eigenvector associated with the top eigenvalue $\lambda_{\max}(A)$ as the top eigenvector. We use \rho({\mbox{\boldmathA}}) to denote the spectral norm of $A$, i.e. \max\{|\lambda_{\max}({\mbox{\boldmathA}})|,|\lambda_{\min}({\mbox{\boldmathA}})|\}. In the rest of the paper, we define Q({\mbox{\boldmathx}})\coloneqq{\mbox{\boldmathx}}^{\top}{\mbox{\boldmathA}}{\mbox{\boldmathx}}, $P_{U}$ as the projection matrix onto the column space of matrix $U$, and {\mbox{\boldmathA}}_{UV}=P_{V\otimes U}{\mbox{\boldmathA}}P_{V\otimes U}. We next define vectorization of a matrix and matricization of a vector in Definition 1.1 and 1.2 respectively.

Definition 1.1** (Vectorization).**

For a matrix {\mbox{\boldmathX}}\in\mathbb{R}^{p\times p}, we define the vectorization of the matrix {\textsc{vec}}\left({\mbox{\boldmathX}}\right)=\mathbb{R}^{p^{2}} to be a vector constructed by stacking columns of $X$ together.

Definition 1.2** (Matricization).**

Given a vector {\mbox{\boldmathx}}\in\mathbb{R}^{p^{2}}, we define the matricization of a vector by {\textsc{mat}}\left({\mbox{\boldmathx}}\right):\mathbb{R}^{p^{2}}\mapsto\mathbb{R}^{p\times p} where the first $p$ coordinates of the vectors form the first column of the matrix and second $p$ indices form the second column and so on. Hence {\mbox{\boldmathX}}_{i,j}={\mbox{\boldmathx}}_{(j-1)p+i}.

Note that the above operations could also be defined based on rows of the matrix $X$. In the rest of the paper, we only refer to Definition 1.1 and 1.2 for the above operations. Furthermore, in the rest of the paper, we fix $d=p^{2}$, for the sake of simplicity. We note that there is nothing special about $A$ being square and all of our discussion would apply to arbitrary rectangular matrices as well. The advantage of focusing on square matrices is a simplified exposition and reduction in the number of parameters of which we need to keep track. In the case of $d$ not being a perfect square, Definition 1.2 could be changed to handle rectangular matrices.

2 Low-Rank Principal Eigenmatrix Analysis

In this section, we first introduce the Low-rank Principal Eigenmatrix Analysis problem in §2.1 and then introduce the matricized rank-truncated power method for solving the problem in §2.2.

2.1 The Model and The Problem

Consider the following noisy ($d\times d$) matrix model:

[TABLE]

where $A$ is a noisy observation of the true signal matrix {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}} with the noise matrix represented by the matrix $E$. As a motivating application, we consider $A$ to be the empirical covariance matrix and {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}} to the true covariance matrix. But the above model is general and applies to several other situations, for example, $A$ could be the noisy version of the true image {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}}, represented in term of pixel arrays. The problem of PCA considers estimating the top eigenvector, \bar{\mbox{\boldmathx}}, of the matrix {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}} based on the observed matrix $A$. Sparse PCA assumes \bar{\mbox{\boldmathx}} is $s$-sparse and estimates it by following constrained quadratic maximization problem:

[TABLE]

Low-rank Principal Eigenmatrix Analysis assumes that the top eigenvector \bar{\mbox{\boldmathx}} of $A$, when matricized as {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}, according to the Definition 1.2, has low-rank. Note that, we are still estimating the top eigenvector. But the low-rank structure is made on the matricized eigenvector {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}. Hence we call it the eigenmatrix, i.e., {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}={\textsc{mat}}\left(\bar{\mbox{\boldmathx}}\right) is the top eigenmatrix of $A$. Furthermore, though {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}} is a low-rank matrix, the eigenvector \bar{\mbox{\boldmathx}}={\textsc{vec}}\left({}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}\right) can be dense, different from sparse PCA. Before proceeding, we make the following remarks. First, as discussed in §1, our assumption $d=p^{2}$ is done only for the sake of convenience. In the case of general $d$, one could assume that there exist a $p_{1}\times p_{2}$ rectangular matricization for which such a low-rank structure exists. Next, note that one could assume a more general model, where the vector \bar{\mbox{\boldmathx}}, after a permutation of its indices followed by the matricization operation has a low-rank. Further, the permutation could be estimated as well. In this work, we do not concentrate on the general model, but we plan to address this in the future. Our goal in this work is to provide a deeper understanding of the simpler model.

Based on our assumption, the Low-Rank Principal Eigenmatrix Estimator is given by the following optimization problem (2):

[TABLE]

Note that the above optimization problem is the natural maximum likelihood solution for estimating the top eigenmatrix {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}} along with the constraints motivated by our assumption. The above problem is a non-convex problem – indeed the presence of the rank constraint along with the normalization constraint makes it highly non-convex. In this next section, we propose an iterative method for solving the above optimization problem.

2.2 smart-pm for Solving Problem 2

In general, problem (2) is non-convex and NP-hard. In order to solve it efficiently, we propose a Sequentially MATricized Rank Truncated-Power Method (smart-pm) outlined in Algorithm 1. It is an iterative procedure based on the standard power method for eigenvalue problems, while maintaining the desired matricized low-rank structure for the intermediate solutions.

Definition 2.1**.**

We define Matricized Rank-Truncation operator as following,

[TABLE]

for a matrix $X$ and

[TABLE]

for a vector $x$, where $U$ and $V$ consists of the first $k$ columns of the left and right singular matrices of $X$, respectively. $P_{U}$ and $P_{V}$ are rank-$k$ projection matrix onto $U$ and $V$, respectively. For a matrix, operation RankTrunc({\mbox{\boldmathX}},k) is basically the truncated SVD that retains top $k$ singular values and sets the rest to zero.

Given an initial approximation {\mbox{\boldmathX}}_{0} and rank $k$, Algorithm 1 generates a sequence of intermediate low-rank eigenmatrices {\mbox{\boldmathX}}_{1},{\mbox{\boldmathX}}_{2},\ldots satisfying \lVert{\mbox{\boldmathX}}_{t}\rVert_{F}=1 and {\rm rank}\left({\mbox{\boldmathX}}_{t}\right)\leq k. At each iteration, the computational complexity is in $O(p^{4}+p^{2}k)$ which is $O(d^{2})=O(p^{4})$ for matrix-vector product {\mbox{\boldmathA}}{\textsc{vec}}\left({\mbox{\boldmathX}}_{t-1}\right) and $O(p^{2}k)$ for getting largest $k$ singular values and corresponding singular vectors. The extra computation time $O(p^{2}k)$ is negligible since $k$ is often small and $O(p^{4})$ is often the dominant term. In the next section we analyze the smart-pm algorithm and discuss its rates of convergence.

3 Theoretical Analysis of smart-pm

In this section, we provide a computational and statistical convergence analysis of the smart-pm algorithm. Our first result is on the algorithmic convergence of smart-pm. Specifically, we show that when all the rank $4k^{2}$ projections of the matrix $A$ are positive semi-definite, the iterates of the smart-pm algorithm are monotonically increasing in terms of the objective value it is optimizing.

Proposition 1**.**

If all the rank $4k^{2}$ projections of the matrix $A$ are positive semi-definite, then the sequence \{Q({\textsc{vec}}\left({\mbox{\boldmathX}}_{t}\right))\}_{t\geq 1} is a monotonically increasing sequence, where {\mbox{\boldmathX}}_{t} is obtained from the SMART-PM algorithm.

Proof.

Note that the iterate {\mbox{\boldmathx}}_{t}={\textsc{vec}}\left({\mbox{\boldmathX}}_{t}\right) in the smart-pm algorithm solves the following constrained optimization problem:

[TABLE]

For any \textsc{rank}({\mbox{\boldmathX}})\leq k and \textsc{rank}({\mbox{\boldmathX}}_{t-1})\leq k, we have \textsc{rank}({\mbox{\boldmathX}}-{\mbox{\boldmathX}}_{t-1})\leq 2k. Suppose {\mbox{\boldmathX}}-{\mbox{\boldmathX}}_{t-1} assumes the following form of SVD decomposition {\mbox{\boldmathX}}-{\mbox{\boldmathX}}_{t-1}={\mbox{\boldmathU}}{\mbox{\boldmathD}}_{2k}{\mbox{\boldmathV}}^{\top} where $U$ and $V$ are left and right singular matrix, respectively, then {\mbox{\boldmathX}}-{\mbox{\boldmathX}}_{t-1}=P_{U}({\mbox{\boldmathX}}-{\mbox{\boldmathX}}_{t-1})P_{V} and {\mbox{\boldmathx}}-{\mbox{\boldmathx}}_{t-1}=P_{V\otimes U}({\mbox{\boldmathx}}-{\mbox{\boldmathx}}_{t-1}), where $P_{U}$ is the projection matrix onto the column space of matrix $U$. Since all the rank $4k^{2}$ projections of the matrix $A$ are positive semi-definite, we have

[TABLE]

Clearly,

[TABLE]

The second inequality holds because {\mbox{\boldmathx}}_{t} by definition maximizes L({\mbox{\boldmathx}},{\mbox{\boldmathx}}_{t-1}) at step $t$. ∎

The above results is a purely computational result justifying the use of smart-pm method as a heuristic for solving the non-convex problem in Equation 2. We next consider the general noisy matrix perturbation model (1) and show that if matrix {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}} has a unique low-rank (or approximately low-rank) top eigenmatrix, then under suitable condition, smart-pm can estimate this eigenmatrix from the noisy observation $A$, under certain conditions on the initial matrix {\mbox{\boldmathX}}_{0}. In order to establish such a result, we first list the set of assumptions required, precisely.

Assumption 1**.**

Assume that the largest eigenvalue of {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}}\in\mathbb{S}^{p^{2}} is \lambda=\lambda_{\max}({}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}})>0 that is non-degenerate, with a gap \triangle\lambda=\lambda-\max_{j>1}|\lambda_{j}({}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}})| between the largest and the remaining eigenvalues. Moreover, assume that the eigenvector \bar{\mbox{\boldmathx}}\in\mathbb{R}^{d} corresponding to the dominant eigenvalue $\lambda$ is of rank $k$ when matricized, i.e. {\textsc{mat}}\left(\bar{\mbox{\boldmathx}}\right)\in\mathbb{R}^{p\times p} is of rank $k\leq p$.

Assumption 2**.**

Assume that {\mbox{\boldmathE}}\in\mathbb{S}^{p^{2}\times p^{2}} does not have low-rank eigenmatrices. Mathematically, let $w$ be any eigenvector of $E$, \|{\mbox{\boldmathw}}\|_{\ell_{2}}^{2}=1, and {\mbox{\boldmathW}}={\textsc{mat}}\left({\mbox{\boldmathw}}\right)\in\mathbb{R}^{p\times p} be the corresponding eigenmatrix. $W$ is of rank $p$ and \max\{|\sigma_{1}({\mbox{\boldmathW}})|,|\sigma_{p}({\mbox{\boldmathW}})|\}\asymp\sqrt{1/p}.

Remark 1**.**

Assumption 1 posits that the top eigenvector (or the eigenmatrix) of {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}} is isolated. Such an analysis is common in the analysis of power method. Assumption 2 is crucial and captures a particular conditions on the error matrix $E$ under which the smart-pm has efficient estimation rates. It posits that the eigenmatrices of $E$ have are high-rank with a flat spectrum. This assumption is complementary to the structural assumption made on the true signal to be estimated, i.e., top eigenmatrix {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}} is low-rank. For example, if the true signal matrix {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}} is corrupted by an additive standard Gaussian iid noise matrix $E$, then with high probability such an assumption holds naturally.

Now we first show that, under Assumption 2 when the noise matrix $E$ does not have low-rank eigenmatrices, the spectral norm of the low rank projection \rho(P_{V\otimes U}{\mbox{\boldmathE}}P_{V\otimes U}) can be small, even when \rho({\mbox{\boldmathE}}) is large.

Lemma 1**.**

Consider low-rank projection $P_{U}$ and $P_{V}$ where ${\rm rank}\left(P_{U}\right)={\rm rank}\left(P_{V}\right)=k$ and $1\leq k\leq p$. For any matrix {\mbox{\boldmathE}}\in\mathbb{S}^{p^{2}\times p^{2}} satisfying Assumption 2, \rho({\mbox{\boldmathE}}_{UV})=O(\sqrt{k/p})\cdot\rho\left({\mbox{\boldmathE}}\right), where {\mbox{\boldmathE}}_{UV}=P_{V\otimes U}{\mbox{\boldmathE}}P_{V\otimes U}.

Proof.

Suppose {\mbox{\boldmathE}}=\sum_{i=1}^{p^{2}}\lambda_{i}{\mbox{\boldmathw}}_{i}{\mbox{\boldmathw}}_{i}^{\top} and \lVert{\mbox{\boldmathw}}_{i}\rVert^{2}=1. For $1\leq i\leq p^{2}$, let {\mbox{\boldmathW}}_{i}={\textsc{mat}}\left({\mbox{\boldmathw}}_{i}\right) be the $i$-th eigenmatrices of $E$. Then

[TABLE]

Then \lVert P_{V\otimes U}{\mbox{\boldmathw}}_{i}\rVert_{\ell_{2}}=O(\sqrt{k/p}).

Let {\mbox{\boldmathW}}=\begin{pmatrix}{\mbox{\boldmathw}}_{1},\cdots,{\mbox{\boldmathw}}_{p^{2}}\end{pmatrix}, {\mbox{\boldmathE}}={\mbox{\boldmathW}}{\mbox{\boldmath\Lambda}}{\mbox{\boldmathW}}^{\top} and {\mbox{\boldmathQ}}=P_{V\otimes U}{\mbox{\boldmathW}}, we have \lVert{\mbox{\boldmathQ}}\rVert_{F}^{2}=\lVert{\mbox{\boldmathQ}}{\mbox{\boldmathQ}}^{\top}\rVert_{\star}=O(kp), \lVert{\mbox{\boldmathQ}}\rVert^{2}=O(k/p) and P_{V\otimes U}{\mbox{\boldmathE}}{\mbox{\boldmathW}}={\mbox{\boldmathQ}}{\mbox{\boldmath\Lambda}}. Thus,

[TABLE]

Thus, we have

[TABLE]

∎

We now state our main result as below, which shows that under appropriate conditions, the smart-pm algorithm can recover the low rank eigenmatrix.

Theorem 1**.**

Assume that Assumption 1 and 2 hold. Let $1\leq\bar{k}\leq k$. Assume that \kappa\coloneqq\triangle\lambda/\rho({\mbox{\boldmathE}}_{UV})>2 where {\mbox{\boldmathE}}_{UV} is defined in Lemma 1 for any rank-$k$ matrices $U$ and $V$. Define

[TABLE]

and

[TABLE]

If |\langle{\mbox{\boldmathX}}_{0},{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}\rangle|\geq\theta+\delta for some {\mbox{\boldmathX}}_{0} with {\rm rank}\left({\mbox{\boldmathX}}_{0}\right)=k and \lVert{\mbox{\boldmathX}}_{0}\rVert_{F}=1 and $\theta\in(0,1)$ such that

[TABLE]

Then we either have

[TABLE]

or for all $t\geq 0$

[TABLE]

Proof.

We sketch the proof here, while the details are relegated to Appendix A in supplementary material. The proof is carried out in the following steps.

Lemma 2 establishes the perturbation theory of symmetric eigenvalue problem for rank truncated {\mbox{\boldmathA}}_{UV}={}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}}_{UV}+{\mbox{\boldmathE}}_{UV}, i.e. bounds \lVert{\mbox{\boldmathx}}_{1}({\mbox{\boldmathA}}_{UV})-\bar{\mbox{\boldmathx}}\rVert_{\ell_{2}}, where {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}}_{UV}=\lambda_{1}\bar{\mbox{\boldmathx}}\bar{\mbox{\boldmathx}}^{\top} and {\mbox{\boldmathx}}_{1}({\mbox{\boldmathA}}_{UV}) is the eigenvector corresponding to the largest eigenvalue of {\mbox{\boldmathA}}_{UV}. 2. 2.

Lemma 4 quantifies the error introduced by the rank-truncation step in the SMART-PM, i.e. |\langle RankTrunc({\mbox{\boldmathX}}^{\prime}_{t},k),{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}\rangle|-|\langle{\mbox{\boldmathX}}^{\prime},{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}\rangle|. 3. 3.

Lemma 5 establishes that each step of the SMART-PM improves eigenvector estimation, i.e. \sqrt{1-|\langle{\mbox{\boldmathX}}_{t},{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}\rangle|} decreases geometrically with factor $\mu$. 4. 4.

Based on results from Lemma 2, 4, and 5, this step is similar to that of Theorem 4 in Yuan and Zhang (2013) and thus omitted in the Appendix A.

∎

Remark 2**.**

For any fixed eigen-gap $\triangle\lambda$, if \rho({\mbox{\boldmathE}}_{UV})<\triangle\lambda/2, then $\gamma<1$ and \delta=O(1/\kappa)=O(\rho({\mbox{\boldmathE}}_{UV})). If $\gamma$ is sufficiently small, then the requirement that $\mu<1$ can be satisfied for a sufficiently small $\theta$ of the order $\sqrt{\bar{k}/k}$. Theorem 1 shows that under appropriate conditions, as long as we can find an initial {\mbox{\boldmathX}}_{0} such that

[TABLE]

for some constant $C$, then \sqrt{1-|\langle{\mbox{\boldmathX}}_{t},{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}\rangle|} converges geometrically until

[TABLE]

Theorem 1 thus provides strong theoretical justification of the proposed smart-pm algorithm. Specifically, the replacement of the full matrix perturbation error \rho({\mbox{\boldmathE}}) with \sqrt{k/p}\,\rho({\mbox{\boldmathE}}) gives theoretical insights on why SMART-PM gives superior results in Section 4. To illustrate this further, we briefly describe a consequence of Theorem 1 in the case, when {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}} is corrupted by an iid standard Gaussian noise matrix $E$. In this case, \rho({\mbox{\boldmathE}})=O(\sqrt{p^{2}}) with high-probability. The spectral norm of the low rank projection \rho(P_{V\otimes U}{\mbox{\boldmathE}}P_{V\otimes U})=O(\sqrt{pk}), which grows linearly in $\sqrt{pk}$ instead of $p$. Hence, we can run smart-pm with an appropriate initial vector to obtain an approximate solution {\mbox{\boldmathX}}_{t} with error \|{\mbox{\boldmathX}}_{t}-{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}\|_{F}=O\left(\sqrt{pk}\right). Note that for estimating a $p\times p$ matrix with rank $k$, a simple parameter counting argument shows that the minimax lower bound is of the order $O\left(\sqrt{pk}\right)$. Thus the smart-pm algorithm achieves the minimax lower bound provided an initializer satisfying the conditions of Theorem 1 is used. The question of obtaining such an initializer is more delicate. In the case of sparse PCA, an interesting statistical-computational trade-off exists; see for example Berthet and Rigollet (2013); Ma and Wigderson (2015); Ma and Wu (2015); Wang et al. (2016); Brennan et al. (2018). We conjecture that such a phenomenon exists in our problem as well.

Finally, we remark that while the steps of the proof of Theorem 1 are motivated by the proof of Theorem 4 in Yuan and Zhang (2013) on the analysis of truncated power method for sparse PCA, the details of our proof are very different because of our matricized low-rank structure on the eigenvectors and the involved non-commutativity issues involved. Also, note that Lemma 1 qualifies the largest absolute eigenvalue of a low rank projection of noise matrix \rho\left({\mbox{\boldmathE}}_{UV}\right)=O(\sqrt{k/p})\cdot\rho\left({\mbox{\boldmathE}}\right), which gives a different rate from that of the sparse case.

4 Experiments

In this section, we show simulation results that confirm the relevance of our theoretical findings regarding smart-pm algorithm. The data \{{\mbox{\boldmathy}}_{1},\ldots,{\mbox{\boldmathy}}_{n}\} are generated from \mathcal{N}(0,{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}}) with true covariance

[TABLE]

where \lVert{\mbox{\boldmath\Sigma}}_{\epsilon}\rVert=\lambda_{2} is set to $1$ throughout experiments. The empirical covariance is

[TABLE]

Theorem 1 implies that under appropriate conditions, the estimation error \|{\mbox{\boldmathX}}_{t}-{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}\|_{F} is proportional to $\delta$ which is an increasing function of \rho({\mbox{\boldmathE}}_{UV}) and decreasing function of eigen-gap $\triangle\lambda=\lambda_{1}-\lambda_{2}$. For this model, we have \rho({\mbox{\boldmathE}}_{UV})=O(\sqrt{pk/n}). Therefore, for fixed dimensionality $p^{2}$, the error bound is controlled by the triplet $\{n,k,\lambda_{1}-\lambda_{2}\}$.

In this study, we consider a setup with dimensionality $d=p^{2}=32^{2}=1024$, and the eigenmatrix {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}={\textsc{mat}}\left(\bar{\mbox{\boldmathx}}\right) a rank $\bar{k}$ uniform random matrix with $\bar{k}=1$ and \lVert{\textsc{mat}}\left({}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}\right)\rVert_{F}=1. We compare the convergence trajectory, sample efficiency and effect of eigen-gap between smart-pm and vanilla power method. We also show how different value of $k$ affects the performance of smart-pm.

4.1 Convergence Trajectory

Figure 2(a), 2(b) and 2(c) show the convergence trajectory for smart-pm and Power method for $n=800$, $k=2$, and different $\lambda_{1}\in\{5,10,100\}$. The y-axis corresponds to the log of estimation error \|{\mbox{\boldmathX}}_{t}-{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}\|_{F} and the x-axis corresponds to the number of iterations. We tried three different methods of initialization: (1) Randomly generated {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}_{0}; (2) Rank-$k$ randomly generated {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}_{0}; (3) Maximum eigenmatrix of $A$. It can be seen from Figure 2 that all methods converge faster when the eigen-gap is larger. smart-pm with maximum eigenmatrix of $A$ as initial value (smart-pm 3) converges faster than the other methods. smart-pm with random initialization (smart-pm 1) converge faster than simple power method under all settings. However, random rank-$k$ initialization (smart-pm 2) does not work as well as random initialization (smart-pm 1) perhaps due to $k$ being set at 2, while $\bar{k}=1$.

4.2 Sample Efficiency and Eigen-Gap Effect

In this case, we work with $n\in\{100,200,400,800,1600\}$, $\lambda_{1}\in\{5,10,100\}$ and $\lambda_{2}=1$ given $k=2\geq\bar{k}$. $k$ can also be tuned and selected using cross-validation as in Yuan and Zhang (2013). For each pair $\{n,\lambda_{1}\}$, we generate $100$ data sets, calculate the empirical covariance matrices and employ smart-pm and Power method to compute a rank-$k$ eigenmatrix \widehat{\mbox{\boldmathX}}.

Figure 3(a), 3(b) and 3(c) shows the log of estimation curves as functions of sample size $n$ increases with $\lambda_{1}\in\{5,10,100\}$, respectively. For fixed eigen-gap, the estimation error decreases as sample size $n$ increases. For fixed sample size, the estimation error is smaller with larger eigen-gap. Again, estimation error by smart-pm 3 with maximum eigenmatrix of $A$ as initiate {\mbox{\boldmathX}}_{0} has smaller mean and variance. smart-pm 1 with random initialization performs as well as smart-pm 3 under almost all setting except in the beginning of Figure 3(a) when both sample size and eigen-gap are very small. smart-pm 2 with rank-$k$ random initialization performs well when sample size or eigen-gap is large enough but may gets unstable with small eigen-gap and sample size. In all settings, smart-pm outperforms the Power method.

4.3 Varying Input Rank $k$

In this case, we fix sample size $n=100$, $\lambda_{1}=100$ and $\lambda_{2}=1$, and test the values of $k\in\{1,2,4,8,16,32\}$. We generate $100$ empirical covariance matrices and employ smart-pm and power method to compute a rank-$k$ eigenmatrix. From experiments in the previous sections, the random initialization of {\mbox{\boldmathX}}_{0} performs well. Hence, we employ the random initialization here. Recall that $p=32$ in our setting. When $k=p=32$, smart-pm degenerates to the power method. Figure 4 shows the log of estimation error curve as function of $k$. The power method is not relevant to $k$ so the corresponding curve is a horizontal line. The error curve of smart-pm stays below while approaching that of the power method as $k$ increases. Also, the variance of the error becomes larger as input rank $k$ increases. Similar observations are also made for other fixed pairs $\{n,\lambda_{2}-\lambda_{1}\}$.

5 Discussion

In this paper, we propose a novel Low-rank Principal Eigenmatrix Analysis technique for data analysis which is formulated as a non-convex optimization problem. We propose the smart-pc algorithm for solving the above formulation and established its computational and statistical properties.

Our method opens up a lot of potential directions for follow-up work. First, the question of obtaining an initializer satisfying the condition of Theorem 1 or proving the nonexistence of a obtaining such an initializer in polynomial time is extremely interesting. It is also interesting to develop convex relaxations for the proposed non-convex problem that works provably well. Next, analyzing the general permuted low-rank eigenmatrix model, discussed in §2.1 is challenging. Finally, note that in this work we only considered matricization of the top-eigenvector. One could potentially also consider tensorizations of the top-eigenvector and enforce structures on that. Indeed, we also noted empirically that the eigenvectors when tensorized also have a certain low-rank property. We plan to investigate these directions for future work.

Appendix A Proof of Theorem 1

This section provides detailed proof of Theorem 1. It is carried out in the following steps: Lemma 2 establishes the perturbation theory of symmetric eigenvalue problem for rank truncated {\mbox{\boldmathA}}_{UV}={}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}}_{UV}+{\mbox{\boldmathE}}_{UV}, i.e. bounds \lVert{\mbox{\boldmathx}}_{1}({\mbox{\boldmathA}}_{UV})-\bar{\mbox{\boldmathx}}\rVert_{\ell_{2}}, where {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}}_{UV}=\lambda_{1}\bar{\mbox{\boldmathx}}\bar{\mbox{\boldmathx}}^{\top} and {\mbox{\boldmathx}}_{1}({\mbox{\boldmathA}}_{UV}) is the eigenvector corresponding to the largest eigenvalue of {\mbox{\boldmathA}}_{UV}. Lemma 3 conducts convergence analysis of traditional power method. Lemma 4 quantifies the error introduced by the rank-truncation step in the smart-pm, i.e. |\langle RankTrunc({\mbox{\boldmathX}}^{\prime}_{t},k),{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}\rangle|-|\langle{\mbox{\boldmathX}}^{\prime},{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}\rangle|. Lemma 5 establishes that each step of the smart-pm improves eigenvector estimation, i.e. \sqrt{1-|\langle{\mbox{\boldmathX}}_{t},{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}\rangle|} decreases geometrically with factor $\mu$. In the proof, we denote the pair of bold-faced ({\mbox{\boldmathx}}_{t},{\mbox{\boldmathX}}_{t}) as the pair of vectorization and matricization and use frequently the fact that {\mbox{\boldmathx}}_{t}^{\top}\bar{\mbox{\boldmathx}}=\langle{\mbox{\boldmathX}}_{t},{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}\rangle.

A.1 The Perturbation Theorem of Rank-Truncated Symmetric Eigenvalue Problem

In this section, main lemma 2 states the perturbation theory of symmetric eigenvalue problem for rank truncated {\mbox{\boldmathA}}={}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}}+{\mbox{\boldmathE}}.

Lemma 2**.**

Consider projection $P_{U}$ and $P_{V}$ with P_{U}{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}P_{V}={}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}} where ${\rm rank}\left(P_{U}\right)={\rm rank}\left(P_{V}\right)=k$ and $1\leq k\leq p$. Under Assumption 1 and 2, if \kappa\coloneqq\triangle\lambda/\rho({\mbox{\boldmathE}}_{UV})>2, then the ratio of the second largest (in absolute value) to the largest eigenvalue of matrix {\mbox{\boldmathA}}_{UV} is no more than

[TABLE]

Moreover,

[TABLE]

Proof.

Define {\mbox{\boldmathA}}_{UV}\coloneqq P_{V\otimes U}{\mbox{\boldmathA}}P_{V\otimes U}, {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}}_{UV}\coloneqq P_{V\otimes U}{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}}P_{V\otimes U}, and {\mbox{\boldmathE}}_{UV}\coloneqq P_{V\otimes U}{\mbox{\boldmathE}}P_{V\otimes U}. Using Weyl’s inequality, we obtain

[TABLE]

and $\forall j\geq 2$,

[TABLE]

Then,

[TABLE]

Now write the eigenvector w.r.t the largest eigenvalue \lambda_{1}({\mbox{\boldmathA}}_{UV}) as {\mbox{\boldmathx}}_{1}({\mbox{\boldmathA}}_{UV})=\alpha\bar{\mbox{\boldmathx}}+\beta{\mbox{\boldmathz}}, where \lVert\bar{\mbox{\boldmathx}}\rVert_{\ell_{2}}=\lVert{\mbox{\boldmathz}}\rVert_{\ell_{2}}=1, \bar{\mbox{\boldmathx}}^{\top}{\mbox{\boldmathz}}=0 and $\alpha^{2}+\beta^{2}=1$. This implies that

[TABLE]

Then right multiplying the above equation by {\mbox{\boldmathz}}^{\top}, we have

[TABLE]

This leads to,

[TABLE]

where t\coloneqq{\rho\left({\mbox{\boldmathE}}_{UV}\right)}/{(\triangle\lambda-2\rho\left({\mbox{\boldmathE}}_{UV}\right))} under the condition that \triangle\lambda>2\rho\left({\mbox{\boldmathE}}_{UV}\right). Then we have $1=\alpha^{2}+\beta^{2}\leq\alpha^{2}(1+t^{2})$ and $\alpha^{2}\geq\frac{1}{1+t^{2}}$. Without loss of generality, we may assume that $\alpha>0$ because otherwise we can replace \bar{\mbox{\boldmathx}} with -\bar{\mbox{\boldmathx}}. It follows that

[TABLE]

We have

[TABLE]

∎

A.2 The Error Analysis of SMART-PM

The following lemma from Lemma 11 in Yuan and Zhang (2013) conducts convergence analysis of traditional power method. We restate it here for completeness.

Lemma 3**.**

Let \bar{\mbox{\boldmathx}} be the eigenvector with the largest (in absolute value) eigenvalue of a symmetric matrix $A$, and let $\gamma<1$ be the ratio of the second largest to largest eigenvalue in absolute values. Given any $y$ such that \lVert{\mbox{\boldmathy}}\rVert_{\ell_{2}}=1 and \bar{\mbox{\boldmathx}}^{\top}{\mbox{\boldmathy}}>0; Let {\mbox{\boldmathy}}^{\prime}={\mbox{\boldmathA}}{\mbox{\boldmathy}}/\lVert{\mbox{\boldmathA}}{\mbox{\boldmathy}}\rVert_{\ell_{2}}, then

[TABLE]

The following lemma quantifies the error introduced by the rank-truncation step in smart-pm algorithm.

Lemma 4**.**

Consider \bar{\mbox{\boldmathx}} with {\rm rank}\left({\textsc{mat}}\left(\bar{\mbox{\boldmathx}}\right)\right)=\bar{k}. For {\mbox{\boldmathy}}\in\mathbb{R}^{p^{2}} and $1\leq k<p$, let the singular value decomposition of {\mbox{\boldmathY}}={\textsc{mat}}\left({\mbox{\boldmathy}}\right) be {\mbox{\boldmathU}}{\mbox{\boldmathD}}{\mbox{\boldmathV}}^{T} where the diagonal of $D$ contains singular values in (absolute) non-increasing order. Define Matricized Rank-Truncation operator RankTrunc({\textsc{mat}}\left({\mbox{\boldmathy}}\right),k)=P_{U_{k}}{\mbox{\boldmathY}}P_{V_{k}} and RankTrunc({\mbox{\boldmathy}},k)={\textsc{vec}}\left(RankTrunc({\textsc{mat}}\left({\mbox{\boldmathy}}\right),k)\right), where $P_{U_{k}}$ and $P_{V_{k}}$ are projection matrix onto the first k columns of $U$ and $V$,respectively. If \|\bar{\mbox{\boldmathx}}\|_{\ell_{2}}=\|{\mbox{\boldmathy}}\|_{\ell_{2}}=1, then

[TABLE]

Proof.

Define {\mbox{\boldmathY}}={\textsc{mat}}\left({\mbox{\boldmathy}}\right), {\mbox{\boldmathy}}={\textsc{vec}}\left({\mbox{\boldmathY}}\right) and similar for \bar{\mbox{\boldmathx}} and \bar{\mbox{\boldmathX}}. Suppose the SVD decomposition of $Y$ (\bar{\mbox{\boldmathX}}) assume the form {\mbox{\boldmathU}}{\mbox{\boldmathD}}{\mbox{\boldmathV}}^{\top} ({\mbox{\boldmathP}}{\mbox{\boldmath\Lambda}}{\mbox{\boldmathQ}}^{\top}) where $U$, $V$, $P$ and $Q$ are orthonormal matrices and {\mbox{\boldmathD}}={\rm diag}\left(d_{1},\ldots,d_{p}\right) and {\mbox{\boldmath\Lambda}}={\rm diag}\left(\lambda_{1},\ldots,\lambda_{k}^{\star}\right) with singular values in non-increasing order. For simplicy, we work with non-negative singular value because otherwise we may just change the sign of one singular vector. Then,

[TABLE]

We also have,

[TABLE]

Define $P_{1}=P_{P_{\bar{k}}\otimes Q_{\bar{k}}}P_{V_{k}\otimes U_{k}}^{\perp}$, $P_{2}=P_{P_{\bar{k}}\otimes Q_{\bar{k}}}P_{V_{k}\otimes U_{k}}$, $P_{3}=P_{P_{\bar{k}}\otimes Q_{\bar{k}}}^{\perp}P_{V_{k}\otimes U_{k}}$ and $P_{4}=P_{P_{\bar{k}}\otimes Q_{\bar{k}}}^{\perp}P_{V_{k}\otimes U_{k}}^{\perp}$. Let $k_{i}={\rm rank}\left(P_{i}\right)$. Since \|\bar{\mbox{\boldmathx}}\|_{\ell_{2}}=\|{\mbox{\boldmathy}}\|_{\ell_{2}}=1 and {\rm rank}\left(\bar{\mbox{\boldmathx}}\right)=\bar{k}, we have \sum_{i=1}^{4}\lVert P_{i}{\mbox{\boldmathy}}\rVert^{2}_{\ell_{2}}=1 and \sum_{i=1}^{2}\lVert P_{i}{\mbox{\boldmathx}}\rVert^{2}_{\ell_{2}}=1. Then

[TABLE]

where the first equation follows from the fact that $P_{i}^{\top}P_{j}=0$ for $i\neq j$ and P_{3}\bar{\mbox{\boldmathx}}=P_{4}\bar{\mbox{\boldmathx}}=0. The first inequality follows from Cauchy-Schwarz inequality and the third inequality follow from 2\lVert P_{1}{\mbox{\boldmathy}}\rVert_{\ell_{2}}\lVert P_{1}\bar{\mbox{\boldmathx}}\rVert_{\ell_{2}}\lVert P_{2}{\mbox{\boldmathy}}\rVert_{\ell_{2}}\lVert P_{2}\bar{\mbox{\boldmathx}}\rVert_{\ell_{2}}\leq\lVert P_{1}{\mbox{\boldmathy}}\rVert_{\ell_{2}}^{2}\lVert P_{2}\bar{\mbox{\boldmathx}}\rVert^{2}_{\ell_{2}}+\lVert P_{2}{\mbox{\boldmathy}}\rVert^{2}_{\ell_{2}}\lVert P_{1}\bar{\mbox{\boldmathx}}\rVert^{2}_{\ell_{2}} and \lVert P_{1}\bar{\mbox{\boldmathx}}\rVert^{2}_{\ell_{2}}+\lVert P_{2}\bar{\mbox{\boldmathx}}\rVert^{2}_{\ell_{2}}=1. The last inequality comes from \lVert P_{3}{\mbox{\boldmathy}}\rVert^{2}_{\ell_{2}}/k_{3}\geq\lVert P_{1}{\mbox{\boldmathy}}\rVert^{2}_{\ell_{2}}/k_{1}, since $d_{1}\geq\cdots\geq d_{p}$.

This implies that

[TABLE]

where the second inequality follows from $\bar{k}\leq k$. Then we obtain that

[TABLE]

Finally from (8) and (9),

[TABLE]

∎

Next lemma says that each step of the smart-pm algorithm improves eigenvector estimation.

Lemma 5**.**

Assume that $k\geq\bar{k}$. If |\langle{\mbox{\boldmathX}}_{t-1},{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathX}}}\rangle|>\theta+\delta for $\theta\in(0,1)$, then the iterate {\mbox{\boldmathx}}_{t} and {\mbox{\boldmathx}}_{t-1} satisfy the following inequality,

[TABLE]

where $\gamma$ and $\delta$ are defined in (3) and (4), respectively, and

[TABLE]

Proof.

Let {\mbox{\boldmathU}}={\mbox{\boldmathU}}_{t-1}\cup{\mbox{\boldmathU}}_{t}\cup\bar{\mbox{\boldmathU}} and {\mbox{\boldmathV}}={\mbox{\boldmathV}}_{t-1}\cup{\mbox{\boldmathV}}_{t}\cup\bar{\mbox{\boldmathV}}. Consider the following vector

[TABLE]

where {\mbox{\boldmathA}}_{UV}\coloneqq P_{V\otimes U}{}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}}P_{V\otimes U} denotes the matrix after project the rows and columns of {}\mkern 3.0mu\overline{\mkern-3.0mu{\mbox{\boldmathA}}} on the spaces spanned by {\mbox{\boldmathV}}\otimes{\mbox{\boldmathU}}, respectively. We note that replacing {\mbox{\boldmathX}}_{t} with {\textsc{mat}}\left({\mbox{\boldmathy}}_{t}\right) in Algorithm 1 does not affect the output iteration sequence \{{\mbox{\boldmathX}}_{t}\} because of the low-rankness of {\mbox{\boldmathX}}_{t-1} and the fact that the truncation operation is invariant to scaling. Therefore for notation simplicity, in the following proof we will simply assume that {\mbox{\boldmathX}}_{t} is redefined as {\mbox{\boldmathX}}_{t}={\textsc{mat}}\left({\mbox{\boldmathy}}_{t}\right) according to (10).

Let {\mbox{\boldmathx}}_{1}({\mbox{\boldmathA}}_{UV}) be the eigenvector with the largest (in absolute value) eigenvalue of {\mbox{\boldmathA}}_{UV}. Without loss of generality and for simplicity, we may assume that {\mbox{\boldmathy}}_{t}{\mbox{\boldmathx}}_{1}({\mbox{\boldmathA}}_{UV})\geq 0 and {\mbox{\boldmathy}}_{t-1}\bar{\mbox{\boldmathx}}\geq 0 because otherwise we can simply do appropriate sign changes in the proof. From Lemma 3, we have

[TABLE]

By assumption that |{\mbox{\boldmathx}}_{t-1}^{\top}\bar{\mbox{\boldmathx}}|>\theta+\delta and Lemma 2, we have

[TABLE]

Thus we have

[TABLE]

where $\gamma$ is defined in Lemma 2. Using the fact that \lVert{\mbox{\boldmathy}}_{t}\rVert_{\ell_{2}}=\lVert{\mbox{\boldmathx}}_{1}({\mbox{\boldmathA}}_{UV})\rVert_{\ell_{2}}=1, the above result is equivalent to

[TABLE]

Using Lemma 2,

[TABLE]

This is equivalent to

[TABLE]

Applying Lemma 4 and using $k>\bar{k}$, if (1+(\bar{k}/k)^{1/2})\sqrt{1-({\mbox{\boldmathy}}_{t}^{\top}\bar{\mbox{\boldmathx}})^{2}}\leq 1, we have

[TABLE]

where

[TABLE]

∎