Spectral Methods for Data Science: A Statistical Perspective

Chapter 1 Introduction

In contemporary science and engineering applications, the volume of available data is growing at an enormous rate. The emergence of this trend is due to recent technological advances that have enabled the collection, transmission, storage and processing of data from every corner of our life, in the forms of images, videos, network traffic, email logs, electronic health records, genomic and genetic measurements, high-frequency financial trades, grocery transactions, online exchanges, and so on. In the meantime, modern applications often require reasonings about an unprecedented scale of features or parameters of interest. This gives rise to the pressing demand of developing low-complexity algorithms that can effectively distill actionable insights from large-scale and high-dimensional data. In addition to the curse of dimensionality, the challenge is further compounded when the data in hand are noisy, messy, and contain missing features.

Towards addressing the above challenges, spectral methods have emerged as a simple yet surprisingly effective approach to information extraction from massive and noisy data. In a nutshell, spectral methods refer to a collection of algorithms built upon the eigenvectors (resp. singular vectors) and eigenvalues (resp. singular values) of some properly designed matrices generated from data. Remarkably, spectral methods lend themselves to a diverse array of applications in practice, including community detection in networks [newman2006finding, abbe2017community, rohe2011spectral, mcsherry2001spectral], angular synchronization in cryo-EM [singer2011three, singer2011angular], joint image alignment [chen2016projected], clustering [von2007tutorial, ng2002spectral], ranking [negahban2016rank, chen2015spectral, chen2017spectral], dimensionality reduction [belkin2003laplacian], low-rank matrix estimation [achlioptas2007fast, keshavan2010matrix], tensor estimation [montanari2016spectral, cai2019tensor], covariance and precision matrix estimation [fan2013large, fan2020robust], shape reconstruction [li2004fast], econometric and financial modeling [fan2021recent], among others. Motivated by their applicability to numerous real-world problems, this monograph seeks to offer a unified and comprehensive treatment towards establishing the theoretical underpinnings for spectral methods, particularly through a statistical lens.

At the heart of spectral methods is the idea that the eigenvectors or singular vectors of certain data matrices reveal crucial information pertaining to the targets of interest. We single out a few examples that epitomize this idea.

Clustering corresponds to the grouping of individuals based on their mutual similarities, which constitutes a fundamental task in unsupervised learning and spans numerous applications such as image segmentation (e.g., grouping pixels based on the objects they represent in an image) [browet2011community] and community detection (e.g., grouping users on the basis of their social circles) [fortunato2016community]. For concreteness, let us take a look at a simple scenario with $n$ individuals such that: (1) there exists a latent partitioning that divides all individuals into two groups, with the first $n/2$ individuals belonging to the first group and the rest belonging to the second group (without loss of generality); and (2) we observe pairwise similarity measurements generated based on their group memberships. Ideally, if we know whether any two individuals belong to the same group or not, then we can form an adjacency matrix $\bm{A}=[A_{i,j}]_{1\leq i,j\leq n}$ such that

As a key observation, this matrix $\bm{A}$, as illustrated in Figure 1.1(a), turns out to be a rank-2 matrix

where $\bm{1}_{n}$ represents an $n$-dimensional all-one vector. After subtracting $\frac{1}{2}\bm{1}_{n}\bm{1}^{\top}_{n}$ from $\bm{A}$, the eigenvector \bm{u}_{2}\coloneqq[\begin{array}[]{cc}\bm{1}^{\top}_{n/2}&-\bm{1}^{\top}_{n/2}\end{array}] of the remaining component uncovers the underlying group structure; namely, all positive entries of $\bm{u}_{2}$ represent one group, with all negative entries of $\bm{u}_{2}$ reflecting another group. In reality, however, we typically only get to collect imprecise information about whether two individuals belong to the same group, thus resulting in a corrupted version of $\bm{A}$ (see Figure 1.1(b)). Fortunately, the eigenvector (the one corresponding to $\bm{u}_{2}$ above) of the observed data matrix (with proper arrangement) might continue to be informative, as long as the noise level is not overly high. To illustrate the practical applicability, we plot in Figure 1.1(c) the numerical performance of this approach, which allows for perfect clustering of all individuals for a wide range of noisy scenarios. Similar ideas continue to fare well on the clustering of real data, where we illustrate in Figure 1.2 that the penultimate eigenvector of a Laplacian matrix (also known as the Fiedler vector) of an undirected social network reveals two communities of 62 dolphins residing in Doubtful Sound, New Zealand.

If all sample vectors approximately lie within $\bm{U}^{\star}$, then one might be able to infer $\bm{U}^{\star}$ by inspecting the rank-$r$ leading eigenspace of $\bm{M}$ (or its variants), provided that the signal-to-noise ratio exceeds some reasonable level. This reflects the role of spectral methods in enabling meaningful dimensionality reduction and factor analysis.

In practice, a key benefit of PCA is its ability to remove nuance factors in, and extract out salient features from, each data point. As an illustration, the first four images of Figure 1.3 are representative ones sampled from a face dataset [georghiades2001few], which correspond to faces of the same person under different illumination and occlusion conditions. In contrast, the “eigenface” [turk1991face] depicted in the last image of Figure 1.3 corresponds to the first principal component (i.e., $r=1$), which effectively removes the nuance factors and highlights the feature of the face.

A proliferation of big-data applications has to deal with matrix estimation in the presence of missing data, either due to the infeasibility to acquire complete observations of a massive data matrix [davenport2016overview] such as the Netflix problem in recommender systems (as users only watch and rate a small fraction of movies), or because of the incentive to accelerate computation by means of sub-sampling [mahoney2016lecture]. Imagine that we are asked to estimate a large matrix $\bm{M}^{\star}=[M_{i,j}^{\star}]_{1\leq i,j\leq n}$, even though a dominant fraction of its entries are unseen. While in general we cannot predict anything about the missing entries, reliable estimation might become possible if $\bm{M}^{\star}$ is known a priori to enjoy a low-rank structure, as is the case in many applications like structure from motion [tomasi1992shape] and sensor network localization [javanmard2013localization]. This low-rank assumption motivates the use of spectral methods. More specifically, suppose the entries of $\bm{M}^{\star}$ are randomly sampled such that each entry is observed independently with probability $p\in(0,1]$. An unbiased estimate $\bm{M}=[M_{i,j}]_{1\leq i,j\leq n}$ of $\bm{M}^{\star}$ can be readily obtained via rescaling and zero filling (also called the inverse probability weighting method):

To capture the assumed low-rank structure of $\bm{M}^{\star}$, it is natural to resort to the best rank-$r$ approximation of $\bm{M}$ (with $r$ the true rank of $\bm{M}^{\star}$), computable through the rank-$r$ singular value decomposition of $\bm{M}$. Given its (trivial) success when $p=1$, we expect the algorithm to perform well when $p$ is close to 1. The key question, however, is where the algorithm stands if the vast majority of the entries is missing. While we shall illuminate this in Chapters 3 and 4, Figure 1.4 provides some immediate numerical assessment, which demonstrates the appealing performance of spectral methods—in terms of both Euclidean and entrywise estimation errors—even when the missing rate is quite high.

Another important application of spectral methods arises from the context of ranking, a task of central importance in, say, web search and recommendation systems. In a variety of scenarios, humans find it difficult to simultaneously rank many items, but relatively easier to express pairwise preferences. This gives rise to the problem of ranking based on pairwise comparisons. More specifically, imagine we are given a collection of $n$ items, and wish to identify top-ranked items based on pairwise preferences (with uncertainties in comparison outcomes) between observed pairs of items. A classical statistical model proposed by bradley1952rank, luce2012individual postulates the existence of a set of latent positive scores $\{w_{i}^{\star}\}_{1\leq i\leq n}$—each associated with an item—that determines the ranks of these items. The outcome of the comparison between items $i$ and $j$ is generated in a way that

As it turns out, the preference scores are closely related to the stationary distribution of a Markov chain associated with the above probability kernel, thus forming the basis of spectral ranking algorithms. To elucidate it in a little more detail, let us construct a probability transition matrix $\bm{P}^{\star}=[P_{i,j}^{\star}]_{1\leq i,j\leq n}$ with

Clearly, it forms a probability transition matrix since each element is nonnegative and the entries in each row add up to one. It is straightforward to verify that the score vector $\bm{w}^{\star}\coloneqq[w_{i}^{\star}]_{1\leq i\leq n}$ satisfies $\bm{w}^{\star\top}=\bm{w}^{\star\top}\bm{P}^{\star}$, namely $\bm{w}^{\star}$ is a left eigenvector of $\bm{P}^{\star}$ associated with eigenvalue one. A candidate method then consists of (i) forming an unbiased estimate of $\bm{P}^{\star}$ (which can be easily obtained using pairwise comparison outcomes), (ii) computing its left eigenvector (in fact, the leading left eigenvector), and (iii) reporting the ranking result in accordance with the order of the elements in this eigenvector. This spectral ranking scheme, which shares similar spirit with the celebrated PageRank algorithm [page1999pagerank], exhibits intriguing performance when identifying the top-ranked items, as showcased in the numerical experiments in Figure 1.5(b).

In all preceding applications, the core ideas underlying the development of spectral methods can be described in a unified fashion:

Identify a key matrix $\bm{M}^{\star}$—which is typically unobserved—whose eigenvectors or singular vectors disclose the information being sought after;

Construct a surrogate matrix $\bm{M}$ of $\bm{M}^{\star}$ using the data samples in hand, and compute the corresponding eigenvectors or singular vectors of this surrogate matrix.

Viewed in this light, this monograph aims to identify key factors—e.g., certain spectral structure of $\bm{M}^{\star}$ as well as the size of the approximation error $\bm{M}-\bm{M}^{\star}$—that exert main influences on the efficacy of the resultant spectral methods.

2 A modern statistical perspective

The idea of spectral methods can be traced back to early statistical literature on methods of moments (e.g., pearson1894contributions, hansen1982large), where one seeks to extract key parameters of the probability distributions of interest by examining the empirical moments of data. While classical matrix perturbation theory lays a sensible foundations for the analysis of spectral methods [stewart1990matrix], the theoretical understanding can be considerably enhanced through the lens of statistical modeling—a way of thinking that has flourished in the past decade. To the best of our knowledge, however, a systematic and comprehensive introduction to the modern statistical foundation of spectral methods, as well as an overview of recent advances, is previously unavailable.

