An $\ell_{\infty}$ Eigenvector Perturbation Bound and Its Application to Robust Covariance Estimation

Introduction

The perturbation of matrix eigenvectors (or singular vectors) has been well studied in matrix perturbation theory (Wedin, 1972; Stewart, 1990). The best known result of eigenvector perturbation is the classic Davis-Kahan theorem (Davis and Kahan, 1970). It originally emerged as a powerful tool in numerical analysis, but soon found its widespread use in other fields, such as statistics and machine learning. Its popularity continues to surge in recent years, which is largely attributed to the omnipresent data analysis, where it is a common practice, for example, to employ PCA (Jolliffe, 2002) for dimension reduction, feature extraction, and data visualization.

The eigenvectors of matrices are closely related to the underlying structure in a variety of problems. For instance, principal components often capture most information of data and extract the latent factors that drive the correlation structure of the data (Bartholomew et al., 2011); in classical multidimensional scaling (MDS), the centered squared distance matrix encodes the coordinates of data points embedded in a low dimensional subspace (Borg and Groenen, 2005); and in clustering and network analysis, spectral algorithms are used to reveal clusters and community structure (Ng et al., 2002; Rohe et al., 2011). In those problems, the low dimensional structure that we want to recover, is often ‘perturbed’ by observation uncertainty or statistical errors. Besides, there might be a sparse pattern corrupting the low dimensional structure, as in approximate factor models (Chamberlain and Rothschild, 1982; Stock and Watson, 2002) and robust PCA (De La Torre and Black, 2003; Candès et al., 2011).

A general way to study these problems is to consider

where $A$ is a low rank matrix, $S$ is a sparse matrix, and $N$ is a random matrix regarded as random noise or estimation error, all of which have the same size $d_{1}\times d_{2}$. Usually $A$ is regarded as the ‘signal’ matrix we are primarily interested in, $S$ is some sparse contamination whose effect we want to separate from $A$, and $N$ is the noise (or estimation error in covariance matrix estimation).

The decomposition (1) forms the core of a flourishing literature on robust PCA (Chandrasekaran et al., 2011; Candès et al., 2011), structured covariance estimation (Fan et al., 2008, 2013), multivariate regression (Yuan et al., 2007) and so on. Among these works, a standard condition on $A$ is matrix incoherence (Candès et al., 2011). Let the singular value decomposition be

where $U_{ij}$ and $V_{ij}$ are the $(i,j)$ entry of $U$ and $V$, respectively. It is usually expected that $\mu_{0}:=\max\{\mu(U),\mu(V)\}$ is not too large, which means the singular vectors $u_{i}$ and $v_{i}$ are incoherent with the standard basis. This incoherence condition (3) is necessary for us to separate the sparse component $S$ from the low rank component $A$; otherwise $A$ and $S$ are not identifiable. Note that we do not need any incoherence condition on $UV^{T}$, which is different from Candès et al. (2011) and is arguably unnecessary (Chen, 2015).

Now we denote the eigengap $\gamma_{0}=\min\{\sigma_{i}-\sigma_{i+1}:i=1,\ldots,r\}$ where $\sigma_{r+1}:=0$ for notational convenience. Also we let $E=S+N$, and view it as a perturbation matrix to the matrix $A$ in (1). To quantify the perturbation, we define a rescaled measure as $\tau_{0}:=\max\{\sqrt{d_{2}/d_{1}}\|E\|_{1},\sqrt{d_{1}/d_{2}}\|E\|_{\infty}\}$, where

which are commonly used norms gauging sparsity (Bickel and Levina, 2008). They are also operator norms in suitable spaces (see Section 2). The rescaled norms $\sqrt{d_{2}/d_{1}}\|E\|_{1}$ and $\sqrt{d_{1}/d_{2}}\|E\|_{\infty}$ are comparable to the spectral norm $\|E\|_{2}:=\max_{\|u\|_{2}=1}\|Eu\|_{2}$ in many cases; for example, when $E$ is an all-one matrix, $\sqrt{d_{2}/d_{1}}\|E\|_{1}=\sqrt{d_{1}/d_{2}}\|E\|_{\infty}=\|E\|_{2}$.

Suppose the perturbed matrix $\widetilde{A}$ also has the singular value decomposition:

where $\widetilde{\sigma}_{i}$ are nonnegative and in the decreasing order, and the notation $\wedge$ means $a\wedge b=\min\{a,b\}$. Denote $\widetilde{U}=[\widetilde{u}_{1},\ldots,\widetilde{u}_{r}],V=[\widetilde{v}_{1},\ldots,\widetilde{v}_{r}]$, which are counterparts of top $r$ singular vectors of $A$.

Let $\widetilde{A}=A+E$ and suppose the singular decomposition in (2) and (5). Denote $\gamma_{0}=\min\{\sigma_{i}-\sigma_{i+1}:i=1,\ldots,r\}$ where $\sigma_{r+1}:=0$. Then there exists $C(r,\mu_{0})=O(r^{4}\mu_{0}^{2})$ such that, if $\gamma_{0}>C(r,\mu_{0})\tau_{0}$, up to sign,

where $\mu_{0}=\max\{\mu(U),\mu(V)\}$ is the coherence given after (3) and $\tau_{0}:=\max\{\sqrt{d_{2}/d_{1}}\|E\|_{1},\sqrt{d_{1}/d_{2}}\|E\|_{\infty}\}$.

When $A$ is symmetric, the condition on the eigengap is simply $\gamma_{0}>C(r,\mu_{0})\|E\|_{\infty}$. It naturally holds for a variety of applications, where the low rank structure emerges as a consequence of a few factors driving the data matrix. For example, in Fama-French factor models, the excess returns in a stock market are driven by a few common factors (Fama and French, 1993); in collaborative filtering, the ratings of users are mostly determined by a few common preferences (Rennie and Srebro, 2005); in video surveillance, $A$ is associated with the stationary background across image frames (Oliver et al., 2000). We will have a detailed discussion in Section 2.3.

The eigenvector perturbation was studied by Davis and Kahan (1970), where Hermitian matrices were considered, and the results were extended by Wedin (1972) to general rectangular matrices. To compare our result with these classical results, assuming $\gamma_{0}\geq 2\|E\|_{2}$, a combination of Wedin’s theorem and Mirsky’s inequality (Mirsky, 1960) (the counterpart of Weyl’s inequality for singular values) implies

To understand how matrix incoherence helps, let us consider a simple example with no matrix incoherence, in which (7) is tight up to a constant. Let $A=d(1,0,\ldots,0)^{T}(1,0,\ldots,0)$ be a $d$-dimensional square matrix, and $E=d(0,1/2,0,\ldots,0)^{T}(1,0,\ldots,0)$ of the same size. It is apparent that $\gamma_{0}=d,\tau_{0}=d/2$, and that $v_{1}=(1,0,\ldots,0)^{T},\widetilde{v}_{1}=(2/\sqrt{5},1/\sqrt{5},0,\ldots,0)^{T}$ up to sign. Clearly, the perturbation $\|\widetilde{v}_{1}-v_{1}\|_{\infty}$ is not vanishing as $d$ tends to infinity in this example, and thus, there is no hope of a strong upper bound as in (6) without the incoherence condition.