The current monograph aims to fill this gap by developing a coherent and accessible treatment of spectral methods from a modern statistical perspective. Highlighting algorithmic implications that inform practice, our exposition gravitates around the following central questions: how to characterize the sample efficiency of spectral methods in reaching a prescribed accuracy level, and how to assess the stability of spectral methods in the face of random noise, missing data, and adversarial corruptions? We underscore several distinguishing features of our treatment compared to prior studies:

In comparison to the worst-case performance guarantees derived solely based on classical matrix perturbation theory, our statistical treatment emphasizes the benefit of harnessing the “typical” behavior of data models, which offers key insights into how to harvest performance gains by leveraging intrinsic properties of data generating mechanisms.

In contrast to classical asymptotic theory [van2000asymptotic], we adopt a non-asymptotic (or finite-sample) analysis framework that draws on tools from recent developments of concentration inequalities [tropp2015introduction] and high-dimensional statistics [wainwright2019high]. This framework accommodates the scenario where both the sample size and the number of features are enormous, and unveils a clearer and more complete picture about the interplay and trade-off between salient model parameters.

3 Organization

We now present a high-level overview of the structure of this monograph.

Chapter 5 concludes this monograph by identifying a few directions that are worthy of future investigation.

While this monograph pursues a coherent and accessible treatment that might appeal to a broad audience, it does not necessarily deliver the sharpest possible results for the applications discussed herein in terms of the logarithmic terms and/or pre-constants. The bibliographic notes at the end of each chapter contain information about the state-of-the-art theory for each application as a pointer to further readings.

4 What is not here and complementary readings

The topics presented in this monograph do not cover the tensor decomposition methods studied in another recent strand of work [anandkumar2014tensor]. While such tensor-based methods are also sometimes referred to as spectral methods, their primary focus is to invoke tensor decomposition to learn latent variables, based on higher-order moments estimated from data samples. We elect not to discuss this class of methods but instead refer the interested reader to the recently published monograph by MAL-057. Another monograph by kannan2009spectral provides an in-depth computational and algorithmic treatment of spectral methods from the perspective of theoretical computer science. The applications and results covered therein (e.g., fast matrix multiplication) complement the ones presented in the current monograph. In addition, spectral methods have been frequently employed to initialize nonconvex optimization algorithms. We will not elaborate on the nonconvex optimization aspect here but instead recommend the reader to the recent overview article by chi2019nonconvex. Finally, spectral methods are widely adopted to estimate high-dimensional covariance and precision matrices, and extract latent factors for econometric and statistical modeling. This topic alone has a huge literature, and we refer the interested reader to fan2020statistical for in-depth discussions.

5 Notation

Before moving forward, let us introduce some notation that will be used throughout this monograph.

When it comes to diagonal matrices, we employ $\mathsf{diag}([\theta_{1},\theta_{2},\cdots,\theta_{r}])$ to abbreviate the diagonal matrix with diagonal elements $\theta_{1},\cdots,\theta_{r}$. For any diagonal matrix $\bm{\Theta}=\mathsf{diag}([\theta_{1},\theta_{2},\cdots,\theta_{r}])$, we adopt the shorthand notation $\sin\bm{\Theta}\coloneqq\mathsf{diag}([\sin\theta_{1},\sin\theta_{2},\cdots,\sin\theta_{r}])$; the notation $\sin^{2}\bm{\Theta}$, $\cos\bm{\Theta}$, and $\cos^{2}\bm{\Theta}$ is defined analogously.

Finally, this monograph makes heavy use of the following standard notation: (1) $f(n)=O\left(g(n)\right)$ or $f(n)\lesssim g(n)$ means that there exists a universal constant $c>0$ such that $\left|f(n)\right|\leq c|g(n)|$ holds for all sufficiently large $n$; (2) $f(n)\gtrsim g(n)$ means that there exists a universal constant $c>0$ such that $|f(n)|\geq c\left|g(n)\right|$ holds for all sufficiently large $n$; (3) $f(n)\asymp g(n)$ means that there exist universal constants $c_{1},c_{2}>0$ such that $c_{1}|g(n)|\leq|f(n)|\leq c_{2}|g(n)|$ holds for all sufficiently large $n$; and (4) $f(n)=o(g(n))$ indicates that $f(n)/g(n)\rightarrow 0$ as $n\rightarrow\infty$. Additionally, we sometimes use $f(n)\gg g(n)$ (resp. $f(n)\ll g(n)$) to indicate that there exists some sufficiently large (resp. small) universal constant $c>0$ such that $|f(n)|\geq c\left|g(n)\right|$ (resp. $|f(n)|\leq c\left|g(n)\right|$).

Characterizing the performance of spectral methods requires understanding the perturbation of eigenspaces and/or that of singular subspaces. Classical matrix perturbation theory (e.g., stewart1990matrix) offers elementary toolkits that prove effective for this purpose, which we review in this chapter.

Setting the stage, consider a real-valued matrix $\bm{M}^{\star}$ and its perturbed version as follows

where $\bm{E}=\bm{M}-\bm{M}^{\star}$ denotes a real-valued perturbation or error matrix. In statistical applications, $\bm{M}$ can be an observed or estimated data matrix such as the sample covariance matrix, and $\bm{M}^{\star}$ is the target matrix such as the population covariance matrix. This chapter primarily aims to address the following questions by means of elementary linear algebra:

For a symmetric matrix $\bm{M}^{\star}$, how does the eigenspace change in response to a symmetric perturbation matrix $\bm{E}$?

For a general matrix $\bm{M}^{\star}$, how is the singular subspace affected as a result of the perturbation matrix $\bm{E}$?

We shall also explore eigenvector perturbation for a special class of asymmetric matrices: probability transition matrices.

We begin this chapter by gathering a few elementary materials in matrix analysis that prove useful for our theoretical development. The readers familiar with matrix analysis can proceed directly to Section 2.2.

This class of matrix norms enjoys several useful properties, as summarized in the following lemma. The proof can be found in stewart1990matrix.

For any unitarily invariant norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$, one has

Next, we review classical perturbation bounds for eigenvalues of symmetric matrices and for singular values of general matrices.

An immediate implication of Lemma 2.1.3 (resp. Lemma 2.1.5) is that the eigenvalues of a real symmetric matrix (resp. the singular values of a general matrix) are stable vis-à-vis small perturbations.

2 Preliminaries: Distance and angles between subspaces

In order to develop perturbation theory for eigenspaces and singular subspaces, we first need to delineate a metric that quantifies the proximity of two subspaces in a meaningful way.

For the sake of convenience, we further introduce two $n\times(n-r)$ matrices $\bm{U}^{\star}_{\perp}$ and $\bm{U}_{\perp}$, such that $[\bm{U}^{\star},\bm{U}_{\perp}^{\star}]$ and $[\bm{U},\bm{U}_{\perp}]$ are both $n\times n$ orthonormal matrices. In other words, $\bm{U}_{\perp}^{\star}$ and $\bm{U}_{\perp}$ represent the orthogonal complement of $\bm{U}^{\star}$ and $\bm{U}$, respectively.

2.2 Distance metrics and principal angles

To measure the distance between the two subspaces $\mathcal{U}$ and ${\mathcal{U}}^{\star}$, a naive idea is to employ the “metric” ${|\kern-1.07639pt|\kern-1.07639pt|\bm{U}-{\bm{U}}^{\star}|\kern-1.07639pt|\kern-1.07639pt|}$, where ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ is a certain norm of interest (e.g., the spectral norm or the Frobenius norm). An immediate drawback arises, however, since this “metric” does not take into account the global rotational ambiguity—namely, for any rotation matrix $\bm{R}\in\mathcal{O}^{r\times r}$, the columns of the matrix $\bm{U}\bm{R}$ also form a valid orthonormal basis of $\mathcal{U}$. This means that even when the two subspaces $\mathcal{U}$ and ${\mathcal{U}}^{\star}$ coincide, one might still have ${|\kern-1.07639pt|\kern-1.07639pt|\bm{U}-{\bm{U}^{\star}}|\kern-1.07639pt|\kern-1.07639pt|}\neq 0$, depending on how we rotate these matrices.

The takeaway of the above discussion is that any meaningful metric employed to measure the proximity of two subspaces should account for the rotational ambiguity properly. In what follows, we single out a few widely used metrics that meet such a requirement.

Distance with optimal rotation. Given the global rotational ambiguity, it is natural to first adjust the rotation matrix suitably before computing the distance. One choice is to measure the distance upon optimal rotation, namely,

where ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ is a certain norm to be chosen (e.g., the spectral norm or the Frobenius norm).

Distance between projection matrices. As an established fact, the projection matrix onto a subspace $\mathcal{U}$—given by $\bm{U}\bm{U}^{\top}$—is unique and unaffected by how $\bm{U}$ is rotated (since $\bm{U}\bm{U}^{\top}=\bm{U}\bm{R}\bm{R}^{\top}\bm{U}^{\top}$ for any rotation matrix $\bm{R}\in\mathcal{O}^{r\times r}$). The rotational invariance of the projection matrix motivates us to define the distance between $\mathcal{U}$ and ${\mathcal{U}}^{\star}$ as follows

Geometric construction via principal angles. Let $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{r}\geq 0$ be the singular values of $\bm{U}^{\top}{\bm{U}}^{\star}$, arranged in descending order. Given that $\|\bm{U}^{\top}{\bm{U}}^{\star}\|\leq\|\bm{U}\|\,\|{\bm{U}}^{\star}\|=1$, all the singular values $\{\sigma_{i}\}_{i=1}^{r}$ fall within the interval $$. Therefore, one can define the principal angles (or canonical angles) between the two subspaces of interest as

To see why this definition makes sense, consider the simplest example where $r=1$. In this case, the principal angle $\theta_{1}$ coincides with the conventionally defined angle between two unit vectors $\bm{U}$ and ${\bm{U}}^{\star}$. Armed with these angles, one might measure the distance between the subspaces $\mathcal{U}$ and ${\mathcal{U}}^{\star}$ through the following metric

where ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ is again some matrix norm to be selected, and

With slight abuse of notation, we can define other diagonal matrices such as $\cos\bm{\Theta}$ analogously, where $\cos(\cdot)$ is applied in an entrywise manner to the diagonal elements of $\bm{\Theta}$. Such matrices will be useful for future discussions.

2.3 Intimate connections between the distance metrics

It turns out that the metrics (2.3), (2.4) and (2.7) introduced above are tightly related, as we shall explain in this subsection. The proofs of all the results in this subsection are deferred to Section 2.6.

To begin with, we take a look at the relation between $\mathsf{dist}_{\mathsf{p},{\left|\kern-0.75346pt\left|\kern-0.75346pt\left|\cdot\right|\kern-0.75346pt\right|\kern-0.75346pt\right|}}(\cdot,\cdot)$ and $\mathsf{dist}_{\mathsf{sin},{\left|\kern-0.75346pt\left|\kern-0.75346pt\left|\cdot\right|\kern-0.75346pt\right|\kern-0.75346pt\right|}}(\cdot,\cdot)$, which is perhaps best illuminated by the following lemma.

Consider the settings of Section 2.2.1. If $2r\leq n$, then the singular values of $\bm{U}\bm{U}^{\top}-{\bm{U}}^{\star}{\bm{U}}^{\star\top}$ (including zeros) are given by

In a nutshell, Lemma 2.2.1 establishes an explicit link between (a) the difference of the projection matrices and (b) the principal angles between the two subspaces of interest. This lemma and its analysis unveil the following crucial equivalence relation under two of our favorite norms—the spectral norm and the Frobenius norm; in light of this, we might refer to these metrics as the $\sin{\bm{\Theta}}$ distances from time to time.

Consider the settings of Section 2.2.1, and recall the definition of $\sin\bm{\Theta}$ in (2.8). For any $1\leq r\leq n$, one has

Next, we move on to demonstrate the (near) equivalence of $\mathsf{dist}_{{\left|\kern-0.75346pt\left|\kern-0.75346pt\left|\cdot\right|\kern-0.75346pt\right|\kern-0.75346pt\right|}}(\cdot,\cdot)$ and $\mathsf{dist}_{\mathsf{p},{\left|\kern-0.75346pt\left|\kern-0.75346pt\left|\cdot\right|\kern-0.75346pt\right|\kern-0.75346pt\right|}}(\cdot,\cdot)$ under the above-mentioned two norms.

In words, $\mathsf{dist}_{{\left|\kern-0.75346pt\left|\kern-0.75346pt\left|\cdot\right|\kern-0.75346pt\right|\kern-0.75346pt\right|}}(\cdot,\cdot)$ and $\mathsf{dist}_{\mathsf{p},{\left|\kern-0.75346pt\left|\kern-0.75346pt\left|\cdot\right|\kern-0.75346pt\right|\kern-0.75346pt\right|}}(\cdot,\cdot)$ are equivalent up to a factor of $\sqrt{2}$, when ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ is the spectral norm or the Frobenius norm.

2.4 The distance metrics of choice in this monograph

In conclusion, the following metrics, which are seemingly distinct at first glance, are (nearly) equivalent in measuring the distance between two subspaces $\bm{U}$ and $\bm{U}^{\star}$:

when ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ represents either the spectral norm or the Frobenius norm. Viewed in this light, we shall mainly concentrate on the following metrics throughout the rest of this monograph:

3 Perturbation theory for eigenspaces

Armed with the above metrics for subspace distances, we are in a position to identify key factors that affect the perturbation of eigenvectors and eigenspaces.

Let $\bm{M}^{\star}$ and $\bm{M}=\bm{M}^{\star}+\bm{E}$ be two $n\times n$ real symmetric matrices. We express the eigendecomposition of $\bm{M}^{\star}$ and $\bm{M}$ as follows

Here, $\{\lambda_{i}\}$ (resp. $\{\lambda_{i}^{\star}\}$) denote the eigenvalues of $\bm{M}$ (resp. $\bm{M}^{\star}$), and $\bm{u}_{i}$ (resp. $\bm{u}_{i}^{\star}$) stands for the eigenvector associated with the eigenvalue $\lambda_{i}$ (resp. $\lambda_{i}^{\star}$). Additionally, we take

The matrices $\bm{U}^{\star}$, $\bm{U}_{\perp}^{\star}$, $\bm{\Lambda}^{\star}$, and $\bm{\Lambda}_{\perp}^{\star}$ are defined analogously.

3.2 A warm-up example

In general, the eigenvector/eigenspace of a real symmetric matrix might change drastically even upon a small perturbation. To understand this, consider the following toy example borrowed from hsu2016notes:

where $0<\epsilon<1$ can be arbitrarily small. It is straightforward to check that the leading eigenvectors of $\bm{M}^{\star}$ and $\bm{M}$ are given respectively by

which are both quite large regardless of the size of $\epsilon$ or the size of the perturbation $\|\bm{E}\|$.

On closer inspection, this “pathological” behavior comes up due to the fact that perturbation size $\epsilon$ is comparable to the eigengap of $\bm{M}^{\star}$ (namely, $\lambda_{1}(\bm{M}^{\star})-\lambda_{2}(\bm{M}^{\star})=2\epsilon$). This hints at the important role played by the eigengap in influencing eigenspace perturbation.

3.3 The Davis-Kahan sin𝚯𝚯\bm{\Theta} theorem

At the core of classical eigenspace perturbation theory lies the landmark result of davis1970rotation, which delivers powerful eigenspace perturbation bounds in terms of the size of the perturbation matrix as well as the associated eigengap. Here and throughout, for any symmetric matrix $\bm{A}$, we denote by $\mathsf{eigenvalues}(\bm{A})$ the set of eigenvalues of $\bm{A}$.

Consider the settings in Section 2.3.1. Assume that

This conclusion remains valid if Assumption (2.24) is replaced by

In fact, Theorem 2.3.1 can be generalized to accommodate any unitarily invariant norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$, in the sense that

The proof of Theorem 2.3.1 and Remark 2.3.2 is quite elementary and can be found in Section 2.3.4.

Theorem 2.3.1 is commonly referred to as the Davis-Kahan sin$\bm{\Theta}$ theorem, given that it concerns the sin${\bm{\Theta}}$ distance between subspaces. Both bounds scale linearly with the perturbation size, and are inversely proportional to the eigengap $\Delta$. Informally, if we view $\|\bm{E}\|$ as the noise size and interpret the eigengap as the “signal strength” (which dictates how easy it is to distinguish the $r$ eigenvalues of interest from the remaining spectrum), then Theorem 2.3.1 asserts that the eigenspace perturbation degrades gracefully as the signal-to-noise-ratio decreases.

The careful reader might notice that Theorem 2.3.1 stays silent on the allowable size $\|\bm{E}\|$ of the perturbation. Note, however, that a restriction on $\|\bm{E}\|$ is somewhat hidden in Assumptions (2.24) and (2.26). When the eigenvalues in $\bm{\Lambda}^{\star}$ (resp. $\bm{\Lambda}$) and $\bm{\Lambda}^{\star}_{\perp}$ (resp. $\bm{\Lambda}_{\perp}$) are suitably ordered, it is oftentimes more convenient to work with the following corollary, which makes apparent the constraint on the size $\|\bm{E}\|$ with regard to the eigengap of $\bm{M}^{\star}$.

Consider the settings in Section 2.3.1. Suppose that $|\lambda_{1}^{\star}|\geq|\lambda_{2}^{\star}|\geq\cdots\geq|\lambda_{r}^{\star}|>|\lambda_{r+1}^{\star}|\geq\cdots\geq|\lambda_{n}^{\star}|$ and $|\lambda_{1}|\geq|\lambda_{2}|\geq\cdots\geq|\lambda_{n}|$ (i.e., the eigenvalues are sorted by their magnitudes). If $\|\bm{E}\|<(1-1/\sqrt{2})(|\lambda_{r}^{\star}|-|\lambda_{r+1}^{\star}|)$, then

The proof of Corollary 2.3.4 is also given in Section 2.3.4.

3.4 Proof of the Davis-Kahan sin𝚯𝚯\bm{\Theta} theorem

The proof proceeds by controlling the distance metric {\big{|}\kern-1.07639pt\big{|}\kern-1.07639pt\big{|}\bm{U}_{\perp}^{\top}\bm{U}^{\star}\big{|}\kern-1.07639pt\big{|}\kern-1.07639pt\big{|}}, where ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ denotes a unitarily invariant norm.

We start by proving the theorem under Assumption (2.24), and claim that it suffices to consider the case where

In fact, if this condition is violated, then one can employ a “centering” trick by enforcing global offset to $\bm{M}^{\star}$ and $\bm{M}$ as follows

It is straightforwardly seen that (a) $\bm{M}^{\star}_{\mathsf{c}}$ (resp. $\bm{M}_{\mathsf{c}}$) and $\bm{M}^{\star}$ (resp. $\bm{M}$) share the same eigenvectors; (b) the eigenvalues of $\bm{M}^{\star}_{\mathsf{c}}$ (resp. $\bm{M}_{\mathsf{c}}$) associated with $\bm{U}^{\star}$ (resp. $\bm{U}_{\perp}$) reside within $[-\gamma,\gamma]$ (resp. $(-\infty,-\gamma-\Delta]\cup[\gamma+\Delta,\infty)$), where $\gamma=\frac{\beta-\alpha}{2}\geq 0$. Consequently, this reduces to a scenario that resembles (2.29). In addition, we isolate two immediate consequences of Assumptions (2.24) and (2.29) that prove useful:

where we recall that $\sigma_{\min}(\bm{\Lambda}_{\perp})$ is the minimal singular value of $\bm{\Lambda}_{\perp}$.

Armed with the above spectral conditions, we are prepared to study $\bm{U}_{\perp}^{\top}\bm{U}^{\star}$. This is controlled through the following identity (obtained by the definition of eigenvectors):

Let ${\bm{R}}\coloneqq\left(\bm{M}-\bm{M}^{\star}\right)\bm{U}^{\star}=\bm{E}\bm{U}^{\star}$. The triangle inequality then tells us that

Here, the middle line follows from Lemma 2.1.2 in Section 2.1, whereas the last inequality arises from the properties (2.30). As a consequence,