Our result is very different from the sparse PCA literature, in which it is usually assumed that the leading eigenvectors are sparse. In Johnstone and Lu (2009), it is proved that there is a threshold for $p/n$ (the ratio between the dimension and the sample size), above which PCA performs poorly, in the sense that $\langle\widetilde{v}_{1},v_{1}\rangle$ is approximately . This means that the principal component computed from the sample covariance matrix reveals nothing about the true eigenvector. In order to mitigate this issue, in Johnstone and Lu (2009) and subsequent papers (Vu and Lei, 2012; Ma, 2013; Berthet and Rigollet, 2013), sparse leading eigenvectors are assumed. However, our result is different, in the sense that we require a stronger eigengap condition $\gamma_{0}>C(r,\mu_{0})\|E\|_{\infty}$ (i.e. stronger signal), whereas in Johnstone and Lu (2009), the eigengap of the leading eigenvectors is a constant times $\|E\|_{2}$. This explains why it is plausible to have a strong uniform eigenvector perturbation bound in this paper.

We will illustrate the power of this perturbation result using robust covariance estimation as one application. In the approximate factor model, the true covariance matrix admits a decomposition into a low rank part $A$ and a sparse part $S$. Such models have been widely applied in finance, economics, genomics, and health to explore correlation structure.

However, in many studies, especially financial and genomics applications, it is well known that the observations exhibit heavy tails (Gupta et al., 2013). This problem can be resolved with the aid of recent results of concentration bounds in robust estimation (Catoni, 2012; Hsu and Sabato, 2014; Fan et al., 2017), which produces the estimation error $N$ in (1) with an optimal entry-wise bound. It nicely fits our perturbation result, and we can tackle it easily by following the ideas in Fan et al. (2013).

where $\Lambda_{1}=\text{diag}\{\lambda_{1},\ldots,\lambda_{r}\}$, $\Lambda_{2}=\mbox{diag}\{\lambda_{r+1},\ldots,\lambda_{n}\}$, and where $|\lambda_{1}|\geq|\lambda_{2}|\geq\ldots\geq|\lambda_{n}|$. Note the best rank-$r$ approximation of $A$ under the Frobenius norm is $A_{r}:=\sum_{i\leq r}\lambda_{i}v_{i}v_{i}^{T}$.This is a consequence of Wielandt-Hoffman theorem. Analogously, the spectral decomposition of $\widetilde{A}$ is

We will use notations $O(\cdot)$ and $\Omega(\cdot)$ to hide absolute constants.We write $a=O(b)$ if there is a constant $C>0$ such that $a<Cb$; and $a=\Omega(b)$ if there is a constant $C^{\prime}>0$ such that $a>C^{\prime}b$. The next theorem bounds the perturbation of eigenspaces up to a rotation.

This result involves an unspecified rotation $R$, due to the possible presence of multiplicity of eigenvalues. In the case where $\lambda_{1}=\cdots=\lambda_{r}>0$, the individual eigenvectors of $V$ are only identifiable up to rotation. However, assuming an eigengap (similar to Davis-Kahan theorem), we are able to bound the perturbation of individual eigenvectors (up to sign).

Assume the conditions in Theorem 2.1. In addition, suppose $\delta$ satisfies $\delta>\|E\|_{2}$, and for any $i\in[r]$, the interval $[\lambda_{i}-\delta,\lambda_{i}+\delta]$ does not contain any eigenvalues of $A$ other than $\lambda_{i}$. Then, up to sign,

To understand the above two theorems, let us consider the case where $A$ has exactly rank $r$ (i.e., $\varepsilon=0$), and $r$ and $\mu$ are not large (say, bounded by a constant). Theorem 2.1 gives a uniform entrywise bound $O(\|E\|_{\infty}/|\lambda_{r}|\sqrt{d})$ on the eigenvector perturbation. As a comparison, the Davis–Kahan $\sin\Theta$ theorem (Davis and Kahan, 1970) gives a bound $O(\|E\|_{2}/|\lambda_{r}|)$ on $\|\widetilde{V}R-V\|_{2}$ with suitably chosen rotation $R$.To see how the Davis-Kahan $\sin\Theta$ theorem relates to this form, we can use the identity $\|\sin\Theta(\widetilde{V},V)\|_{2}=\|\widetilde{V}\widetilde{V}^{T}-VV^{T}\|_{2}$ (Stewart, 1990), and the (easily verifiable) inequality $2\min_{R}\|\widetilde{V}R-V\|_{2}\geq\|\widetilde{V}\widetilde{V}^{T}-VV^{T}\|_{2}\geq\min_{R}\|\widetilde{V}R-V\|_{2}$ where $R$ is an orthogonal matrix. This is an order of $\sqrt{d}$ larger than the bound given in Theorem 2.1 when $\|E\|_{\infty}$ is of the same order as $\|E\|_{2}$. Thus, in scenarios where $\|E\|_{2}$ is comparable to $\|E\|_{\infty}$, this is a refinement of Davis-Kahan theorem, because the max-norm bound in Theorem 2.1 provides an entry-wise control of perturbation. Although $\|E\|_{\infty}\geq\|E\|_{2}$,Since $\|E\|_{1}\|E\|_{\infty}\leq\|E\|_{2}^{2}$ (Stewart, 1990), the inequality follows from $\|E\|_{1}=\|E\|_{\infty}$ by symmetry. there are many settings where the two quantities are comparable; for example, if $E$ has a submatrix whose entries are identical and has zero entries otherwise, then $\|E\|_{\infty}=\|E\|_{2}$.

Theorem 2.2 provides the perturbation of individual eigenvectors, under a usual eigengap assumption. When $r$ and $\mu$ are not large, we incur an additional term $O(\|E\|_{2}/\delta\sqrt{d})$ in the bound. This is understandable, since $\|\widetilde{v}_{i}-v_{i}\|_{2}$ is typically $O(\|E\|_{2}/\delta)$.

We do not pursue the optimal bound in terms of $r$ and $\mu(V)$ in this paper, as the two quantities are not large in many applications, and the current proof is already complicated.

2 Rectangular matrices

Define $\mu_{0}=\max\{\mu(V),\mu(U)\}$, where $\mu(U)$ (resp. $\mu(V)$) is the coherence of $U$ (resp. $V$). This $\mu_{0}$ will appear in the statement of our results, as it controls both the structure of left and right singular spaces. When, specially, $A$ is a symmetric matrix, the spectral decomposition of $A$ is also the singular value decomposition (up to sign), and thus $\mu_{0}$ coincides with $\mu$ defined in Section 2.1.