Next, we turn to the scenario where Assumption (2.26) is in effect; it can be analyzed in a similar manner and hence we remark only on the difference. Assuming (2.29) holds without loss of generality, we have

Applying the triangle inequality to (2.31) in a different way yields

a conclusion that coincides with (2.32). The rest of the proof is the same as the one in the previous case.

Before concluding, we remark that {\big{|}\kern-1.07639pt\big{|}\kern-1.07639pt\big{|}\bm{U}_{\perp}^{\top}\bm{U}^{\star}\big{|}\kern-1.07639pt\big{|}\kern-1.07639pt\big{|}}={\big{|}\kern-1.07639pt\big{|}\kern-1.07639pt\big{|}\sin\bm{\Theta}\big{|}\kern-1.07639pt\big{|}\kern-1.07639pt\big{|}} holds for any unitarily invariant norm ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$; see li1998relative. This together with the above analysis leads to Remark 2.3.2.

We first examine the spectral ranges of $\bm{\Lambda}_{\perp}$ and $\bm{\Lambda}^{\star}$. Let $\lambda_{i}(\bm{M}^{\star})$ (resp. $\lambda_{i}(\bm{M})$) be the $i$-th largest eigenvalue of $\bm{M}^{\star}$ (resp. $\bm{M}$), sorted by their values (as opposed to their magnitudes). Then Weyl’s inequality (cf. Lemma 2.1.3 in Section 2.1) asserts that

Suppose that $\bm{M}^{\star}$ has $r_{1}$ positive (resp. $r_{2}=r-r_{1}$ negative) eigenvalues whose magnitudes exceed $|\lambda_{r+1}^{\star}|$. Then for any $i$ obeying $1\leq i\leq r_{1}$ or $i>n-r_{2}$, the triangle inequality gives

where the last inequality arises from our assumption on $\|\bm{E}\|$. On the contrary, if $r_{1}<i\leq n-r_{2}$, then one has

As a consequence, there are exactly $r$ (resp. $n-r$) eigenvalues of $\bm{M}$ whose magnitudes exceed (resp. lie below) $|\lambda_{r+1}^{\star}|+\|\bm{E}\|$.

The above observation together with the ordering $|\lambda_{1}|\geq|\lambda_{2}|\geq\cdots\geq|\lambda_{n}|$ implies

In addition, the assumption that $|\lambda_{1}^{\star}|\geq|\lambda_{2}^{\star}|\geq\cdots\geq|\lambda_{n}^{\star}|$ tells us that

Taking $\beta=-\alpha=|\lambda_{r+1}^{\star}|+\|\bm{E}\|$ and $\Delta=|\lambda_{r}^{\star}|-|\lambda_{r+1}^{\star}|-\|\bm{E}\|>(|\lambda_{r}^{\star}|-|\lambda_{r+1}^{\star}|)/\sqrt{2}$, we can invoke Theorem 2.3.1 under Assumption (2.26) to establish the advertised results.

4 Perturbation theory for singular subspaces

There is no shortage of scenarios where the data matrices under consideration are asymmetric or rectangular. In these cases, one is often asked to study singular value decomposition (SVD) rather than eigendecomposition. Fortunately, the eigenspace perturbation theory can be naturally extended to accommodate perturbation of singular subspaces.

Here, $\sigma_{1}\geq\cdots\geq\sigma_{n_{1}}$ (resp. $\sigma_{1}^{\star}\geq\cdots\geq\sigma_{n_{1}}^{\star}$) stand for the singular values of $\bm{M}$ (resp. $\bm{M}^{\star}$) arranged in descending order, $\bm{u}_{i}$ (resp. $\bm{u}_{i}^{\star}$) denotes the left singular vector associated with the singular value $\sigma_{i}$ (resp. $\sigma_{i}^{\star}$), and $\bm{v}_{i}$ (resp. $\bm{v}_{i}^{\star}$) represents the right singular vector associated with $\sigma_{i}$ (resp. $\sigma_{i}^{\star}$). In addition, we denote

The matrices $\bm{\Sigma}^{\star},\bm{\Sigma}_{\perp}^{\star},\bm{U}^{\star},\bm{U}_{\perp}^{\star},\bm{V}^{\star},\bm{V}_{\perp}^{\star}$ are defined analogously.

4.2 Wedin’s sin𝚯𝚯\bm{\Theta} theorem

wedin1972perturbation developed a perturbation bound for singular subspaces that parallels the Davis-Kahan sin$\bm{\Theta}$ theorem for eigenspaces. In what follows, we present a version that is convenient for subsequent discussions in this monograph.

Consider the settings in Section 2.4.1. If $\|\bm{E}\|<\sigma_{r}^{\star}-\sigma_{r+1}^{\star}$, then one has

This theorem simultaneously controls the perturbation of left and right singular subspaces. As a worthy note, both the interaction between $\bm{E}$ and $\bm{U}^{\star}$, and that between $\bm{E}$ and $\bm{V}^{\star}$, come into play in determining the perturbation bounds. In particular, if $\|\bm{E}\|<(1-1/\sqrt{2})(\sigma_{r}^{\star}-\sigma_{r+1}^{\star})$, then one can apply Lemma 2.1.2 in Section 2.1 to obtain

akin to the eigenspace perturbation bounds (2.28).

4.3 Proof of the Wedin sin𝚯𝚯\bm{\Theta} theorem

We now present a proof of the Wedin theorem. Similar to the proof of the Davis-Kahan theorem, we start by bounding {\big{|}\kern-1.07639pt\big{|}\kern-1.07639pt\big{|}\bm{U}_{\perp}^{\top}\bm{U}^{\star}\big{|}\kern-1.07639pt\big{|}\kern-1.07639pt\big{|}}, where ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\cdot\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$ stands for any unitarily invariant norm. To this end, it is seen that

Here, the first identity is valid as long as $\bm{\Sigma}^{\star}$ is invertible (which is guaranteed since $\sigma_{\min}(\bm{\Sigma}^{\star})=\sigma_{r}^{\star}>\sigma_{r+1}^{\star}+\|\bm{E}\|>0$ under our assumption), the second line follows from the identities $\bm{M}-\bm{E}=\bm{M}^{\star}$ and $\bm{M}^{\star}=\bm{U}^{\star}\bm{\Sigma}^{\star}\bm{V}^{\star\top}+\bm{U}_{\perp}^{\star}\bm{\Sigma}_{\perp}^{\star}\bm{V}_{\perp}^{\star\top}$, the third line holds since $\bm{M}=\bm{U}\bm{\Sigma}\bm{V}^{\top}+\bm{U}_{\perp}\bm{\Sigma}_{\perp}\bm{V}_{\perp}^{\top}$, whereas the last identity exploits the property

Applying the triangle inequality and Lemma 2.1.2 in Section 2.1 to the identity (2.47) yields

Here, the second line uses the properties $\|\bm{\Sigma}^{\star-1}\|=1/\sigma_{r}^{\star}$ and $\|\bm{\Sigma}_{\perp}\|=\sigma_{r+1}$, while the last inequality follows from Weyl’s inequality $\sigma_{r+1}\leq\sigma_{r+1}^{\star}+\|\bm{E}\|$ (cf. Lemma 2.1.5 in Section 2.1). Repeating the same argument yields

To finish up, combine the inequalities (2.48) and (2.49) to obtain

When $\|\bm{E}\|<\sigma_{r}^{\star}-\sigma_{r+1}^{\star}$, we can rearrange terms to arrive at

The proof is then completed by invoking Lemmas 2.2.2 and 2.2.3.

5 Eigenvector perturbation for probability transition matrices

Thus far, our eigenvector perturbation analysis has been constrained to the set of symmetric matrices. Note, however, that the utility of eigenvectors is by no means confined to symmetric matrices. In fact, eigenvector analysis plays a vital role in studying asymmetric matrices as well, most notably the family of probability transition matrices of Markov chains. This section explores how to extend eigenvector perturbation theory to accommodate an important class of probability transition matrices associated with reversible Markov chains.

In words, the distribution $\bm{\pi}$ is invariant with respect to $\bm{P}$. Clearly, $\bm{\pi}$ is the left eigenvector of $\bm{P}$ associated with eigenvalue $1$, with the corresponding right eigenvector given by $\bm{1}$. By the Gershgorin circle theorem (see, e.g., olver2006applied), the modulus of all eigenvalues must be bounded by the maximum of the row sum, which is 1. Given that $1$ is an eigenvalue of $\bm{P}$, the largest modulus of the eigenvalues of $\bm{P}$ is precisely 1, and therefore $\bm{\pi}$ is the leading left eigenvector of $\bm{P}$. In addition, a Markov chain is said to be reversible when the following detailed balance equations are satisfied:

where $\bm{\pi}=[\pi_{i}]_{1\leq i\leq n}$ is the stationary distribution obeying (2.50). It will be seen in the proof of Theorem 2.5.1 that all eigenvalues of such a matrix ${\bm{P}}$ are real. For readers who wish an introduction to the basics of Markov chains, we recommend the monograph by bremaud2013markov.

The leading left eigenvectors of $\bm{P}^{\star}$ and $\bm{P}$—or equivalently, the vectors representing their stationary distributions—are denoted by $\bm{\pi}^{\star}$ and $\bm{\pi}$, respectively. Here, we allow $\bm{E}$ to be fairly general, meaning that $\bm{P}$ does not necessarily represent a reversible Markov chain. The question is: how does the matrix $\bm{E}$ affect the perturbation $\bm{\pi}-\bm{\pi}^{\star}$ of the leading left eigenvector of interest?

5.2 Perturbation of the leading eigenvector

Now we are ready to present the perturbation bound for the leading left eigenvector of a probability transition matrix, a result originally developed in chen2017spectral.

Consider the settings in Section 2.5.1. Suppose that $\bm{P}^{\star}$ represents a reversible Markov chain, whose stationary distribution vector $\bm{\pi}^{\star}$ is strictly positive. Assume that

The similarity between Theorem 2.5.1 and Corollary 2.3.4 is noteworthy. Indeed, recalling that the largest eigenvalue of the probability transition matrix $\bm{P}^{\star}$ is precisely 1, one might view $1-\max\left\{\lambda_{2}(\bm{P}^{\star}),-\lambda_{n}(\bm{P}^{\star})\right\}$ as the gap between the first and the second largest eigenvalues of $\bm{P}^{\star}$ (in magnitude), akin to the eigengap $|\lambda_{r}^{\star}|-|\lambda_{r+1}^{\star}|$ in Corollary 2.3.4 (with $r=1$). In words, Theorem 2.5.1 guarantees that as long as the size of the perturbation matrix $\bm{E}$ is not too large, the perturbation of the leading left eigenvector—or equivalently, the perturbation of the stationary distribution of the associated Markov chain—is proportional to the size of the noise when projected onto the direction $\bm{\pi}^{\star}$, as measured by $\|\bm{\pi}^{\star\top}\bm{E}\|_{\bm{\pi}^{\star}}$. As we shall demonstrate in Section 3.6, this perturbation theory delivers powerful techniques for analyzing the ranking problem described previously in Chapter 1.

Sensitivity and perturbation analyses for the steady-state distributions of Markov chains have been studied in the literature; see, e.g., mitrophanov2005sensitivity, liu2012perturbation, jiang2017unified, rudolf2018perturbation and the references therein.

5.3 Proof of Theorem 2.5.1

Since $\bm{\pi}^{\star}$ and $\bm{\pi}$ denote respectively the leading left eigenvectors of $\bm{P}^{\star}$ and $\bm{P}$, namely,

the perturbation $\bm{\pi}-\bm{\pi}^{\star}$ admits the following decomposition

Here, the last relation hinges upon the fact that $\bm{\pi}$ and $\bm{\pi}^{\star}$ are probability vectors, and hence $\left(\bm{\pi}-\bm{\pi}^{\star}\right)^{\top}\bm{1}=1-1=0$. Apply the triangle inequality with respect to the norm $\|\cdot\|_{\bm{\pi}^{\star}}$ to obtain

where the last line relies on the definition of the matrix norm $\|\cdot\|_{\bm{\pi}^{\star}}$. Rearranging terms, we are left with

with the proviso that $\|\bm{P}-\bm{P}^{\star}\|_{\bm{\pi}^{\star}}+\|\bm{P}^{\star}-\bm{1}\bm{\pi}^{\star\top}\|_{\bm{\pi}^{\star}}<1$. The proof would then be completed as long as one could justify that

with $\|\cdot\|$ the usual spectral norm, where the last line replaces \big{(}\bm{\Pi}^{\star}\big{)}^{1/2}\bm{x} with $\bm{v}$. Consequently, we obtain

where we define $\bm{S}^{\star}\coloneqq\left(\bm{\Pi}^{\star}\right)^{1/2}\bm{P}^{\star}\left(\bm{\Pi}^{\star}\right)^{-1/2}$ and \bm{\pi}^{\star}_{1/2}\coloneqq\big{[}\sqrt{\pi_{i}^{\star}}\,\big{]}_{1\leq i\leq n}. Several basic properties regarding $\bm{S}^{\star}$ are in order; see bremaud2013markov.

Since $\bm{P}^{\star}$ represents a reversible Markov chain with stationary distribution $\bm{\pi}^{\star}$, the matrix $\bm{S}^{\star}$ is symmetric, whose eigenvalues are real-valued. This can be verified by the detailed balance equations (2.51).

Given that $\bm{S}^{\star}$ is obtained via a similarity transformation of $\bm{P}^{\star}$, we see that $\bm{S}^{\star}$ and $\bm{P}^{\star}$ share the same set of eigenvalues. This can easily be verified from the definition of eigenvectors:

In particular, $\lambda_{1}(\bm{S}^{\star})=\lambda_{1}(\bm{P}^{\star})=1$, and $\bm{\pi}^{\star}_{1/2}$ is precisely the eigenvector of $\bm{S}^{\star}$ associated with $\lambda_{1}(\bm{S}^{\star})=1$. Thus, from the eigendecomposition of the symmetric matrix $\bm{S}^{\star}$, it is easy to see that the eigenvalues of $\bm{S}^{\star}-\bm{\pi}^{\star}_{1/2}(\bm{\pi}^{\star}_{1/2})^{\top}$ are $0,\lambda_{2}(\bm{S}^{\star}),\cdots,\lambda_{n}(\bm{S}^{\star})$.

Taking the preceding facts collectively, we reach

Here, (i) relies on Property (c), while (ii) follows from Property (b). This concludes the proof.

6 Appendix: Proofs of auxiliary lemmas in Section 2.2

Given that singular values are unitarily invariant, it suffices to look at the singular values of the following matrix

Consequently, the singular values of $\bm{U}\bm{U}^{\top}-{\bm{U}}^{\star}{\bm{U}}^{\star\top}$ are composed of those of $\bm{U}^{\top}{\bm{U}}_{\perp}^{\star}$ and those of $\bm{U}_{\perp}^{\top}{\bm{U}}^{\star}$ combined. It then boils down to characterizing the spectrum of $\bm{U}^{\top}{\bm{U}}_{\perp}^{\star}$ and $\bm{U}_{\perp}^{\top}{\bm{U}}^{\star}$.

To pin down the singular values of $\bm{U}^{\top}{\bm{U}}_{\perp}^{\star}$, we first turn attention to the eigenvalues of $\bm{U}^{\top}\bm{U}_{\perp}^{\star}\bm{U}_{\perp}^{\star\top}\bm{U}$. Assuming that the SVD of $\bm{U}^{\top}\bm{U}^{\star}$ is given by $\bm{X}\bm{\Sigma}\bm{Y}^{\top}$ (where $\bm{X}$ and $\bm{Y}$ are $r\times r$ orthonormal matrices, and $\bm{\Sigma}$ is diagonal), we can derive

Here, the penultimate identity follows from our construction (cf. (2.5)), where we define $\cos\bm{\Theta}\coloneqq\mathsf{diag}([\cos\theta_{1},\cdots,\cos\theta_{r}])$. Therefore, for any $1\leq i\leq r$, the $i$-th largest singular value of $\bm{U}^{\top}\bm{U}_{\perp}^{\star}$ obeys

which results from the ordering in (2.6). This means that, if $r\leq n-r$, then the singular values of $\bm{U}^{\top}{\bm{U}}_{\perp}^{\star}$ are precisely given by $\{\sin\theta_{i}\}_{1\leq i\leq r}$. Repeating this argument reveals that the singular values of $\bm{U}_{\perp}^{\top}{\bm{U}}^{\star}$ are also $\{\sin\theta_{i}\}_{1\leq i\leq r}$ if $r\leq n-r$.

Combining the above observations thus completes the proof.

6.2 Proof of Lemma 2.2.2

A closer inspection of the proof of Lemma 2.2.1 (in particular, (2.60) and the orthonormality of $\bm{X}$) reveals that

where we have used the basic property $\mathsf{Tr}(\bm{A}\bm{B})=\mathsf{Tr}(\bm{B}\bm{A})$. Similarly,

Note that the above identities hold for all $1\leq r\leq n$. In addition, the relation (2.59) tells us that

Putting the above identities together immediately establishes the advertised results.

6.3 Proof of Lemma 2.2.3

Here, the penultimate line relies on the singular value decomposition $\bm{U}^{\top}\bm{U}^{\star}=\bm{X}\bm{\Sigma}\bm{Y}^{\top}$, while the two identities in the last line result from the orthonormality of $\bm{X}$ and $\bm{Y}$, respectively. In addition, note that

where the first inequality holds since $\bm{X}$ and $\bm{Y}$ are both orthonormal matrices and hence $\bm{X}\bm{Y}^{\top}$ is also orthonormal.

On the other hand, we make the observation that

where the last relation holds since $\bm{X}\bm{\Sigma}\bm{Y}^{\top}$ is the SVD of $\bm{U}^{\top}\bm{U}^{\star}$. Continue the derivation to obtain

Here, (i) follows by setting $\bm{Q}=\bm{R}^{\top}\bm{X}$ (since both $\bm{X}$ and $\bm{R}$ are orthonormal matrices), (ii) results from the unitary invariance of the spectral norm, whereas (iii) holds by setting $\bm{O}=\bm{Y}^{\top}\bm{Q}$. Moreover, recognizing that $\|\bm{\Sigma}\bm{O}\|\leq\|\bm{\Sigma}\|\cdot\|\bm{O}\|\leq 1$ (and hence $2\bm{I}_{r}-\bm{\Sigma}\bm{O}-\bm{O}^{\top}\bm{\Sigma}\succeq\bm{0}$), one can obtain

Here, the inequality follows by taking $\bm{u}$ to be $\bm{e}_{r}$ (recall that by construction, $\sigma_{r}=\cos\theta_{r}\geq 0$ is the smallest singular value of $\bm{\Sigma}$), and the penultimate line holds by combining the facts $|\bm{e}_{r}^{\top}\bm{O}\bm{e}_{r}|\leq\|\bm{O}\|=1$ and $\bm{e}_{r}^{\top}\bm{e}_{r}=1$. Putting (2.65) and (2.64) together yields

where we again use the inequality $2\sin(\theta/2)\geq\sin\theta$ for all $\theta\in[0,\pi/2]$.

Finally, invoking the relation $\|\sin\bm{\Theta}\|=\|\bm{U}\bm{U}^{\top}-\bm{U}^{\star}\bm{U}^{\star\top}\|$ (see Lemma 2.2.2) establishes the claimed spectral norm bounds.

Regarding the Frobenius norm upper bound, one sees that

where (i) holds since $\bm{U}$ and $\bm{U}^{\star}$ are both $n\times r$ matrices with orthonormal columns, and (ii) follows since $\bm{X}^{\top}\bm{X}=\bm{Y}^{\top}\bm{Y}=\bm{I}$ (and hence $\mathsf{Tr}(\bm{Y}\bm{X}^{\top}\bm{X}\bm{\Sigma}\bm{Y}^{\top})=\mathsf{Tr}(\bm{Y}^{\top}\bm{Y}\bm{X}^{\top}\bm{X}\bm{\Sigma})=\mathsf{Tr}(\bm{\Sigma})$). Furthermore,

where (iii) holds by construction (cf. (2.5)), and the last identity results from Lemma 2.2.2. This taken collectively with (2.66) reveals that