Let $A_{r}=\sum_{i\leq r}\sigma_{i}u_{i}v_{i}^{T}$ be the best rank-$r$ approximation of $A$ under the Frobenius norm, and let $\varepsilon_{0}=\sqrt{d_{1}/d_{2}}\|A-A_{r}\|_{\infty}\vee\sqrt{d_{2}/d_{1}}\|A-A_{r}\|_{1}$, which also balances the two dimensions. Note that in the special case where $A$ is symmetric, this approximation error $\varepsilon_{0}$ is identical to $\varepsilon$ defined in Section 2.1. The next theorem bounds the perturbation of singular spaces.

Similar to Theorem 2.2, under an assumption of gaps between singular values, the next theorem bounds the perturbation of individual singular vectors.

Suppose the same assumption in Theorem 2.3. In addition, suppose $\delta_{0}$ satisfies $\delta_{0}>\|E\|_{2}$, and for any $i\in[r]$, the interval $[\sigma_{i}-\delta_{0},\sigma_{i}+\delta_{0}]$ does not contain any eigenvalues of $A$ other than $\sigma_{i}$. Then, up to sign,

As mentioned in the beginning of this section, we will use dilation to augment all $d_{1}\times d_{2}$ matrices into symmetric ones with size $d_{1}+d_{2}$. In order to balance the possibly different scales of $d_{1}$ and $d_{2}$, we consider a weighted max-norm. This idea will be further illustrated in Section 5.

3 Examples: which matrices have such structure?

In many problems, low-rank structure naturally arises due to the impact of pervasive latent factors that influence most observed data. Since observations are imperfect, the low-rank structure is often ‘perturbed’ by an additional sparse structure, gross errors, measurement noises, or the idiosyncratic components that can not be captured by the latent factors. We give some motivating examples with such structure.

Panel data in stock markets. Consider the excess returns from a stock market over a period of time. The driving factors in the market are reflected in the covariance matrix as a low rank component $A$. The residual covariance of the idiosyncratic components is often modeled by a sparse component $S$. Statistical analysis including PCA is usually conducted based on the estimated covariance matrix $\widetilde{A}=\widehat{\Sigma}$, which is perturbed from the true covariance $\Sigma=A+S$ by the estimation error $N$ (Stock and Watson, 2002; Fan et al., 2013). In Section 3.1, we will develop a robust estimation method in the presence of heavy-tailed return data.

Video surveillance. In image processing and computer vision, it is often desired to separate moving objects from static background before further modeling and analysis (Oliver et al., 2000; Hu et al., 2004). The static background corresponds to the low rank component $A$ in the data matrix, which is a collection of video frames, each consisting of many pixels represented as a long vector in the data matrix. Moving objects and noise correspond to the sparse matrix $S$ and noise matrix $N$. Since the background is global information and reflected by many pixels of a frame, it is natural for the incoherence condition to hold.

In our theorems, we require that the coherence $\mu$ is not too large. This is a natural structural condition associated with the low rank matrices. Consider the following very simple example: if the eigenvectors $v_{1},\ldots,v_{r}$ of the low rank matrix $A$ are uniform unit vectors in a sphere, then with high probability, $\max_{i}\|v_{i}\|_{\infty}=O(\sqrt{\log n})$, which implies $\mu=O(\log n)$. An intuitive way to understand the incoherence structure is that no coordinates of $v_{1}$ (or $v_{2},\ldots v_{r}$) are dominant. In other words, the eigenvectors are not concentrated on a few coordinates.

4 Other perturbation results

Although the eigenvector perturbation theory is well studied in numerical analysis, there is a renewed interest among statistics and machine learning communities recently, due to the wide applicability of PCA and other eigenvector-based methods. In Cai and Zhang (2016); Yu et al. (2015), they obtained variants or improvements of Davis-Kahan theorem (or Wedin’s theorem), which are user-friendly in the statistical contexts. These results assume the perturbation is deterministic, which is the same as Davis-Kahan theorem and Wedin’s theorem. In general, these results are sharp, even when the perturbation is random, as evidenced by the BBP transition (Baik et al., 2005).

However, these classical results can be suboptimal, when the perturbation is random and the smallest eigenvalue gap $\lambda_{1}-\lambda_{2}$ does not capture particular spectrum structure. For example, Vu (2011); O’Rourke et al. (2013) showed that with high probability, there are bounds sharper than the Wedin’s theorem, when the signal matrix is low-rank and satisfies certain eigenvalue conditions.

In this paper, our perturbation results are deterministic, thus the bound can be suboptimal when the perturbation is random with certain structure (e.g. the difference between sample covariance and population one for i.i.d. samples). However, the advantage of a deterministic result is that it is applicable to any random perturbation. This is especially useful when we cannot make strong random assumptions on the perturbation (e.g., the perturbation is an unknown sparse matrix). In Section 3, we will see examples of this type.

Application to robust covariance estimation

We will study the problem of robust estimation of covariance matrices and show the strength of our perturbation result. Throughout this section, we assume both rank $r$ and the coherence $\mu(V)$ are bounded by a constant, though this assumption can be relaxed. We will use $C$ to represent a generic constant, and its value may change from line to line.

To initiate our discussions, we first consider sub-Gaussian random variables. Let $X=(X_{1},\dots,X_{d})$ be a random $d$-dimensional vector with mean zero and covariance matrix

To apply Theorem 2.2, we treat $\Sigma_{1}$ as $A$ and $\widehat{\Sigma}-\Sigma_{1}$ as $E$. If the conditions in Theorem 2.2 are satisfied, we will obtain

Note there are simple bounds on $\|E\|_{\infty}$ and $\|E\|_{2}$:

By assuming a strong uniform eigengap, the conditions in Theorem 2.2 are satisfied, and the bound in (13) can be simplified. Define the uniform eigengap as

Note that $\gamma\leq\min\{\lambda_{r},\delta\}$, so if \gamma>C(1+\big{(}d\sigma^{2}+\lambda_{1}\big{)}\sqrt{\log d/n}), we have

In particular, when $\lambda_{1}\asymp\gamma$ and $\gamma\gg\max\{1,\sigma^{2}d\sqrt{\log d/n}\}$, we have

The above analysis pertains to the structure of sample covariance matrix. In the following subsections, we will estimate the covariance matrix using more complicated robust procedure. Our perturbation theorems in Section 2 provide a fast and clean approach to obtain new results.

2 PCA for robust covariance estimation

The usefulness of Theorem 2.2 is more pronounced when the random variables are heavy-tailed. Consider again the covariance matrix $\Sigma$ with structure (11). Instead of assuming sub-Gaussian distribution, we assume there exists a constant $C>0$ such that $\max_{j\leq d}EX_{j}^{4}<C$, i.e. the fourth moments of the random variables are uniformly bounded.

Unlike sub-Gaussian variables, there is no concentration bound similar to (12) for the empirical covariance matrix. Fortunately, thanks to recent advances in robust statistics (e.g., Catoni (2012)), robust estimate of $\Sigma$ with guaranteed concentration property becomes possible. We shall use the method proposed in Fan et al. (2017). Motivated by the classical $M$-estimator of Huber (1964), Fan et al. (2017) proposed a robust estimator for each element of $\widehat{\Sigma}$, by solving a Huber loss based minimization problem