where the first inequality holds since $\bm{X}$ and $\bm{Y}$ are both orthonormal matrices and hence $\bm{X}\bm{Y}^{\top}$ is also orthonormal.

With regards to the Frobenius norm lower bound, it is seen that

Here, (iii) sets $\bm{Q}=\bm{X}^{\top}\bm{R}\bm{Y}$ and identifies $\bm{\Sigma}$ as $\cos\bm{\Theta}$, (iv) comes from the elementary inequality $\langle\bm{A},\bm{B}\rangle\leq\|\bm{A}\|\,\|\bm{B}\|_{*}$, whereas the last line follows since $\cos\theta_{i}\geq 0$. Additionally, it is easily seen that

where the penultimate relation follows from the elementary inequality $2\sin(\theta/2)\geq\sin\theta$ (which holds for any $0\leq\theta\leq\pi/2$), and the last line invokes Lemma 2.2.2. Combining the inequalities (2.68) and (2.69), we establish the claimed lower bound.

7 Notes

Matrix perturbation theory is a firmly established topic that has been extensively studied in the past several decades. Two classic books that offer in-depth discussions of perturbation theory for eigenspaces and singular subspaces are stewart1990matrix, sun1987perturbation. Other valuable resources on this topic include bhatia2013matrix, horn2012matrix. The exposition herein is largely influenced by the excellent lecture notes by montanari2011notes, hsu2016notes. In addition, the book [kato2013perturbation] offers a more abstract treatment of perturbation theory from the viewpoint of linear operators. Several variants of the sin$\bm{\Theta}$ theorem amenable to statistical analysis are available in the statistics literature as well (e.g., MR3371006, vu2013minimaxSPCA, cai2018rate, zhang2018heteroskedastic).

We point out several well-known extensions of the theorems presented in this chapter. To begin with, the current exposition restricts attention to the real case for simplicity, while in fact all results herein generalize to the complex-valued case [stewart1990matrix]. In addition, Theorem 2.4.1 together with Lemma 2.2.3 reveals the existence of two rotation matrices $\bm{R}_{U}$ and $\bm{R}_{V}$ obeying

but falls short of illuminating the connection between $\bm{R}_{U}$ and $\bm{R}_{V}$. An extension derived in dopico2000sintheta establishes a similar perturbation bound even when $\bm{R}_{U}$ and $\bm{R}_{V}$ are taken to be the same rotation matrix.

In order to invoke the sin$\bm{\Theta}$ theorems (Theorems 2.3.1 and 2.4.1), an important ingredient lies in developing a tight upper bound on the spectral norm $\|\bm{E}\|$ of the perturbation matrix $\bm{E}$. This is where statistical/probabilistic tools play a major role. Rather than presenting an encyclopedia of probabilistic techniques (which can be gleaned from tropp2015introduction, vershynin2016high, boucheron2013concentration, wainwright2019high, tropp2011freedman, raginsky2013concentration, howard2020time), this monograph singles out only two useful matrix concentration inequalities that suffice for the applications considered herein.

The first result is an extension of the celebrated matrix Bernstein inequality [oliveira2009concentration, tropp2012user, hopkins2016fast]. This is an elegant and convenient tail bound for the sum of independent random matrices, resulting in effective performance guarantees for a diverse array of statistical applications. We refer the interested reader to tropp2015introduction for a highly accessible introduction of the classical matrix Bernstein inequality, and hopkins2016fast for a proof of the truncated variant stated in Theorem 3.1.1.

Let $\{\bm{X}_{i}\}_{1\leq i\leq m}$ be a sequence of independent real random matrices with dimension $n_{1}\times n_{2}$. Suppose that for all $1\leq i\leq m$,

hold for some quantities $0\leq q_{0}\leq 1$ and $q_{1}\geq 0$. In addition, define the matrix variance statistic $v$ as

Note that when the $\bm{X}_{i}$’s are i.i.d. zero-mean random matrices, the matrix variance statistic simplifies to

To make it more user-friendly, we record a straightforward consequence of Theorem 3.1.1 as follows.

Suppose the assumptions of Theorem 3.1.1 hold, and set $n\coloneqq\max\{n_{1},n_{2}\}$. For any $a\geq 2$, with probability exceeding $1-2n^{-a+1}-mq_{0}$ one has

Let $\{\bm{X}_{i}\}_{1\leq i\leq m}$ be a set of independent real random matrices with dimension $n_{1}\times n_{2}$. Suppose that

Set $n\coloneqq\max\{n_{1},n_{2}\}$, and recall the definition of variance statistic in (3.2). For any $a\geq 2$, with probability exceeding $1-2n^{-a+1}$ one has

By virtue of the above inequalities, the key to bounding $\left\|\sum\nolimits_{i}\bm{X}_{i}\right\|$ largely lies in controlling the following two crucial quantities:

where the former depends on the number $m$ of random matrices involved.

An important family of random matrices that merits special attention comprises the ones with independent random entries, that is, matrices of the form $\bm{X}=[X_{i,j}]_{1\leq i,j\leq n}$ with independent $X_{i,j}$’s. While the spectral norm of such a matrix can also be analyzed via matrix Bernstein (by treating $\bm{X}$ as the sum of independent random matrices $X_{i,j}\bm{e}_{i}\bm{e}_{j}^{\top}$), this approach is typically loose in terms of the logarithmic factor. Motivated by the abundance of such random matrices in practice, we record below a strengthened non-asymptotic spectral norm bound, which is of significant utility and is tighter than what matrix Bernstein has to offer for this case.

Then there exists some universal constant $c>0$ such that for any $t\geq 0$,

This result, which appeared in bandeira2016sharp, can be established via tighter control of the expected spectral norm in conjunction with Talagrand’s concentration inequality. Two remarks are in order.

which enjoys the desired symmetry and can be analyzed directly using Theorem 3.1.5. The resulting bound on $\|\mathcal{S}(\bm{X})\|$ can be translated back to $\|\bm{X}\|$ via the elementary identity $\|\bm{X}\|=\|\mathcal{S}(\bm{X})\|$. For conciseness, we will occasionally apply Theorem 3.1.5 directly to asymmetric matrices without invoking the dilation trick.

with probability at least $1-n^{-8}$ for some constant $\widetilde{c}>0$. To see this, it suffices to set $\widetilde{c}=\sqrt{9c}$ and take $t=B\sqrt{9c\log n}$ in (3.8).

The inequality (3.9) continues to hold if we replace $n^{-8}$ with $n^{-\alpha}$ for any positive constant $\alpha>0$. Here and below, we often go with the artificial choice like $n^{-8}$ since it is small enough for our purpose.

2 Low-rank matrix denoising

To catch a glimpse of the effectiveness of the approach we have introduced, let us start by trying it out on a warm-up example: low-rank matrix denoising.

where $\bm{E}=[E_{i,j}]_{1\leq i,j\leq n}$ is a symmetric noise matrix. It is assumed that the entries $\{{E}_{i,j}\}_{i\geq j}$ are independently generated obeying

The aim is to estimate the eigenspace $\bm{U}^{\star}$ from the data matrix $\bm{M}$. Despite its simplicity, this problem has been extensively studied in the literature [koltchinskii2016perturbation, bao2018singular, ding2020high, xia2019normal, li2021minimax]. It also bears close relevance to the famous angular/phase synchronization problem [singer2011angular, bandeira2017tightness].

In order to estimate the low-rank factors specified by $\bm{U}^{\star}$, a natural scheme is to resort to the rank-$r$ leading eigenspace of the data matrix $\bm{M}$. More precisely, denote by $\lambda_{1},\cdots,\lambda_{n}$ the eigenvalues of $\bm{M}$ sorted by their magnitudes, i.e.,

2.2 Performance guarantees

We now examine the accuracy of the above spectral estimate. Towards this, a key step lies in bounding the spectral norm of the noise matrix $\bm{E}$. We claim for the moment that (which will be established in Section 3.2.3)

with probability at least $1-O(n^{-8})$. Armed with this claim and the fact $\lambda_{r+1}^{\star}=0$, we are in a situation where it is quite easy to see how the Davis-Kahan theorem applies. According to Corollary 2.3.4, with probability greater than $1-O(n^{-8})$ one has

provided that the noise variance is sufficiently small obeying $\sigma\sqrt{n}\leq\frac{1-1/\sqrt{2}}{5}|\lambda_{r}^{\star}|$ so that $\|\bm{E}\|\leq(1-1/\sqrt{2})|\lambda_{r}^{\star}|$.

The tightness of the statistical guarantee (3.13) can be assessed when compared with the minimax lower bound. For instance, it is well-known in the literature (e.g., cheng2020tackling) that: even for the case with $r=1$, one cannot hope to achieve \mathsf{dist}\big{(}\widehat{\bm{U}},\bm{U}^{\star}\big{)}=o(\sigma\sqrt{n}/|\lambda_{r}^{\star}|)—in a minimax sense—regardless of the estimator $\widehat{\bm{U}}$ in use. Consequently, the spectral method turns out to be orderwise statistically optimal for low-rank matrix denoising.

Before concluding, we record several immediate consequences of the above analysis that will be useful later on. Specifically, assuming that $\sigma\sqrt{n}\leq\frac{1-1/\sqrt{2}}{5}|\lambda_{r}^{\star}|$, we see from Weyl’s inequality (cf. Lemma 2.1.3) that

We further remark on the Euclidean statistical accuracy when estimating the unknown matrix $\bm{M}^{\star}$ using $\widehat{\bm{M}}\coloneqq\bm{U}\bm{\Lambda}\bm{U}^{\top}$, where \bm{\Lambda}\coloneqq\mathsf{diag}\big{(}[\lambda_{1},\cdots,\lambda_{r}]\big{)}. It is seen from the triangle inequality that

where the last inequality relies on (3.14). Since the rank of $\bm{U}\bm{\Lambda}\bm{U}^{\top}-\bm{M}^{\star}$ is at most $2r$, with probability at least $1-O(n^{-8})$ one has

2.3 Proof of the inequality (3.12) on ‖𝑬‖norm𝑬\|\bm{E}\|

We plan to employ Theorem 3.1.5. Given that Gaussian entries are unbounded, we introduce a truncated copy $\widetilde{\bm{E}}=[\widetilde{E}_{i,j}]_{1\leq i,j\leq n}$ defined as follows

It is readily seen from the property of Gaussian distributions that

which combined with the union bound leads to

Given that $B\coloneqq\max_{i,j}|\widetilde{E}_{i,j}|\leq 5\sigma\sqrt{\log n}$, we can invoke Theorem 3.1.5 (or more directly, (3.9)) to demonstrate that

Combining the above two observations implies that

with probability exceeding $1-O(n^{-8})$, as claimed.

3 Principal component analysis and factor models

Principal component analysis (PCA) and factor models [jolliffe1986principal, lawley1962factor, fan2020statistical]—which serve as an effective unsupervised learning tool for exploring and understanding data—arise frequently in data-intensive applications in economics, finance, psychology, signal processing, speech, neuroscience, traffic data analysis, among other things [stock2002forecasting, mccrae1992introduction, scharf1991svd, chen2015reduced, balzano2018streaming, fan2020robust]. PCA and factor models not only allow for dimensionality reduction, but also provide intermediate means for data visualization, noise removal, anomaly detection, and other downstream tasks. In this section, we investigate a simple, yet broadly applicable, factor model.

In this monograph, we concentrate on the following tractable statistical model for pedagogical reasons. See fan2020statistical for more general settings (including, say, heavy-tailed distributions and non-isotropic noise covariance matrices).

The vectors $\bm{f}_{i}$ and $\bm{\eta}_{i}$ ($1\leq i\leq n$) are all independently generated according to

the condition number of the low-rank matrix $\bm{L}^{\star}\bm{L}^{\star\top}=\bm{U}^{\star}\bm{\Lambda}^{\star}\bm{U}^{\star\top}$.

3.2 Algorithm

As a starting point, it is readily seen under Assumption 3.1 that

In brief, the covariance matrix $\bm{M}^{\star}$ is a low-rank matrix superimposed by a scaled identity matrix; for this reason, this model is also frequently referred to as the spiked covariance model [johnstone2001distribution]. The key takeaway is that the top-$r$ eigenspace of the covariance matrix $\bm{M}^{\star}$ in (3.21) coincides with the $r$-dimensional principal subspace being sought after (i.e., the one spanned by $\bm{L}^{\star}$ or $\bm{U}^{\star}$).

The above observation motivates a simple spectral algorithm, which begins by computing a sample covariance matrix

In the presence of missing data or heteroskedastic noise (meaning that the variance of the noise entries varies across different entries), the second part of the covariance matrix $\bm{M}^{\star}$ (i.e., $\sigma^{2}\bm{I}_{p}$ in (3.21)) might no longer be a scaled identity. Under such circumstances, one might need to carefully adjust the diagonal entries of $\bm{M}$ in order for the algorithm to succeed; see, e.g., lounici2014high, loh2012high, zhang2018heteroskedastic, cai2019subspace, zhu2019high, yan2021inference. The reader might consult Section 3.9 for an introduction to a commonly adopted diagonal deletion idea to address the aforementioned issue.

3.3 Performance guarantees

where $\bm{M}^{\star}$ is defined in (3.21), and

To apply the Davis-Kahan theorem, we are in need of controlling the size of the perturbation matrix $\bm{E}$. This is achieved by the following lemma, whose proof is deferred to Section 3.3.4.

Consider the settings in Section 3.3.1. Suppose that $n\geq cr\log^{3}(n+p)$ for some sufficiently large constant $c>0$. Then with probability exceeding $1-O((n+p)^{-10})$, one has

With Lemma 3.3.2 in place, we are ready to present the following theorem that controls the estimation error of the spectral algorithm.

Consider the settings in Section 3.3.1. Suppose that n\geq C\big{(}\kappa^{2}r+r\log^{2}(n+p)+\frac{\kappa\sigma^{2}p}{\lambda_{r}^{\star}}+\frac{\sigma^{4}p}{(\lambda_{r}^{\star})^{2}}\big{)}\log^{3}(n+p) for some sufficiently large constant $C>0$. Then with probability at least $1-O((n+p)^{-10})$, the following holds:

The third term $\kappa\sqrt{(r\log(n+p))/n}$ on the right-hand side of (3.25) arises due to the randomness of $\{\bm{f}_{i}\}$ but not that of $\{\bm{\eta}_{i}\}$. If our goal is instead to estimate the eigenspace of $\bm{L}^{\star}(\frac{1}{n}\sum_{i}\bm{f}_{i}\bm{f}_{i}^{\top})\bm{L}^{\star\top}$ as opposed to that of $\bm{L}^{\star}\bm{L}^{\star\top}$, then this term can be erased.

To interpret what Theorem 3.3.3 conveys, we include a few remarks in the sequel, focusing on the simple scenario where $\kappa=O(1)$. In view of Remark 3.3.4, we shall ignore the term $\kappa\sqrt{(r\log(n+p))/n}$ in the discussion below.

In comparison to the matrix denoising task (cf. Section 3.2.2) where \mathsf{dist}\big{(}\bm{U},\bm{U}^{\star}\big{)} scales linearly with the noise level $\sigma$ (cf. (3.13)), the above performance guarantees for PCA exhibit contrasting behavior in two different regimes depending on the strength of the signal-to-noise ratio (SNR), measured in terms of $\lambda_{r}^{\star}/\sigma^{2}$:

When the SNR is sufficiently large with $\lambda_{r}^{\star}/\sigma^{2}\gtrsim 1$, then the dominant factor in (3.25) is the term \sigma\big{(}{\frac{p\log(n+p)}{\lambda_{r}^{\star}n}}\big{)}^{1/2}, which scales linearly with the noise level.

When the SNR drops below the threshold $\lambda_{r}^{\star}/\sigma^{2}\lesssim 1$, then the term $\frac{\sigma^{2}}{\lambda_{r}^{\star}}\sqrt{\frac{p\log(n+p)}{n}}$—which scales quadratically with the noise level—enters the picture and becomes the dominant effect.

In truth, the quadratic term emerges since our spectral method operates upon the sample covariance matrix, which inevitably contains second moments of the noise components.

Natural questions arise as to whether the performance guarantees in Theorem 3.3.3 are tight, and whether the statistical accuracy can be further improved by designing more intelligent algorithms. These questions can be addressed by looking into the fundamental statistical limits. As established in the literature [zhang2018heteroskedastic, cai2019subspace], one cannot hope to achieve

in a minimax sense, regardless of the choice of the estimator $\widehat{\bm{U}}$; see, e.g., zhang2018heteroskedastic for a precise statement. Comparing (3.26) with Theorem 3.3.3 reveals the near statistical optimality of the spectral method (modulo some log factor), and confirms the tightness of the eigenspace perturbation theory when applied to this problem.

The Davis-Kahan sin$\bm{\Theta}$ theorem (cf. Corollary 2.3.4) thus implies that: if the perturbation size obeys $\|\bm{E}\|\leq(1-1/\sqrt{2})\lambda_{r}^{\star}$, then one has

Here, the penultimate inequality results from Lemma 3.3.2; the last line is valid as long as $n\gtrsim(\sigma^{2}/\lambda_{1}^{\star})p\log^{3}(n+p)$—a condition that would hold under the assumption of this theorem—so that the fourth term is dominated by the second one in the parenthesis of the penultimate line. Finally, it is immediately seen from Lemma 3.3.2 that the condition $\|\bm{E}\|\leq(1-1/\sqrt{2})\lambda_{r}^{\star}$ would hold under the assumption of this theorem.

3.4 Proof of Lemma 3.3.2

We start by applying the triangle inequality to (3.24) as follows

In order to develop an upper bound on this quantity, one needs to control the spectral norm of $\frac{1}{n}\bm{F}\bm{F}^{\top}-\bm{I}_{r}$, $\frac{1}{n}\bm{F}\bm{Z}^{\top}$, $\frac{1}{n}\bm{Z}\bm{F}^{\top}$ and $\frac{1}{n}\bm{Z}\bm{Z}^{\top}-\sigma^{2}\bm{I}_{p}$. All of these terms share similar randomness structure, namely, they are all averages of independent zero-mean random matrices. As a result, the truncated matrix Bernstein inequality in Corollary 3.1.3 becomes applicable. In what follows, we shall only demonstrate how to control the size of $\frac{1}{n}\bm{F}\bm{Z}^{\top}$; the other terms can be bounded similarly.

Write $\bm{F}\bm{Z}^{\top}=\sum_{i=1}^{n}\bm{f}_{i}\bm{\eta}_{i}^{\top}$. Since the entries of $\bm{F}\bm{Z}^{\top}$ might be unbounded, we start by identifying an appropriate truncation level. From standard properties about Gaussian distributions and the union bound, it is straightforward to verify that

with probability greater than $1-(n+p)^{-11.5}$. In other words, with the choice $L\coloneqq 25\sqrt{rp}\sigma\log(n+p)$ one has

Additionally, the symmetry of Gaussian distributions implies

To invoke the truncated Bernstein inequality, it remains to determine the variance statistic. Towards this end, letting ${\bm{B}}_{i}=\bm{f}_{i}\bm{\eta}_{i}^{\top}$, we observe that

where we use the fact that $r\leq p$. Taking these bound together and applying the truncated matrix Bernstein theorem (see Corollary 3.1.3) demonstrate that if $n\gtrsim r\log^{3}(n+p)$, one has

Repeating the above analysis yields that: if $n\gtrsim r\log^{3}(n+p)$, with probability at least $1-O((n+p)^{-10})$ one has

Note that we do not get rid of the second term on the right-hand side of (3.28c) since we do not assume $n\gtrsim p\log^{3}(n+p)$.

Substituting the above results (3.28) into (3.27) and recognizing the basic fact $\|\bm{L}^{\star}\|=\|\bm{U}^{\star}(\bm{\Lambda}^{\star})^{1/2}\|\leq\|(\bm{\Lambda}^{\star})^{1/2}\|=\sqrt{\lambda_{1}^{\star}}$, we conclude that

4 Graph clustering and community recovery