where $l_{\alpha}$ is the Huber loss defined as

The parameter $\alpha$ is suggested to be $\alpha=\sqrt{nv^{2}/\log(\epsilon^{-1})}$ for $\epsilon\in(0,1)$, where $v$ is assumed to satisfy $v\geq\max_{ij}\sqrt{\text{Var}(X_{i}X_{j})}$. If $\log(\epsilon^{-1})\leq n/8$, Fan et al. (2017) showed

From this result, the next proposition is immediate by taking $\epsilon=d^{-3}$.

Suppose that there is a constant $C$ with $\max_{j\leq d}EX_{j}^{4}<C$. Then with probability greater than $1-d^{-1}(1+d^{-1})$, the robust estimate of covariance matrix with $\alpha=\sqrt{3nv^{2}\log(d)}$ satisfies

where $v$ is a pre-determined parameter assumed to be no less than $\max_{ij}\sqrt{\text{Var}(X_{i}X_{j})}$.

3 Robust covariance estimation via factor models

In this subsection, we will apply Theorem 2.2 to robust large covariance matrix estimation for approximate factor models in econometrics. With this theorem, we are able to extend the data distribution in factor analysis beyond exponentially decayed distributions considered by Fan et al. (2013), to include heavy-tailed distributions.

Suppose the observation $y_{it}$, say, the excess return at day $t$ for stock $i$, admits a decomposition

where $\Sigma_{u}:=\text{Cov}(u_{t})$. To circumvent the identifiability issue common in latent variable models, here we also assume, without loss of generality, $\text{Cov}(f_{t})=I_{r}$ and that $B^{T}B$ is a diagonal matrix, since rotating $B$ will not affect the above decomposition (16).

We will need two major assumptions for our analysis: (1) the factors are pervasive in the sense of Definition 3.1, and (2) there is a constant $C>0$ such that $\|\Sigma_{u}^{-1}\|_{2},\|\Sigma_{u}\|_{2}\leq C$, which are standard assumptions in the factor model literature. The pervasive assumption is reasonable in financial applications, since the factors have impacts on a large fraction of the outcomes (Chamberlain and Rothschild, 1982; Bai, 2003). If the factor loadings $\{b_{i}\}_{i=1}^{d}$ are regarded as random realizations from a bounded random vector, the assumption holds (Fan et al., 2013).

In the factor model (15), the factors are called pervasive if there is a constant $C>0$ such that $\|B\|_{\max}\leq C$ and the eigenvalues of the $r$ by $r$ matrix $B^{T}B/d$ are distinct and bounded away from zero and infinity.

Let $\{\lambda_{i},v_{i}\}_{i=1}^{r}$ be the top $r$ eigenvalues and eigenvectors of $\Sigma$, and similarly, $\{\overline{\lambda}_{i},\overline{v}_{i}\}_{i=1}^{r}$ for $BB^{T}$. In the following proposition, we show that pervasiveness is naturally connected to the incoherence structure. This connects well between the econometrics and machine learning literatures and provide a good interpretation on the concept of the incoherence. Its proof can be found in the appendix.

Our goal is to obtain a good covariance matrix estimator by exploiting the structure (16). Our strategy is to use a generalization of the principal orthogonal complement thresholding (POET) method proposed in Fan et al. (2013). The generic POET procedure encompasses three steps:

Given three pilot estimators $\widehat{\Sigma},\widehat{\Lambda}=\mbox{diag}(\widehat{\lambda}_{1},\dots,\widehat{\lambda}_{r}),\widehat{V}=(\widehat{v}_{1},\dots,\widehat{v}_{r})$ respectively for true covariance $\Sigma$, leading eigenvalues $\Lambda=\mbox{diag}(\lambda_{1},\dots,\lambda_{r})$ and leading eigenvectors $V=(v_{1},\dots,v_{r})$, compute the principal orthogonal complement $\widehat{\Sigma}_{u}$:

Apply the correlation thresholding to $\widehat{\Sigma}_{u}$ to obtain thresholded estimate $\widehat{\Sigma}_{u}^{\top}$ defined as follows:

where $s_{ij}(\cdot)$ is the generalized shrinkage function (Antoniadis and Fan, 2001; Rothman et al., 2009) and $\tau_{ij}=\tau(\hat{\sigma}_{u,ii}\hat{\sigma}_{u,jj})^{1/2}$ is an entry-dependent threshold. $\tau$ will be determined later in Theorem 3.1. This step exploits the sparsity of $\Sigma_{u}$.

Construct the final estimator $\widehat{\Sigma}^{\top}=\widehat{V}\widehat{\Lambda}\widehat{V}^{T}+\widehat{\Sigma}_{u}^{\top}$.

The key feature in the above procedure lies in the flexibility of choosing the pilot estimators in the first step. We will choose $\widehat{\Sigma}$ according to data generating distribution. Typically we can use $\hat{\lambda}_{i},\hat{v}_{i}$ for $i\leq r$ as the eigenvalues/vectors of $\widehat{\Sigma}$. However, $\hat{\Lambda}$ and $\hat{V}$ in general do not have to come from the spectral information of $\widehat{\Sigma}$ and can be obtained separately via different methods.

To guide the selection of proper pilot estimators, Fan et al. (2017+) provided a high level sufficient condition for this simple procedure to be effective, and its performance is gauged, in part, through the sparsity level of $\Sigma_{u}$, defined as $m_{d}:=\max_{i\leq d}\sum_{j\leq d}|\Sigma_{u,ij}|^{q}$. When $q=0$, $m_{d}$ corresponds to the maximum number of nonzero elements in each row of $\Sigma_{u}$. For completeness, we present the theorem given by Fan et al. (2017+) in the following.

Let $w_{n}=\sqrt{\log d/n}+1/\sqrt{d}$. Suppose there exists $C>0$ such that $\|\Sigma_{u}^{-1}\|,\|\Sigma_{u}\|\leq C$ and we have pilot estimators $\widehat{\Sigma},\widehat{\Lambda},\widehat{V}$ satisfying

Under the pervasiveness condition of the factor model (15), with $\tau\asymp w_{n}$, if $m_{d}w_{n}^{1-q}=o(1)$, the following rates of convergence hold with the generic POET procedure:

where $\|A\|_{\Sigma}=d^{-1/2}\|\Sigma^{-1/2}A\Sigma^{-1/2}\|_{F}$ is the relative Frobenius norm.

Let us decompose $\widehat{\Sigma}$ into a form such that Theorem 2.2 can be invoked:

where $\widehat{\Sigma}$ is viewed as $\widetilde{A}$, the low-rank part $\sum_{i=1}^{r}\overline{\lambda}_{i}\overline{v}_{i}\overline{v}_{i}^{T}$, which is also $BB^{T}$, is viewed as $A$, and the remaining terms are treated as $E$. The following results follow immediately.

Assume that there is a constant $C>0$ such that $\|\Sigma_{u}\|\leq C$. If the factors are pervasive, then with probability greater than $1-d^{-1}$, we have (19) – (21) hold with $\widehat{\lambda}_{i},\widehat{v}_{i}$ as the leading eigenvalues/vectors of $\widehat{\Sigma}$ for $i\leq r$. In addition, (22) and (23) hold.

The inequality (19) follows directly from Proposition 3.1 under the assumption of bounded fourth moments. It is also easily verifiable that (20), (21) follow from (19) by Weyl’s inequality and Theorem 2.2 (noting that $\|\Sigma_{u}\|_{\infty}\leq\sqrt{d}\|\Sigma_{u}\|$). See Section 3.2 for more details.

Note that in the case of sub-Gaussian variables, sample covariance matrix and its leading eigenvalues/vectors will also serve the same purpose due to (12) and Theorem 2.2 as discussed in Section 3.1.

Simulations

In this subsection, we implement numerical simulations to verify the perturbation bound in Theorem 2.2. We will show that the error behaves in the same way as indicated by our theoretical bound.

The perturbation of eigenvectors is measured by the element-wise error:

where $\{\widetilde{v}_{i}\}_{i=1}^{r}$ are the eigenvectors of $\widetilde{A}=A+E$ in the descending order.

To investigate how the error depends on $\gamma$ and $d$, we generate $E$ according to mechanism (a) with $s=10,L=3$, and run simulations in different parameter configurations: (1) let the matrix size $d$ range from $200$ to $2000$, and choose the eigengap $\gamma$ in $\{10,50,100,500\}$ (Figure 1); (2) fix the product $\gamma\sqrt{d}$ to be one of $\{2000,3000,4000,5000\}$, and let the matrix size $d$ run from $200$ to $2000$ (Figure 2).

To find how the errors behave for $E$ generated from different methods, we run simulations as in (1) but generate $E$ differently. We construct $E$ through mechanism (a) with $L=10,s=3$ and $L=0.6,s=50$, and also through mechanism (b) with $L^{\prime}=1.5,\rho=0.9$ and $L^{\prime}=7.5,\rho=0.5$ (Figure 3). The parameters are chosen such that $\|E\|_{\infty}$ is about $30$.

2 Simulation: robust covariance esitmation

We consider the performance of the generic POET procedure in robust covariance estimation in this subsection. Note that the procedure is flexible in employing any pilot estimators $\widehat{\Sigma},\widehat{\Lambda},\widehat{V}$ satisfying the conditions (19) – (21) respectively.

We implemented the robust procedure with four different initial trios: (1) the sample covariance $\widehat{\Sigma}^{S}$ with its leading $r$ eigenvalues and eigenvectors as $\widehat{\Lambda}^{S}$ and $\widehat{V}^{S}$; (2) the Huber’s robust estimator $\widehat{\Sigma}^{R}$ given in (14) and its top $r$ eigen-structure estimators $\widehat{\Lambda}^{R}$ and $\widehat{V}^{R}$; (3) the marginal Kendall’s tau estimator $\widehat{\Sigma}^{K}$ with its corresponding $\widehat{\Lambda}^{K}$ and $\widehat{V}^{K}$; (4) lastly, we use the spatial Kendall’s tau estimator to estimate the leading eigenvectors instead of the marginal Kendall’ tau, so $\widehat{V}^{K}$ in (3) is replaced with $\widetilde{V}^{K}$. We need to briefly review the two types of Kendall’s tau estimators here, and specifically give the formula for $\widehat{\Sigma}^{K}$ and $\widetilde{V}^{K}$.

Kendall’s tau correlation coefficient, for estimating pairwise comovement correlation, is defined as

Its population expectation is related to the Pearson correlation via the transform r_{jk}=\sin\Bigl{(}\frac{\pi}{2}\,E[\hat{\tau}_{jk}]\Bigr{)} for elliptical distributions (which are far too restrictive for high-dimensional applications). Then \hat{r}_{jk}=\sin\Bigl{(}\frac{\pi}{2}\hat{\tau}_{jk}\Bigr{)} is a valid estimation for the Pearson correlation $r_{jk}$. Letting $\widehat{R}=(\hat{r}_{jk})$ and $\widehat{D}=\mbox{diag}(\sqrt{\widehat{\Sigma}_{11}^{R}},\dots,\sqrt{\widehat{\Sigma}_{dd}^{R}})$ containing the robustly estimated standard deviations, we define the marginal Kendall’s tau estimator as

In the above construction of $\widehat{D}$, we still use the robust variance estimates from $\widehat{\Sigma}^{R}$.

The spatial Kendall’s tau estimator is a second-order U-statstic, defined as

In summary, Method (1) is designed for the case of sub-Gaussian data; Method (3) and (4) work under the situation of elliptical distribution; while Method (2) is proposed in this paper for the general heavy-tailed case with bounded fourth moments without further distributional shape constraints.

We simulated $n$ samples of $(f_{t}^{T},u_{t}^{T})^{T}$ from two settings: (a) a multivariate t-distribution with covariance matrix diag$\{I_{r},5I_{d}\}$ and various degrees of freedom ($\nu=3$ for very heavy tail, $\nu=5$ for medium heavy tail and $\nu=\infty$ for Gaussian tail), which is one example of the elliptical distribution (Fang et al., 1990); (b) an element-wise iid one-dimensional t distribution with the same covariance matrix and degrees of freedom $\nu=3,5$ and $\infty$, which is a non-elliptical heavy-tailed distribution.

Each row of coefficient matrix $B$ is independently sampled from a standard normal distribution, so that with high probability, the pervasiveness condition holds with $\|B\|_{\max}=O(\sqrt{\log d})$. The data is then generated by $y_{t}=Bf_{t}+u_{t}$ and the true population covariance matrix is $\Sigma=BB^{T}+5I_{d}$.

For $d$ running from $200$ to $900$ and $n=d/2$, we calculated errors of the four robust estimators in different norms. The tuning for $\alpha$ in minimization (14) is discussed more throughly in Fan et al. (2017). For the thresholding parameter, we used $\tau=2\sqrt{\log d/n}$. The estimation errors are gauged in the following norms: $\|\widehat{\Sigma}_{u}^{\top}-\Sigma_{u}\|$, $\|(\widehat{\Sigma}^{\top})^{-1}-{\Sigma}^{-1}\|$ and $\|\widehat{\Sigma}^{\top}-\Sigma\|_{\Sigma}$ as shown in Theorem 3.1. The two different settings are separately plotted in Figures 4 and 5. The estimation errors of applying sample covariance matrix $\widehat{\Sigma}^{S}$ in Method (1) are used as the baseline for comparison. For example, if relative Frobenius norm is used to measure performance, $\|(\widehat{\Sigma}^{\top})^{(k)}-{\Sigma}\|_{\Sigma}/\|(\widehat{\Sigma}^{\top})^{(1)}-{\Sigma}\|_{\Sigma}$ will be depicted for $k=2,3,4$, where $(\widehat{\Sigma}^{\top})^{(k)}$ are generic POET estimators based on Method ($k$). Therefore if the ratio curve moves below $1$, the method is better than naive sample estimator (Fan et al., 2013) and vice versa. The more it gets below $1$, the more robust the procedure is against heavy-tailed randomness.