Next, we move on to a central problem that permeates data science applications: clustering. An important formulation that falls under this category is graph clustering or community recovery, which aims to cluster individuals into different communities based on pairwise measurements of their relationships, each of which reveals information about whether or not two individuals belong to the same community [abbe2017community]; see Figure 3.1 for an illustration. There has been a recent explosion of interest in this problem, due to its wide applicability in, say, social network analysis [azaouzi2019community], image segmentation [browet2011community], shape mapping in computer vision [huang2013consistent], haplotype phasing in genome sequencing [chen2016community], to name just a few. This section explores the capability of spectral methods in application to graph clustering; we will revisit the clustering problem again in Section 3.5 for another common formulation.

In this section, we formulate the graph clustering problem via the well-renowned stochastic block model (SBM) introduced in holland1983stochastic—an idealized generative model that commonly serves as a theoretical benchmark for evaluating community recovery algorithms.

Consider an undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ that comprises $n$ vertices, where $\mathcal{V}$ and $\mathcal{E}$ denote the vertex set and the edge set of $\mathcal{G}$, respectively. The $n$ vertices, labelled by $1,\cdots,n$, exhibit community structures and can be grouped into two non-overlapping communities of equal sizes. Here and throughout, $n$ is assumed to be an even number, so that each community contains exactly $n/2$ vertices. To encode the community memberships, we assign $n$ binary-valued variables $x_{i}^{\star}\in\{1,-1\}$ ($1\leq i\leq n$) to the vertices in a way that

The SBM assumes that the set $\mathcal{E}$ of (undirected) edges is generated randomly based on the community memberships of the incident vertices. To be precise, each pair $(i,j)$ of vertices is connected by an edge independently with probability $p$ (resp. $q$) if $i$ and $j$ belong to the same community (resp. different communities). The resultant connectivity pattern is represented by an adjacency matrix $\bm{A}=[A_{i,j}]_{1\leq i,j\leq n}\in\{0,1\}^{n\times n}$, such that for each pair $(i,j)$,

By convention, we take the diagonal entries to be $A_{i,i}=0$ for all $1\leq i\leq n$. As a remark, the matrix $\bm{A}$ is symmetric since $\mathcal{G}$ is an undirected graph, with upper triangular elements being realizations of independent Bernoulli random variables with mean either $p$ (if two nodes are in the same community) or $q$ (otherwise). In addition, it is assumed throughout that $p>q>0$, implying that there are in expectation more within-community edges than across-community edges.

Based on the adjacency matrix $\bm{A}$ generated by the SBM, the goal is to identify the latent community memberships of the vertices. To phrase it in mathematical terms, the aim is to reconstruct the vector $\bm{x}^{\star}=[x_{i}^{\star}]_{1\leq i\leq n}\in\{1,-1\}^{n}$ modulo the global sign, namely, recovering either $\bm{x}^{\star}$ or $-\bm{x}^{\star}$. This is all one can hope for, as there is absolutely no basis to distinguish the names of two groups.

4.2 Algorithm: spectral clustering

Now we describe a spectral method. To simplify presentation, it is assumed without loss of generality that: $x_{i}^{\star}=1$ for any $1\leq i\leq n/2$, and $x_{i}^{\star}=-1$ for any $i>n/2$.

A starting point for the algorithm design is to examine the mean of the adjacency matrix, given as follows

As revealed by the above calculation, the matrix constructed below

exhibits an approximate rank-1 structure, in the sense that its mean

is a rank-1 matrix. The leading eigenvalue of $\bm{M}^{\star}$ and its associated eigenvector are given respectively by

Crucially, the eigenvector $\bm{u}^{\star}$ encapsulates the precise community structure we seek to recover: all positive entries of $\bm{u}^{\star}$ correspond to vertices from one community, while the remaining ones form another community.

Inspired by the above calculation, a candidate spectral clustering algorithm consists of eigendecomposition followed by entrywise rounding:

Compute the leading eigenvector $\bm{u}$ of $\bm{M}$ (constructed in (3.30));

Compute the estimate $\bm{x}=[x_{i}]_{1\leq i\leq n}$ such that for any $1\leq i\leq n$,

In words, the community memberships are estimated in accordance with the signs of the entries of the leading eigenvector of $\bm{M}$, namely, the entries with the same signs are declared to come from the same cluster.

4.3 Performance guarantees: almost exact recovery

Consider the settings in Section 3.4.1, and suppose that $np\gtrsim\log n$. Then with probability at least $1-O(n^{-8})$, one has

This spectral norm bound, in conjunction with the Davis-Kahan sin$\bm{\Theta}$ theorem, leads to the following theoretical support for the spectral method introduced in Section 3.4.2.

Consider the setting in Section 3.4.1, and suppose that

With probability exceeding $1-O(n^{-8})$, the spectral method achieves

can be understood as the mis-clustering rate. In a nutshell, Theorem 3.4.3 asserts that with the assistance of simple rounding (i.e., the $\mathsf{sgn}(\cdot)$ operation), the spectral method allows for almost exact community recovery—namely, correctly clustering all but a vanishing fraction of the vertices—assuming satisfaction of Condition (3.40). Note that “almost exact recovery” is also referred to as “weak consistency” in the literature [abbe2017community].

Let us take a moment to interpret the recovery condition in (3.40). The first requirement in Condition (3.40) ensures the presence of sufficiently many edges in the observed graph, while still permits the graph to be fairly sparse (with average vertex degrees as low as the order of $\log n$). The second requirement in Condition (3.40)—which imposes a lower bound on the separation between the edge densities $p$ and $q$—guarantees that the within-community edges can be adequately differentiated from across-community edges. As a more concrete example, consider the scenario where $p\asymp(\log n)/{n}$ (so that each vertex is only expected to be incident to $O(\log n)$ edges). In this case, the second requirement in Condition (3.40) can be translated into

This indicates that the separation $p-q$ is allowed to be considerably smaller than the edge densities, even in this low-edge-density regime. In comparison, in another extreme case with $p\asymp 1$ (so that each vertex is likely to be connected with a constant fraction of other vertices), the second requirement in Condition (3.40) reads

thereby allowing the edge density difference to be even $\sqrt{n}$ times smaller than the edge densities themselves.

It is readily seen from Lemma 3.4.2 that with with probability at least $1-O(n^{-8})$,

provided that Condition (3.40) holds. Here, $\lambda^{\star}$ is defined in (3.37). Apply Corollary 2.3.4 to yield that with probability at least $1-O(n^{-8})$,

where the last relation follows from Condition (3.40).

Assume, without loss of generality, that \|\bm{u}-\bm{u}^{\star}\|_{2}=\mathsf{dist}\big{(}\bm{u},\bm{u}^{\star}\big{)}. We shall pay attention to the set

which in turn leads to the advertised result

4.4 Proof of Lemma 3.4.2

In addition, the variance of $E_{i,j}$ is upper bounded by

for any $(i,j)$, where (i) follows since $A_{i,j}$ is a Bernoulli random variable with mean either $p$ or $q$, and (ii) is due to the assumption $p>q$. The bound (3.9) and the condition $np\gtrsim\log n$ thus imply that

with probability exceeding $1-O(n^{-8})$.

5 Clustering in Gaussian mixture models

This section is also concerned with clustering, with the aim of grouping unlabeled data points into a few clusters (so that the data within the same cluster share similar characteristics). In contrast to the graph clustering setting in Section 3.4 where only pairwise measurements are available, this section assumes direct access to data samples for each individual. Spectral methods—possibly with the aid of subsequent refinement like $k$-means—continue to be remarkably effective for this setting, achieving practical success in, say, image segmentation [shi2000normalized], text separation [reynolds1995robust], climate modeling [lin2017statistical], and heterogeneity modeling in precision medicine and marketing [fan2014challenges]. Motivated by the empirical successes, understanding the theoretical properties of spectral clustering has garnered growing attention recently. In particular, Gaussian mixture models emerge as a succinct model of attack, providing elegant yet intuitive abstractions to pivotal quantities that dictate the feasibility of spectral clustering.

where the noise vector $\bm{\eta}_{i}\sim\mathcal{N}(\bm{0},\bm{I}_{p})$ is independently generated across the samples. In words, $\xi_{i}^{\star}$ indicates which Gaussian component a sample is generated from. Clustering in this Gaussian mixture model can, therefore, be posed as recovering the set of cluster membership variables $\{\xi_{i}^{\star}\}_{1\leq i\leq n}$ (modulo the global permutation ambiguity).

To simplify our exposition, we impose the following assumptions throughout this section. As a worthy note, this assumption is often non-essential and can be significantly relaxed, which we shall remark on momentarily in Remark 3.5.5.

The centers are independently generated obeying

hold for any pair $i\neq j$, where $o(1)$ denotes a vanishingly small quantity as $p$ approaches infinity. The indication is that the parameter $\Delta$ reflects (approximately) the separation between any pair of centers.

For simplicity of presentation, it is further assumed that there are exactly $n/r$ samples drawn from each of the $r$ Gaussian components. Without loss of generality, we assume that

5.2 Algorithm and rationale

In order to develop a spectral clustering algorithm, it is instrumental to first examine the spectral feature of the following data matrix

Similarly, the Gram matrix $\bm{X}^{\top}\bm{X}$ also inherits this rank-$r$ structure in the following sense (albeit in the form of a “spiked” structure due to the presence of noise):

Recognizing that $\bm{F}^{\star}$ encodes all the cluster membership information, one is motivated to attempt information extraction from the rank-$r$ eigenspace of $\bm{X}^{\top}\bm{X}$, akin to the PCA algorithm introduced in Section 3.3.2.

With the preceding spectral properties in mind, we are ready to present a spectral clustering algorithm tailored to this Gaussian mixture model. Given that the eigenspace of $\bm{X}^{\top}\bm{X}$ might only approximate $\bm{F}^{\star}$ up to global rotation, we include a follow-up $k$-means scheme [macqueen1967some] to produce a valid clustering outcome based on the spectral estimate.

for any matrix $\bm{Z}=[\bm{z}_{1},\cdots,\bm{z}_{n}]$. As will be discussed below, the projection step is not necessary, and we can also simply take ${\bm{Y}}={\bm{U}}{\bm{U}}^{\top}$.

Let $\bm{y}_{i}$ represent the $i$-th column of $\bm{Y}$, and apply the $k$-means algorithm (with $k=r$) to the vectors $\{{\bm{y}}_{i}\}_{1\leq i\leq n}$ to find the cluster centers and cluster labels for all individuals; namely, we compute