The first setting (Figure 4) represents a heavy-tailed elliptical distribution, where we expect Methods (2), (3), (4) all outperform the POET estimator based on the sample covariance, i.e. Method (1), especially in the presence of extremely heavy tails (solid lines for $\nu=3$). As expected, all three curves under various measures show error ratios visibly smaller than $1$. On the other hand, if data are indeed Gaussian (dotted line for $\nu=\infty$), Method (1) has better behavior under most measures (error ratios are greater than $1$). Nevertheless, our robust Method (2) still performs comparably well with Method (1), whereas the median error ratios for the two Kendall’s tau methods are much worse. In addition, the IQR (interquartile range) plots reveal that Method (2) is indeed more stable than two Kendall’s tau Methods (3) and (4). It is also noteworthy that Method (4), which leverages the advantage of spatial Kendall’s tau, performs more robustly than Method (3), which solely base its estimation of the eigen-structure on marginal Kendall’s tau.

The second setting (Figure 5) provides an example of non-elliptical distributed data. We can see that the performance of the general robust Method (2) dominates the other three methods, which verifies the benefit of robust estimation for a general heavy-tailed distribution. Note that Kendall’s tau methods do not apply to distributions outside the elliptical family, excluding even the element-wise iid $t$ distribution in this setting. Nonetheless, even in the first setting where the data are indeed elliptical, with proper tuning, the proposed robust method can still outperform Kendall’s tau by a clear margin.

Proof Organization of Main Theorems

For shorthand, we write $\tau=\|E\|_{\infty}$, and $\kappa=\sqrt{d}\,\|EV\|_{\max}$. An obvious bound for $\kappa$ is $\kappa\leq\sqrt{r\mu}\,\tau$ (by Cauchy-Schwarz inequality). We will use these notations throughout this subsection.

Note that $E_{12}=E_{21}^{T}$ since $E$ is symmetric. Conceptually, the perturbation results in a rotation of $[V,V_{\bot}]$, and we write a candidate orthogonal basis as follows:

The approach of studying perturbation through a quadratic equation is known (see Stewart (1990) for example). Yet, to the best of our knowledge, existing results study perturbation under orthogonal-invariant norms (or unitary-invariant norms in the complex case), which includes a family of matrix operator norms and Frobenius norm, but excludes the matrix max-norm. The advantages of orthogonal-invariant norms are pronounced: such norms of a symmetric matrix only depend on its eigenvalues regardless of eigenvectors; moreover, with suitable normalization they are consistent in the sense $\|AB\|\leq\|A\|\cdot\|B\|$. See Stewart (1990) for a clear exposition.

The max-norm, however, does not possess these important properties. An imminent issue is that it is not clear how to relate $Q$ to $V_{\bot}Q$, which will appear in (29) after expanding $E$ according to (27), and which we want to control. Our approach here is to study $\overline{Q}:=V_{\bot}Q$ directly through a transformed quadratic equation, obtained by left multiplying $V_{\bot}$ to (29). Denote $H=V_{\bot}E_{21},\overline{Q}=V_{\bot}Q,\overline{L}_{1}=\Lambda_{1}+E_{11},\overline{L}_{2}=V_{\bot}(\Lambda_{2}+E_{22})V_{\bot}^{T}$. If we can find an appropriate matrix $\overline{Q}$ with $\overline{Q}=V_{\bot}Q$, and it satisfies the quadratic equation

then $Q$ also satisfies the quadratic equation (29). This is because multiplying both sides of (30) by $V_{\bot}^{T}$ yields (29), and thus any solution $\overline{Q}$ to (30) with the form $\overline{Q}=V_{\bot}Q$ must result in a solution $Q$ to (29).

Here, $\omega$ is defined as $\omega=8(1+r\mu)\kappa/(|\lambda_{r}|-\varepsilon)$.

The second claim of the lemma (i.e., the bound (31)) is relatively easy to prove once the first claim (i.e., the bound on $\|\overline{Q}\|_{\max}$) is proved. To understand this, note that we can rewrite $\overline{V}$ as $\overline{V}=(V+\overline{Q})(I_{r}+\overline{Q}^{T}\overline{Q})^{-1/2}$, and $\|\overline{Q}^{T}\overline{Q}\|_{\max}$ can be controlled by a trivial inequality $\|\overline{Q}^{T}\overline{Q}\|_{\max}\leq d\|\overline{Q}\|_{\max}^{2}\leq w^{2}$. To prove the first claim, we construct a sequence of matrices through recursion that converges to the fixed point $\overline{Q}$, which is a solution to the quadratic equation (30). For all iterates of matrices, we prove a uniform max-norm bound, which leads to a max-bound on $\|\overline{Q}\|_{\max}$ by continuity. To be specific, we initialize $\overline{Q}^{0}=0$, and given $\overline{Q}^{t}$, we solve a linear equation:

and the solution is defined as $\overline{Q}^{t+1}$. Under some conditions, the iterate $\overline{Q}^{t}$ converges to a limit $\overline{Q}$, which is a solution to (30). The next general lemma captures this idea. It follows from Stewart (1990) with minor adaptations.

Let $T$ be a bounded linear operator on a Banach space $\mathcal{B}$ equipped with a norm $\|\cdot\|$. Assume that $T$ has a bounded inverse, and define $\beta=\|T^{-1}\|^{-1}$. Let $\varphi:\mathcal{B}\to\mathcal{B}$ be a map that satisfies

for some $\eta\geq 0$. Suppose that $\mathcal{B}_{0}$ is a closed subspace of $\mathcal{B}$ such that $T^{-1}(\mathcal{B}_{0})\subseteq\mathcal{B}_{0}$ and $\varphi(\mathcal{B}_{0})\subseteq\mathcal{B}_{0}$. Suppose $y\in\mathcal{B}_{0}$ that satisfies $4\eta\|y\|<\beta^{2}$. Then, the sequence initialized with $x_{0}=0$ and iterated through

converges to a solution $x^{\star}$ to $Tx=y+\varphi(x)$. Moreover, we have $x^{\star}\subseteq\mathcal{B}_{0}$, and $\|x^{\star}\|\leq 2\|y\|/\beta$.

In other words, with a suitable orthogonal matrix $R$, the columns of $\overline{V}R$ are $\widetilde{v}_{1},\ldots,\widetilde{v}_{r}$.

It is easy to check that under the assumption of Theorem 2.1, the conditions required in Lemma 5.1 and Lemma 5.3 are satisfied. Hence, the two lemmas imply Theorem 2.1. ∎

We already provided a bound on $\|\overline{V}-V\|_{\max}$ in Lemma 5.1. By the triangular inequality, we can derive a bound on $\|\overline{V}\|_{\max}$. If we can prove a bound on $\|R-I_{r}\|_{\max}$, it will finally leads to a bound on $\|\widetilde{V}-V\|_{\max}$. In order to do so, we use the Davis-Kahan theorem to obtain an bound on $\langle\widetilde{v}_{i},v_{i}\rangle$ for all $i\in[r]$. This will lead to a max-norm bound on $R-I_{r}$ (with the price of potentially increasing the bound by a factor of $r$). The details about the proof of Theorem 2.2 are in the appendix.

We remark that, we assume conditions on $|\lambda_{r}|-\epsilon$ in Theorem 2.1 and Theorem 2.2, which are only useful in cases where $|\lambda_{r}|>\|A-A_{r}\|_{\infty}$. Ideally, we would like to have results with assumptions only involving $\lambda_{r}$ and $\lambda_{r+1}$, since Davis-Kahan theorem only requires a gap in neighboring eigenvalues. Unfortunately, unlike orthogonal-invariant norms that only depend on the eigenvalues of a matrix, the max-norm $\|\cdot\|_{\max}$ is not orthogonal-invariant, and thus it also depends on the eigenvectors of a matrix. For this reason, it is not clear whether we could obtain a lower bound on $\|T^{-1}\|_{\max}^{-1}$ using only the eigenvalues $\lambda_{r}$ and $\lambda_{r+1}$ so that we could apply Lemma 5.2. The analysis appears to be difficult if we do not have a bound on $\|T^{-1}\|_{\max}^{-1}$, considering that even in the analysis of linear equations, we also need invertibility, condition numbers, etc.

2 Asymmetric Case

Let $A^{d},E^{d}$ be $d_{1}+d_{2}$ square matrices defined as

Through the augmented matrices, we can transfer eigenvector results for symmetric matrices to singular vectors of asymmetric matrices. However, we cannot directly invoke the results proved for symmetric matrices, due to an issue about the coherence of $V^{d}$: when $d_{1}$ and $d_{2}$ are not comparable, the coherence $\mu(V^{d})$ can be very large even when $\mu(V)$ and $\mu(U)$ are bounded. To understand this, consider the case where $r=1$, $d_{1}\gg d_{2}$, and all entries of $U$ are $O(1/\sqrt{d_{1}})$, and all entries of $V$ are $O(1/\sqrt{d_{2}})$. Then, the coherences $\mu(U)$ and $\mu(V)$ are $O(1)$, but $\mu(V^{d})=O((d_{1}+d_{2})/d_{2})\gg 1$.

This unpleasant issue about the coherence, nevertheless, can be tackled if we consider a different matrix norm. In order to deal with the different scales of $d_{1}$ and $d_{2}$, we define the weighted max-norm for any matrix $M$ with $d_{1}+d_{2}$ rows as follows:

In other words, we rescale the top $d_{1}$ rows of $M$ by a factor of $\sqrt{d}_{1}$, and rescale the bottom $d_{2}$ rows by $\sqrt{d}_{2}$. This weighted norm serves to balance the potential different scales of $d_{1}$ and $d_{2}$.

The proofs of theorems in Section 2.2 will be almost the same with those in the symmetric case, with the major difference being the new matrix norm. Because the derivation is slightly repetitive, we will provide concise proofs in the appendix . Similar to the decomposition in (2.1),

where $A_{r}^{d}$ is has rank $2r$. Equivalently,

The next key lemma, which is parallel to Lemma 5.1, provides a bound on the solution $\overline{Q}^{d}$ to the quadratic equation

Here, $\omega_{0}$ is defined as $\omega_{0}=8(1+r\mu_{0})\kappa_{0}/3(\sigma_{r}-\varepsilon_{0})$.

In this lemma, the bound (38) bears a similar form to (31): if we consider the max-norm, the first $d_{1}$ rows of $\overline{V}^{d}-V^{d}$ correspond to the left singular vectors $u_{i}$’s, and they scale with $1/\sqrt{d_{1}}$; and the last $d_{2}$ rows correspond to the right singular vectors $v_{i}$’s, which scale with $1/\sqrt{d_{2}}$. Clearly, the weighted max-norm $\|\cdot\|_{w}$ indeed helps to balance the two dimensions.

Appendix A Proofs for Section 2.1

We have the following bound on $\|H\|_{\max}$:

By Cauchy-Schwarz inequality and the definition of $\mu$, for any $i,j\in[d]$,

Using the identity (40) and the above inequality, we derive

If $|\lambda_{r}|>\kappa r\sqrt{\mu}$, then $\overline{L}_{1}$ is an invertible matrix. Furthermore,

Let $Q_{0}$ be any $d\times r$ matrix with $\|Q_{0}\|_{\max}=1$. Note

We will derive upper bounds on $Q_{0}E_{11}$ and $\overline{L}_{2}Q_{0}$, and a lower bound on $Q_{0}\Lambda_{1}$. Since $E_{11}=V^{T}EV$ by definition, we expand $Q_{0}E_{11}$ and use a trivial inequality to derive

By Cauchy-Schwarz inequality and the definition of $\mu$ in (3), for $i,j\in[d]$,

Substituting $\|EV\|_{\max}=\kappa/\sqrt{d}$ into (42), we obtain an upper bound:

To bound $\overline{L}_{2}Q_{0}=(V_{\bot}E_{22}V_{\bot}^{T}+(A-A_{r}))Q_{0}$, we use the identity (40) and write

Using two trivial inequalities $\|EQ_{0}\|_{\max}\leq\|E\|_{\infty}\|Q_{0}\|_{\max}=\|E\|_{\infty}$ and $\|V^{T}Q_{0}\|_{\max}\leq\|V^{T}\|_{\infty}\|Q_{0}\|_{\max}\leq\sqrt{d}$, we have

In the proof of Lemma A.1, we showed $\|VV^{T}\|_{\max}\leq r\mu/d$. Thus,

Moreover, $\|(A-A_{r})Q_{0}\|_{\max}\leq\|A-A_{r}\|_{\infty}\|Q_{0}\|_{\max}=\varepsilon$. Combining the two bounds,

It is straightforward to obtain a lower bound on $\|Q_{0}\Lambda_{1}\|_{\max}$: since there is an entry of $Q_{0}$, say $(Q_{0})_{ij}$, that has an absolute value of $1$, we have

To show $\overline{L}_{1}$ is invertible, we use (42) and (45) to obtain

When $|\lambda_{r}|-\kappa r\sqrt{\mu}>0$, $\overline{L}_{1}$ must have full rank, because otherwise we can choose an appropriate $Q_{0}$ in the null space of $\overline{L}_{1}^{T}$ so that $Q_{0}\overline{L}_{1}=0$, which is a contradiction. To prove the second claim of the lemma, we combine the lower bound (45) with upper bounds (43) and (44) to derive

which is exactly the desired inequality. ∎

Next we prove Lemma 5.2. This lemma follows from Stewart (1990), with minor changes that involves $\mathcal{B}_{0}$. We provide a proof for the sake of completeness.

Let us write $\alpha=\|y\|$ for shorthand and recall $\beta=\|T^{-1}\|^{-1}$. As the first step, we show that the sequence $\{x_{k}\}_{k=0}^{\infty}$ is bounded. By construction in (34), we bound $\|x_{k+1}\|$ using $\|x_{k}\|$:

We use this inequality to derive an upper bound on $\{x_{k}\}$ for all $k$. We define $\xi_{0}=0$ and

then clearly $\|x_{k}\|\leq\xi_{k}$ (which can be shown by induction). It is easy to check (by induction) that the sequence $\{\xi_{k}\}_{k=1}^{\infty}$ is increasing. Moreover, since $4\alpha\eta<\beta^{2}$, the quadratic function

has two fixed points (namely solutions to $\phi(\xi)=\xi$), and the smaller one satisfies

If $\xi_{k}<\xi_{\star}$, then $\xi_{k+1}=\phi(\xi_{k})\leq\phi(\xi_{\star})=\xi_{\star}$. Thus, by induction, all $\xi_{k}$ are bounded by $\xi_{\star}$. This implies $\|x_{k}\|\leq\xi_{\star}<2\alpha/\beta$. The next step is to show that the sequence $\{x_{k}\}$ converges. Using the recursive definition (34) again, we derive

Since $4\alpha\eta/\beta^{2}<1$, the sequence $\{x_{k}\}_{k=0}^{\infty}$ is a Cauchy sequence, and convergence is secured. Let $x^{\star}\in\mathcal{B}$ be the limit. It is clear by assumption that $x_{k}\in\mathcal{B}_{0}$ implies $x_{k+1}\in\mathcal{B}_{0}$, so $x^{\star}\in\mathcal{B}_{0}$ and $\|x^{\star}\|\leq 2\alpha/\beta$ by continuity.

The final step is to show $x^{\star}$ is a solution to $Tx=y+\phi(x)$. Because $\{x_{k}\}_{k=0}^{\infty}$ is bounded and $\phi$ satisfies (33), the sequence $\{\phi(x_{k})\}_{k=0}^{\infty}$ converges to $\phi(x^{\star})$ by continuity and compactness. The linear operator $T$ is also continuous, so we can take limits on both sides of $Tx_{k+1}=y+\phi(x_{k})$, we conclude that $x^{\star}$ is a solution to $Tx=y+\phi(x)$. ∎

With all the preparations, we are now ready to present the key lemma. As discussed in Section 5, we set

Suppose $|\lambda_{r}|-\varepsilon>4r\mu(\tau+2r\kappa)$. Then there exists a solution $\overline{Q}\in\mathcal{B}_{0}$ to the equation (30) with

We will invoke Lemma 5.2 and apply it to the quadratic equation (30). To do so, we first check the conditions required in Lemma 5.2.

Let the linear operator $\mathcal{T}$ be $\mathcal{T}\overline{Q}=\overline{Q}\overline{L}_{1}-\overline{L}_{2}\overline{Q}$. By Lemma A.2, $\mathcal{T}$ has a bounded inverse, and $\beta:=\|\mathcal{T}^{-1}\|_{\max}^{-1}$ is bounded from below:

Let us define $\varphi$ by $\varphi(\overline{Q})=\overline{Q}H^{T}\overline{Q}$. To check the inequalities in (33), observe that

Thus, if we set $\eta=(1+r\mu)\kappa r\sqrt{d}$, then inequalities required in (33) are satisfied. For any $\overline{Q}$ with $\overline{Q}=V_{\bot}Q\in\mathcal{B}_{0}$, obviously $\varphi(\overline{Q})=V_{\bot}QH^{T}\overline{Q}\in\mathcal{B}_{0}$. To show $\mathcal{T}^{-1}(\overline{Q})\in\mathcal{B}_{0}$, let $Q_{0}=\mathcal{T}^{-1}(\overline{Q})$ and observe that

By definition, we know $\overline{L}_{2}Q_{0}=V_{\bot}(E_{22}+\Lambda_{2})V_{\bot}^{T}Q_{0}\in\mathcal{B}_{0}$, so we deduce $Q_{0}\overline{L}_{1}\in\mathcal{B}_{0}$. Our assumption implies $|\lambda_{r}|>\kappa r\sqrt{\mu}$, so by Lemma A.2, the matrix $\overline{L}_{1}$ is invertible, and thus $Q_{0}\in\mathcal{B}_{0}$. The last condition we check is $4\eta\|H\|_{\max}<\beta^{2}$. By Lemma A.1 and (46), this is true if

The above inequality holds when $|\lambda_{r}|>4r\mu(V)(\tau+2r\kappa)+\varepsilon$. Under this condition, we have, by Lemma 5.2,

where, the second inequality is due to $3r\mu(\tau+r\kappa)\leq 3(|\lambda_{r}|-\varepsilon)/4$. ∎

The next lemma is a consequence of Lemma A.3. We define, as in Lemma 5.1, that $\omega=8(1+r\mu)\kappa/(|\lambda_{r}|-\varepsilon)$.

Note that $\overline{\lambda}_{1}\leq\|\overline{Q}^{T}\overline{Q}\|\leq rd\|\overline{Q}\|_{\max}^{2}\leq r\omega^{2}$, which implies $\overline{\lambda}_{1}<1/2$. It is easy to check that $1+|x|\geq(1+x)^{-1/2}\geq 1-|x|$ whenever $|x|<1/2$. From this fact, we know $|(1+\overline{\lambda}_{i})^{-1/2}-1|\leq\overline{\lambda}_{i}\leq r\omega^{2}$. Using Cauchy-Schwarz inequality, we deduce that for any $j,k\in[d]$,

This leads to the desired max-norm bound. ∎

The first claim of the lemma (the existence of $\overline{Q}$ and its max-norm bound) follows directly from Lemma A.3. To prove the second claim, we split $\overline{V}-V$ into two parts:

where we used identity $V_{\bot}^{T}V_{\bot}=I_{d-r}$. Note that $r\omega<1/2$ implies $r\omega^{2}=(r\omega)^{2}/r<1/(4r)<1/2$. Thus, we can use Lemma A.4 and derive

where we used Cauchy-Schwarz inequality. Using the above inequality and the bound on $\|Q\|_{\max}$ (namely, the first claim in the lemma),

Simplifying the bound using $r\omega\leq 1/2$ and a trivial bound $\mu\geq 1$, we obtain (31). ∎

Using the identity in (39), it follows from Davis-Kahan $\sin\Theta$ theorem (Davis and Kahan, 1970) and Weyl’s inequality that

when $\delta_{r}>\|E\|_{2}$, where $\delta_{r}=|\lambda_{r}|-|\lambda_{r+1}|$. Since $\lambda_{r+1}\leq\|A-A_{r}\|_{2}\leq\varepsilon$ and $\|E\|_{2}\leq\tau$, the condition $|\lambda_{r}|-\varepsilon>3\tau$ implies $\delta_{r}>3\|E\|_{2}$. Hence, we have $d(\widetilde{V},V)<1/2$. Moreover,

This implies $d(\widetilde{V},\overline{V})\geq 1$, which is a contradiction. ∎