Near-Optimal Performance Bounds for Orthogonal and Permutation Group Synchronization via Spectral Methods

Introduction

Suppose there are $n$ group elements $\{g_{i}\}_{i=1}^{n}\in{\cal G}$ and we observe their noisy pairwise measurements

where $w_{ij}$ is the noise and ${\cal E}$ is the edge set of an underlying network. How to recover these elements $g_{i}$ from the noisy observations $\{g_{ij}\}_{(i,j)\in{\cal E}}$? Depending on the specific group type, the group synchronization problem is widely used in many applications including computer vision , robotics , clock synchronization and cryo-electron microscopy . In this paper, we will focus on the synchronization of the orthogonal and permutation group.

The group ${\cal G}$ in (1.1) is the orthogonal group $\text{O}(d)$,

The general $\text{O}(d)$ synchronization includes $\hbox{\msbm{Z}}_{2}$-synchronization ($d=1$), angular synchronization $(d=2)$, SO(3) synchronization as special cases . It often arises in rotation estimation and structure-from-motion , and also plays a significant role in SLAM (simultaneous localization and mapping) in robotics .

Permutation group synchronization:

The underlying group ${\cal G}$ in (1.1) becomes permutation group, which is represented by permutation matrices $\Pi_{d}$:

Essentially, the permutation group synchronization is a special case of the $\text{O}(d)$ synchronization since $\Pi_{d}$ is a subgroup of $\text{O}(d).$ Permutation group is directly related to the multi-way matching problem (map synchronization) in computer vision. Suppose there are $n$ images of the same object and each of them has $d$ features. Given a set of partially known feature correspondence among these $n$ images, how to find the all the correct pairwise bijection? This matching problem is one of the core problems in image registration, structure from motion, and object matching problem . This multi-way matching problem can be reformulated as recovering a set of permutation matrices from their pairwise products where each bijection corresponds to a permutation matrix .

Note that every element $\bm{R}$ in $\text{O}(d)$ or $\Pi(d)$ satisfies $\bm{R}^{-1}=\bm{R}^{\top}$. Therefore, the general synchronization problem reduces to recovering $n$ group elements $\{\bm{G}_{i}\}_{i=1}^{n}$ from its noisy measurements

where the edge set ${\cal E}$ is assumed to be a complete graph throughout this manuscript. From now on, we let $\bm{A}_{G}:=[\bm{G}_{ij}]_{1\leq i,j\leq n}$ be the $nd\times nd$ data matrix.

Given its practical importance, many efforts have been taken to solve the group synchronization problem. In absence of noise, group synchronization is easily solvable by sequentially recovering the group elements. However, this sequential strategy no longer works in presence of noise since the noise will be amplified. One common approach is to find the least squares estimator. However, it is usually an NP-hard problem to obtain the least squares estimator exactly, even for the simplest group $\hbox{\msbm{Z}}_{2}=\{1,-1\}.$ As a result, many optimization approaches, including convex relaxation and nonconvex methods, are developed to tackle various challenging scenarios. In this work, we will instead focus on the spectral methods for orthogonal/permutation group synchronization. There are several variants of spectral methods for $\text{O}(d)$ and $\Pi(d)$ group synchronization which are based on the data matrix $\bm{A}_{G}:=[\bm{G}_{ij}]_{1\leq i,j\leq n}$ or its corresponding (normalized) connection Laplacian matrix . Here we will focus the spectral methods which begin with computing the top $d$ eigenvectors of the observed data $\bm{A}_{G}$ and then approximate each group element by rounding all the $d\times d$ blocks of the eigenvectors. In particular, we will investigate its performance and answer the following questions:

1 Related works and our contribution

Group synchronization has found many applications in signal processing, computer vision, and machine learning. Some prominent examples include community detection ($\hbox{\msbm{Z}}_{2}$ synchronization), joint alignment (finite cyclic group $\hbox{\msbm{Z}}_{n}$), angular synchronization , statistical ranking and phase retrieval (unitary group U(1)), object matching (permutation group), rotation estimation (SO(3) group), clock synchronization (cyclic group on a finite interval), and simultaneous localization and mapping (SLAM) in robotics (special Euclidean group SE$(d)$). There have been many efforts on solving the group synchronization problem in different settings by using various approaches including optimization-based approach , spectral methods , and message-passing type methods . For general group synchronization, one important topic is to determine how the noise strength affects the performance of algorithms and solvability. The fundamental recovery criterion for information recovery from pairwise measurements is studied . The accuracy and noise sensitivity of the spectral method for general compact groups are presented in . From now on, we will briefly review the recent literatures on orthogonal and permutation group synchronization and highlight those works which motivate this work.

Orthogonal group synchronization is often considered in rotation estimation arising from computer vision and robotics. One of the most widely used approaches to tackle general $\text{O}(d)$ synchronization is to find the least squares estimator. As pointed out before, it is often an NP-hard problem to find the least squares estimator since the objective function is usually highly nonconvex and even discrete in some cases. This poses a significant challenge to practical implementation. One important idea to overcome this technical difficulty is to find appropriate relaxations which are solvable within polynomial time. Convex relaxation has proven to be a very powerful method . However, the solution to the convex relaxation program is not necessarily equal to that of the original program, i.e., the tightness does not always hold. The study of the tightness of convex relaxation has been a research focus in orthogonal group synchronization. In , Wang and Singer investigated the semidefinite program (SDP) relaxation of the orthogonal group synchronization under random corruption and characterized the phase transition of group recovery from noisy measurements. The tightness of the SDP relaxation for angular synchronization, as a special case of $\text{O}(d)$ synchronization, is studied in with a near-optimal performance bound on the signal-to-noise ratio introduced in the very inspiring work . Recent works propose suboptimal deterministic conditions which guarantee the tightness of the SDP relaxation for general $\text{O}(d)$ synchronization. A similar route of research can also be found for permutation group synchronization. Huang and Guibas studied the convex relaxation approach of the permutation group synchronization in and provided theoretical guarantees for correct recovery. The work investigated exact and robust object matching via SDP relaxation under partially known similarity between objects and the performance bound is near-optimal up to a log-factor. Despite the usefulness of convex relaxation, it remains highly nontrivial to solve large-scale SDPs. In practice, efficient first-order gradient-based approaches are preferred such as Riemannian optimization , the Burer-Monteiro factorization , and iterative reweighing strategy . The major issue of Riemannian optimization is the inherent nonconvexity of the objective function, which could potentially create local optima. Fortunately, we have seen a surge of research in exploring the provably convergent nonconvex methods in solving $\hbox{\msbm{Z}}_{2}$ synchronization , angular synchronization , permutation group synchronization , and $\text{O}(d)$ or SO$(d)$ synchronization in .

Our contribution consists of several aspects: we first study the spectral methods for $\text{O}(d)$ synchronization under Gaussian noise: namely first computing the top $d$ eigenvectors of $\bm{A}_{G}$ and use them to estimate the group elements. We provide a block-wise near-optimal error bound for each group element (modulo a constant) which justifies the usefulness of the spectral methods in $\text{O}(d)$ synchronization. This analysis can be regarded as a natural generalization from $\hbox{\msbm{Z}}_{2}$-synchronization in and angular synchronization in . Then we study the permutation group synchronization under uniform random corruption. We are interested in when the two-step approach, namely, eigenvectors followed by rounding procedure, can give the exact recovery of the planted permutation matrix. The derived bound is also nearly optimal in terms of information theoretical limits, and improves the bound in and matches the bounds obtained via the SDP relaxation in . It is well worth noting that studies a more general setting of permutation group synchronization and provides a near-optimal performance bound for the spectral methods. However, the technical approach is quite different from ours. Our theory is developed by applying the recent popular leave-one-out technique. However, the block-wise analysis of eigenvectors requires additional technical treatments. Our work resolved one question raised in about the block-wise analysis of eigenvectors for matrices with row/column block-wise independence. This framework is quite flexible and can be applied to other problems which require the block-wise analysis of eigenvectors.

2 Organization

Section 2 introduces the mathematical models and the spectral methods for group synchronization. We present the main results in Section 3 with numerical experiments to support our theory in Section 4. The proofs are provided in Section 5.

3 Notation

Given a matrix $\bm{X}$, $\bm{X}^{\top}$ is the transpose of $\bm{X}$ and $\bm{X}\succeq 0$ means $\bm{X}$ is positive semidefinite. $\bm{I}_{n}$ is the $n\times n$ identity matrix, $\bm{J}_{n}$ is the $n\times n$ “1” matrix, and $\bm{1}_{n}$ is an $n\times 1$ “1” vector. $\|\bm{X}\|$ denotes the operator norm of $\bm{X}$ and $\|\bm{X}\|_{F}$ is the Frobenius norm. For two matrices $\bm{X}$ and $\bm{Y}$, we denote $\bm{X}\otimes\bm{Y}$ their Kronecker product, i.e., the $(i,j)$-block of $\bm{X}\otimes\bm{Y}$ is $X_{ij}\bm{Y}$. For a matrix $\bm{X}$, we let $\sigma_{i}(\bm{X})$ and $\lambda_{i}(\bm{X})$ be the $i$th largest singular value and eigenvalue of $\bm{X}$ respectively. For two nonnegative functions $f(n)$ and $g(n)$, we denote $f(n)\lesssim g(n)$ and $f(n)=O(g(n))$ if there exists an absolute positive constant $C$ such that $f(n)\leq Cg(n)$ for all $n$.

Preliminaries

This paper will study two benchmark models of group synchronization under additive noise and uniform corruption (multiplicative noise).

Orthogonal group synchronization under additive Gaussian noise. The pairwise noisy measurement $\bm{G}_{ij}$ is observed between $\bm{G}_{i}$ and $\bm{G}_{j}$,

where $\bm{G}_{i}\in\text{O}(d)$ and $\bm{W}_{ij}\in\hbox{\msbm{R}}^{d\times d}$ is a Gaussian random matrix.

Permutation group synchronization under uniform random corruption. Consider

where $\{\bm{G}_{i}\}_{i=1}^{n}$ are the hidden permutation matrices and $\bm{P}_{ij}\in\hbox{\msbm{R}}^{d\times d}$ is an independent random permutation uniformly sampled from $d!$ permutation matrices. In other words,

where $X_{ij}\sim$Bernoulli($p$) is independent of $\bm{P}_{ij}.$

For both models, our goal is to recover $\bm{G}_{i}$ from the noisy measurements $\bm{G}_{ij}.$ One common method is to find the least squares estimator by minimizing

whose global minimizer equals the global maximizer of the following generalized quadratic form:

However, it is in general NP-hard to find the global optimizer. Therefore, one wants to find an appropriate relaxation of (2.1). The idea of spectral relaxation uses a simple fact: by letting $\bm{R}$ be an $nd\times d$ matrix whose $i$th block equals $\bm{R}_{i}$, then (2.1) is equivalent to

where $\bm{A}_{G}$ is an $nd\times nd$ symmetric matrix whose $(i,j)$-block is $\bm{G}_{ij}$. Note that all $\bm{R}\in\text{O}(d)^{\otimes n}$ satisfies $\bm{R}^{\top}\bm{R}=n\bm{I}_{d}$. The spectral method simply replaces the constraints $\bm{R}\in\text{O}(d)^{\otimes n}$ by $\bm{R}^{\top}\bm{R}=n\bm{I}_{d}$,

whose global maximizer equals the top $d$ eigenvectors of $\bm{A}_{G}.$

As a result, the spectral method is very convenient to use: simply compute the top $d$ eigenvectors of the matrix $\bm{A}_{G}$, denoted by an $nd\times d$ partial orthogonal matrix $\bm{\Phi}$ where $\bm{\Phi}^{\top}=[\bm{\Phi}_{1}^{\top},\cdots,\bm{\Phi}_{n}^{\top}]$ and $\bm{\Phi}_{i}$ is the $i$th $d\times d$ block. In particular, we normalize $\bm{\Phi}$ to be $\bm{\Phi}^{\top}\bm{\Phi}=n\bm{I}_{d}$, i.e., each column is of norm $\sqrt{n}$. Then we implement a rounding procedure to obtain the estimation of $\bm{G}_{i}.$ We summarize the aforementioned procedures in Algorithm 1.

For permutation matrix, a slight modification of the rounding procedure is implemented. Simply speaking, once we get $\bm{\Phi}_{i}$, we estimate $\bm{G}_{i}$ via

where $\Pi_{d}$ is the set of all $d\times d$ permutation matrices. This linear assignment problem can be solved by the Hungarian algorithm in polynomial time . The entire procedures are summarized in Algorithm 2.

How well do these algorithms work? We consider the $\text{O}(d)$ synchronization under additive Gaussian noise as an example. The data matrix $\bm{A}_{G}$ can be naturally written into a spiked matrix model: $\bm{A}_{G}=\bm{G}\bm{G}^{\top}+\bm{\Delta}$ where $\bm{\Delta}=\sigma\bm{W}$. Note that without any noise, the top $d$ eigenvectors exactly give the group elements. If the noise $\bm{\Delta}$ is small, then one can easily invoke the classical matrix perturbation argument, e.g. Davis-Kahan theorem (Theorem 6.2), to obtain an error bound between the top $d$ eigenvectors and the planted group elements in terms of operator or Frobenius norm. Namely,

which will be derived more carefully later in the proof section.

On the other hand, it is much more appealing to provide an error bound for

for some orthogonal matrix $\bm{Q}\in\hbox{\msbm{R}}^{d\times d}$ since this would provide us an error bound for the recovery of each group element. In other words, we need to control the estimation error for each block $\bm{\Phi}_{i}$, which is essentially a generalization of the entrywise bound for the eigenvector discussed in . However, the Davis-Kahan bound does not immediately yield a tight bound for the deviation of each $\bm{\Phi}_{i}$ from $\bm{G}_{i}\bm{Q}$ for some $\bm{Q}\in\text{O}(d)$. This will be the main focus of our paper: we obtain the block-wise perturbation bound of $\bm{\Phi}$ via the leave-one-out technique. We will introduce this technique briefly in Section 3.2 and provide more details in Section 5.

Main theorem

In this section, we will provide theoretical guarantees for the spectral methods in solving the $\text{O}(d)$ and $\Pi_{d}$ synchronization problem under the statistical models (OD) and (PM) respectively.

Our main contribution is providing a near-optimal block-wise error bound of $\widehat{\bm{G}}_{i}$ for all $1\leq i\leq n$. For the $\text{O}(d)$ synchronization under Gaussian noise, we have the following theorem.

Suppose the parameter $\sigma$ in the model (OD) satisfies

for some small constant $c_{0}>0.$ Then with high probability, the estimation $\widehat{\bm{G}}_{i}$ of $\bm{G}_{i}$ from Algorithm 1 satisfies

by letting $\bm{Q}_{j}=\bm{G}_{j}^{\top}\widehat{\bm{G}}_{j}$ for any $1\leq j\leq n.$

Theorem 3.1 includes $\hbox{\msbm{Z}}_{2}$- and angular synchronization as special cases. In particular, if $d=1$, the problem reduces to $\hbox{\msbm{Z}}_{2}$-synchronization and the bound is equivalent to the one derived in ; for $d=2$, our result is closely related to the angular synchronization explored in since SO(2) is isomorphic to U(1).

Simply speaking, Theorem 3.1 provides a theoretical guarantee for the spectral estimator in the orthogonal group synchronization under additive Gaussian noise: the distance of the spectral estimator from the planted signal is controlled by the noise strength. As discussed before, the spectral methods are viewed as a relaxation of the equivalent least squares objective function (2.1). Therefore, they are unlikely to produce the globally optimal least squares estimator (2.1). However, the proximity of the spectral estimator to the ground truth provides allows nonconvex optimization approaches to have a high-quality initialization and enjoy a global convergence to the globally optimal least squares estimator .

Now we briefly discuss the optimality of our result. Note that the model $\bm{A}_{G}=\bm{G}\bm{G}^{\top}+\sigma\bm{W}$ for the $\text{O}(d)$ synchronization under additive Gaussian is essentially the well-known spiked matrix model or the real deformed Wigner matrices . Note that in random matrix theory, it has been extensively studied when the top eigenvectors of $\bm{A}_{G}$ are correlated with the planted signals (low-rank matrix), see e.g. . For this finite-rank spiked matrix model, it has been shown in if the noise level $\sigma$ is above the threshold $\sigma>\sqrt{n/d}$, the leading $d$ eigenvalues of $\bm{A}_{G}$ fail to exit the limiting semicircle compact support of the GOE (Gaussian orthogonal ensemble) for a sufficiently large $n$. This implies the spectral method (plus rounding) is expected to identify the planted signal only in the regime $\sigma\lesssim\sqrt{n/d}$. Thus our bound in Theorem 3.1 differs from this threshold only by a logarithmic and constant factor. Though not explicitly stated, it is believed that $\sigma=\sqrt{n/d}$ is the threshold above which is information-theoretically possible to detect the spikes .

The theoretical result for permutation group synchronization is summarized as follows.

Suppose the parameter $p$ in the model (PM) satisfies

for some universal large constant $C_{0}>0$. Then with high probability, the estimation $\widehat{\bm{G}}_{i}$ of $\bm{G}_{i}$ from Algorithm 2 satisfies

In particular, if $\|\widehat{\bm{G}}_{i}\widehat{\bm{G}}_{j}^{\top}-\bm{G}_{i}\bm{G}_{j}^{\top}\|<\frac{1}{2}$, then Algorithm 2 recovers the hidden permutation matrices $\bm{G}_{i}$ exactly.

The work provides a block-wise bound for permutation group synchronization on general networks in which $p>C_{0}n^{-\frac{1}{2}}\log^{3}(n)$ is needed for the exact recovery of all the permutation matrices with high probability. shows that the SDP relaxation can recover the underlying hidden permutation matrices with high probability if $p>C_{0}n^{-\frac{1}{2}}\log^{2}(nd)$. The bound (3.1) matches the state-of-the-art performance bound in which considers the general simultaneous mapping and clustering problem. However, as pointed out earlier, our technique is quite different from . Note that the information theoretic limit for the exact recovery in $\Pi_{d}$ synchronization is discussed in [21, Corollary 1]: no method whatsoever is able to recover the ground truth if $p<O(1/\sqrt{n})$. Therefore, our bound differs from the information-theoretic limit by a logarithmic factor.

Another synchronization model which is highly relevant to the two aforementioned models is the $\text{O}(d)$ group synchronization with uniform multiplicative noise :

where $\bm{R}_{ij}$ is sampled from the uniform Haar distribution over $\text{O}(d)$Simply speaking, Haar distribution on $\text{O}(d)$ is the unique invariant probability measure on the compact group $\text{O}(d)$.. Though it is not analyzed in our manuscript, the proof technique for the permutation group synchronization under uniform corruption could be directly modified to tackle this $\text{O}(d)$ synchronization under uniform multiplicative corruption.

2 The sketch of proof: leave-one-out technique

We provide a proof sketch for Theorem 3.1 and 3.3, and will proceed to give more technical details in Section 5. The main idea follows from the leave-one-out technique employed in to study $\hbox{\msbm{Z}}_{2}$-synchronization and community detection under the stochastic block model. The major difference of our setting here is the blockwise independence of the noise matrix as well as the multi-dimensionality of the eigenspace, which requires additional technical treatments.

With a bit of calculation, both (OD) and (PM) can be formulated under the framework of the spiked matrix model. Without loss of generality, we assume each $\bm{G}_{i}$ is an identity matrix $\bm{I}_{d}$ and it suffices to consider

where $\bm{Z}^{\top}=[\bm{I}_{d},\cdots,\bm{I}_{d}]\in\hbox{\msbm{R}}^{d\times nd}$. Here $\bm{\Delta}$ is the random noise matrix. More precisely,

For model (OD), the noise matrix $\bm{\Delta}$ is

where $\bm{W}\in\hbox{\msbm{R}}^{nd\times nd}$ is a symmetric Gaussian random matrix.

For model (PM), the corruption matrix $\bm{\Delta}$ is

where $\bm{P}_{ij}$ is a random permutation matrix drawn uniformly from the set of all $d\times d$ permutation matrices and $\bm{J}_{d}$ is a $d\times d$ matrix whose entries are all equal to 1. In fact, the top $d$ eigenvectors of $\bm{A}$ and $\bm{A}_{G}=[\bm{G}_{ij}]_{1\leq i,j\leq n}$ are the same provided that the noise $\bm{\Delta}$ is small and $\bm{G}_{i}=\bm{I}_{d}$. We will justify this fact in Lemma 5.9 in Section 5.3.

Now we briefly introduce the main idea of the leave-one-out technique in obtaining a block-wise error bound for the top $d$ eigenvectors $\bm{\Phi}$ of $\bm{A}$. Assume $(\bm{\Phi},\bm{\Lambda})$ is the top $d$ leading eigen-pairs of $\bm{A}$, i.e.,

In other words, it holds $\bm{\Phi}=\bm{A}\bm{\Phi}\bm{\Lambda}^{-1}.$ The idea of estimating each $\bm{\Phi}_{i}$ relies on choosing a suitable surrogate which is easy to approximate and also close to $\bm{\Phi}_{i}$. One commonly-used choice is to use one-step fixed point iteration, which is inspired by . By definition, $\bm{\Phi}\in\hbox{\msbm{R}}^{nd\times d}$ is the fixed point of the following map:

where $\bm{\Lambda}\in\hbox{\msbm{R}}^{d\times d}$ consists of the top $d$ eigenvectors of $\bm{A}.$

Note that the recovered orthogonal group $\{\bm{G}_{i}\}_{i=1}^{n}$ is unique modulo a global rotation. Therefore, we initialize this fixed point map (3.4) by choosing $\bm{X}=\bm{Z}\bm{Q}$ where $\bm{Q}$ minimizes the distance $d_{F}(\bm{\Phi},\bm{Z})$ between $\bm{\Phi}$ and $\bm{Z}$ is minimized, i.e.,

We hope $f(\bm{Z}\bm{Q})$ is close to $\bm{\Phi}$ uniformly for each $d\times d$ block. Let’s perform a preliminary analysis for the approximation error bound of $\bm{\Phi}_{i}$ with the $i$th block of $f(\bm{Z}\bm{Q})=\bm{A}\bm{Z}\bm{Q}\bm{\Lambda}^{-1}.$ Let $\bm{\Delta}_{i}$ be the $i$th block column of $\bm{\Delta}$, and then it holds

The Davis-Kahan theorem provides a tight bound of the first term and the goal is to estimate $\|\bm{\Delta}_{i}^{\top}(\bm{\Phi}-\bm{Z}\bm{Q})\|$. Note that $\bm{\Delta}_{i}$ and $\bm{\Phi}-\bm{Z}\bm{Q}$ are not statistically independent. Therefore, despite that each $\bm{\Delta}_{ij}$ is either a Gaussian random matrix or a bounded centered random permutation matrix, we cannot immediately apply concentration inequality to obtain a tight bound of $\|\bm{\Delta}_{i}^{\top}(\bm{\Phi}-\bm{Z}\bm{Q})\|$. The remedy is to use the recently popular leave-one-out trick.

The idea is to replace $\bm{\Phi}$ by $\bm{\Phi}^{(i)}$ which is the top $d$ eigenvectors of the following auxiliary matrix $\bm{A}^{(i)}=\bm{Z}\bm{Z}^{\top}+\bm{\Delta}^{(i)}$:

In other words, $\bm{A}$ and $\bm{A}^{(i)}$ only differ by the $i$th block column and row of $\bm{\Delta}.$ Because of this minor difference, the corresponding eigenspace $\bm{\Phi}$ and $\bm{\Phi}^{(i)}$ are very close. More importantly, $\bm{\Phi}^{(i)}$ is independent of $\bm{\Delta}_{i}$ since $\bm{\Phi}^{(i)}$ only depends on $\bm{\Delta}^{(i)}$ which excludes $\bm{\Delta}_{i}$. This important fact allows one to apply the concentration inequality to get a satisfactory bound of $\|\bm{\Delta}_{i}^{\top}(\bm{\Phi}-\bm{Z}\bm{Q})\|$ which will be discussed in more details in Section 5.

Numerics

In this section, we will provide numerical evidence to show that the bound in the main theorems are near-optimal.

We first investigate the performance of Algorithm 1 under various noise levels. Consider $\bm{A}_{G}=\bm{G}\bm{G}^{\top}+\sigma\bm{W}$ where $\bm{W}$ is a symmetric $nd\times nd$ Gaussian random matrix where $n=1000$ and $d=2,3,5$ and 10. Here we introduce another parameter $\kappa$ such that $\sigma=\kappa\sqrt{n/d}$ because Theorem 3.1 implies

would provide a non-trivial bound for some constant $c_{0}.$ We let $\kappa$ vary from $0.05$ to 0.5 since we know that the spectral method is expected to fail for $\kappa>1$ and to succeed for $\kappa<1$ based on the random matrix theory. Here once we obtain $\widehat{\bm{G}}_{i}$, we calculate the maximum blockwise deviation of $\widehat{\bm{G}}_{i}$ from $\bm{G}_{i}$ by using

For each $\kappa$, we run 25 experiments and obtain the boxplot of error, as is shown in Figure 1. The bottom/top edges of each blue box stand for the 25th and 75th percentile of the estimation error, the central red mark indicates the median error, and each red dot is an outlier. We can see that the median block-wise error grows approximately linearly with respect to $\kappa$ (equivalently, $\sigma$) if $\kappa$ is smaller than some threshold for each fixed $d$. Once $\kappa$ exceeds that threshold, the algorithm cannot provide a nontrivial estimation of $\bm{G}_{i}$ as the error becomes 2. In addition, it is also interesting to see that the range of $\kappa$ for a nontrivial error bound (i.e., the threshold) gets larger as $d$ increases. This may be explained by the inequality above (4.1): as $d$ gets larger, the denominator $1+\sqrt{d^{-1}\log n}$ in (4.1) becomes smaller and thus allows a larger range of $\kappa$ for a non-trivial error bound.

2 Permutation group synchronization

We consider the permutation group synchronization under random corruption. Without loss of generality, we assume $\bm{G}_{i}=\bm{I}_{d}$ and $\bm{G}_{ij}=X_{ij}\bm{I}_{d}+(1-X_{ij})\bm{P}_{ij}$ as introduced in (PM2). Then we compute the top $d$ eigenvectors and estimate $\widehat{\bm{G}}_{i}$ by solving

The goal is to study how the performance of Algorithm 2 depends on the parameters $(n,p)$. Note that our theorem indicates that the algorithm works if $p$ is larger than $\sqrt{n^{-1}\log(nd)}$. Therefore, we introduce the parameter $\kappa$ so that

We let $n$ vary from 50 to 1000, and $\kappa$ between $0.1$ and 2. We use two ways of measuring the recovery performance.

Exact recovery: We compute the total number of instances in which $\widehat{\bm{G}}_{i}$ equals $\bm{G}_{i}$ for all $1\leq i\leq n$. For each pair of $(n,\kappa)$, we run 25 experiments and calculate the proportion of successful instances. Figure 2 implies that for $\kappa>1$, the exact recovery holds with high probability for both $d=5$ and 10. This confirmed the near-optimality of our performance bound in Theorem 3.3.

Weak recovery: Instead of looking at the exact recovery of all the permutation matrices, we compute the average alignment to see if most permutation matrices are recovered even if the assumption of Theorem 3.3 is violated, i.e., adding slightly more corruption to the observed data. The average alignment associated with $\{\widehat{\bm{G}}_{i}\}$ is defined as

which is a number between 0 and 1. We simulate 25 instances and the compute the mean of the average alignment. The phase transition plot is provided in Figure 3, which shows that if $\kappa>0.6$, the recovered permutation matrix is highly aligned with the planted signal. We leave the characterization of the critical threshold for the exact/weak recovery to the future work.

Proofs

This section is devoted to the proof of Theorem 3.1 and 3.3.

As discussed in Section 3.2, the blockwise estimation of $\bm{\Phi}$ is given by (3.6):

where $(\bm{\Phi},\bm{\Lambda})$ are the top $d$ eigen-pairs of $\bm{A}$ satisfying $\bm{A}\bm{\Phi}=\bm{\Phi}\bm{\Lambda}$ and $\bm{A}=\bm{Z}\bm{Z}^{\top}+\bm{\Delta}$ is given in (3.2). We have shown that the key is to control $\bm{\Delta}_{i}^{\top}(\bm{\Phi}-\bm{Z}\bm{Q})$ by using the leave-one-out technique and it is much easier to bound $\|\bm{Z}^{\top}(\bm{\Phi}-\bm{Z}\bm{Q})\|$ as well as the spectra of $\bm{\Lambda}.$

Now we will proceed to estimate $\|\bm{\Delta}_{i}^{\top}(\bm{\Phi}-\bm{Z}\bm{Q})\|$ by approximating it with $\|\bm{\Delta}_{i}^{\top}(\bm{\Phi}^{(i)}-\bm{Z}\bm{Q})\|$ where $\bm{\Phi}^{(i)}$ is the top $d$ eigenvectors from the auxiliary matrix $\bm{A}^{(i)}$ defined in (3.7). In addition, we also need to bound the distance between $\bm{\Phi}$ and $\bm{\Phi}^{(i)}$. To do that, we formally define an operator which will be frequently used in the discussion. One can view it as a generalization of taking the “phase” of a matrix which is also known as the matrix sign function .

For any matrix $\bm{X}\in\hbox{\msbm{R}}^{m\times r}$ with $m\geq r$, we define

where $\bm{X}=\bm{U}\bm{\Sigma}\bm{V}^{\top}$ is the SVD (singular value decomposition) of $\bm{X}$ with $\bm{U}^{\top}\bm{U}=\bm{V}^{\top}\bm{V}=\bm{I}_{r}$. In particular, if $\operatorname{rank}(\bm{X})=r$, i.e., $\bm{\Sigma}$ is invertible, then

Note that if the input matrix $\bm{X}$ is not full rank, i.e., $\operatorname{rank}(\bm{X})<r$, then $\mathcal{P}(\cdot)$ is not unique because the SVD of $\bm{X}$ is not unique. Therefore, it is better to treat $\mathcal{P}(\bm{X})$ as a set-valued operator which outputs one representative from the set $\{\bm{U}\bm{V}^{\top}:\bm{X}=\bm{U}\bm{\Sigma}\bm{V}^{\top}\text{ is the SVD of }\bm{X}\}$.

The main purpose of introducing $\mathcal{P}(\cdot)$ is to correctly define the distance among $\bm{\Phi}$, $\bm{\Phi}^{(i)}$, and $\bm{Z}.$ From now on, we let

where the explicit forms of $\bm{Q}$, $\bm{Q}^{(i)}$, and $\bm{S}^{(i)}$ are easy to derive.

Now we decompose $\|\bm{\Delta}_{i}^{\top}(\bm{\Phi}-\bm{Z}\bm{Q})\|$ into three terms and find an upper bound for each of them:

Now it suffices to control each $T_{k},1\leq k\leq 3.$

One may wonder what the main difference of the blockwise analysis of the eigenvectors from the aforementioned works is. The difference comes from the appearance of $T_{3}$, which does not show up for $d=1$. Therefore, the estimation of $\|\bm{Q}-\bm{Q}^{(i)}\bm{S}^{(i)}\|$ requires additional treatments for $d\geq 2.$

Now we present the final estimation of $\|\bm{\Phi}_{i}-[\bm{A}\bm{Z}\bm{Q}]_{i}\bm{\Lambda}^{-1}\|$ in (3.6), which is used to establish Theorem 3.1 and 3.3. Note that the proof of Theorem 5.1 naturally includes the estimation of (5.2). We leave the estimation of (5.2) in Section 5.2 and 5.3.

Under the assumption of Theorem 3.1 and 3.3, i.e.,

Here $C_{0}>0$ is an absolute large constant. Then with at least $1-O(n^{-1}d^{-1})$, it holds

uniformly for all $1\leq i\leq n$. Moreover, we have

By using the key supporting result Theorem 5.1 above as well as Lemma 5.2 and Claim 5.3 below, we provide the proof of Theorem 3.1 and 3.3.

For two $d\times d$ invertible matrices $\bm{X}$ and $\bm{Y}$, it holds that

where $\sigma_{\min}(\cdot)$ denotes the smallest singular value of a matrix.

We will provide a proof of Lemma 5.2 in the appendix which uses the Davis-Kahan theorem. In fact, we found a different proof of the same result in after we finished this manuscript.

With probability at least $1-O(n^{-1})$, we have the following results under the assumption of Theorem 3.1 and 3.3.

It holds that $\|\bm{\Delta}\|\leq n/2$ and $\bm{\Lambda}\succeq\frac{n}{2}\bm{I}_{n}$ for both (OD) and (PM).

Here $\eta$ is defined in (5.3) and (5.4) for both cases respectively.

The claim will be confirmed in Section 5.2 and 5.3 for two scenarios respectively. This claim ensures that $\bm{\Phi}_{i}$ is well approximated by $n\bm{Q}\bm{\Lambda}^{-1}$ since

where the two terms are bounded by Theorem 5.1 and Claim 5.3 respectively and $[\bm{A}\bm{Z}\bm{Q}]_{i}=(n\bm{I}_{d}+\bm{\Delta}_{i}^{\top}\bm{Z})\bm{Q}$ for all $1\leq i\leq n.$ Then applying Lemma 5.2 immediately gives the main results.

With the results available above, we are ready to provide an upper bound for $\mathcal{P}(\bm{\Phi}_{i})-\mathcal{P}(\bm{\Phi}_{j})$. In fact, as long as $\|\mathcal{P}(\bm{\Phi}_{i})-\mathcal{P}(\bm{\Phi}_{j})\|\lesssim\eta$, we have $\|\mathcal{P}(\bm{\Phi}_{i})\mathcal{P}(\bm{\Phi}_{j})^{\top}-\bm{I}_{d}\|\lesssim\eta,\forall i\neq j.$ This finishes the proof of Theorem 3.1. For Theorem 3.3, it requires one extra step to show that $\|\mathcal{P}(\bm{\Phi}_{i})\mathcal{P}(\bm{\Phi}_{j})^{\top}-\bm{I}_{d}\|<1/2$ holds so that the rounding procedure indeed produces the planted permutation matrices correctly.

It suffices to estimate $\bm{\Phi}_{i}-\bm{\Phi}_{j}$ and apply Lemma 5.2. Under Theorem 5.1 and Claim 5.3, we have

where the third inequality uses Claim 5.3 and $\max_{1\leq i\leq n}\|\bm{\Phi}_{i}\|\geq 1$. In order to apply Lemma 5.2, we need to show $\sigma_{\min}(\bm{\Phi}_{i})$ is away from 0.

Let $i^{\prime}$ be the index with the largest $\|\bm{\Phi}_{i^{\prime}}\|$. Now we will show that $\|\mathcal{P}(\bm{\Phi}_{i^{\prime}})-\mathcal{P}(\bm{\Phi}_{j})\|\lesssim\eta$ for all $1\leq j\leq n$ and then $\|\mathcal{P}(\bm{\Phi}_{i})-\mathcal{P}(\bm{\Phi}_{j})\|\lesssim\eta$ follows from triangle inequality. Note that

where $``\otimes"$ denotes the Kronecker product. Note that all the singular values of $\bm{\Phi}$ are $\sqrt{n}$ and thus it holds that

This gives $|1-\sigma_{\min}(\bm{\Phi}_{i^{\prime}})|<\eta\|\bm{\Phi}_{i^{\prime}}\|$ and $\left|~{}\|\bm{\Phi}_{i^{\prime}}\|-1~{}\right|\lesssim\eta\|\bm{\Phi}_{i^{\prime}}\|$ which imply $\max_{1\leq i\leq n}\|\bm{\Phi}_{i}\|\leq 1+O(\eta)=O(1)$. Therefore,

which means $\sigma_{\min}(\bm{\Phi}_{i^{\prime}})>0.$ Applying Lemma 5.2 gives

where $|1-\sigma_{\min}(\bm{\Phi}_{i^{\prime}})|\lesssim\eta.$ By triangle inequality, $\|\mathcal{P}(\bm{\Phi}_{i})-\mathcal{P}(\bm{\Phi}_{j})\|\lesssim\eta$ holds for all pairs of $i$ and $j.$ This finishes the proof of Theorem 3.1.

For (PM) with $p>C_{0}\sqrt{n^{-1}\log(nd)}$ and a sufficiently large constant $C_{0}$, it holds that

Then the diagonal entries of $\mathcal{P}(\bm{\Phi}_{i})\mathcal{P}(\bm{\Phi}_{j})^{\top}$ are its largest $d$ entries since all the diagonal entries are greater than 1/2 while the off-diagonal entries are smaller than 1/2 in magnitude. Thus the rounding procedure in Algorithm 2 recovers the underlying permutation matrix which is $\bm{I}_{d}$. ∎

2 O​(d)O𝑑\text{O}(d) synchronization under Gaussian noise

This section is devoted to the proof of Theorem 5.1 for (OD). We first introduce all the necessary ingredients of the proof as follows and leave their proofs later. Then by using these facts, we can prove Theorem 5.1 for (OD). The key is to obtain upper bounds for $T_{1}$, $T_{2}$, $T_{3}$, and $\|\bm{\Phi}-\bm{Z}\bm{Q}\|$. We provide the bounds for each term below.

Roadmap: The estimation (5.5) of $\|\bm{\Delta}\|$ simply follows from $\|\bm{W}\|\leq 3\sqrt{nd}$ with high probability for a symmetric Gaussian random matrix $\bm{W}$. In particular, we assume $\|\bm{\Delta}\|\leq n/2$ for $\sigma<c_{0}\sqrt{nd^{-1}}$ and a small constant $c_{0}$. The bound (5.6) for the top $d$ eigenvalues of $\bm{A}=\bm{Z}\bm{Z}^{\top}+\bm{\Delta}$ follows from Weyl’s theorem and (5.5) where the top $d$ eigenvalues of $\bm{Z}\bm{Z}^{\top}$ are $n.$ The inequalities (5.7) and (5.8) follow from Lemma 5.5 and 5.6 respectively by using Davis-Kahan theorem (Theorem 6.2). We prove (5.9) in Lemma 5.7, and Lemma 5.8 implies (5.10) and (5.11).

where $\max_{1\leq i\leq n}\|\bm{\Phi}_{i}\|\geq 1.$ The estimation is bounded by

where $\sigma<c_{0}\sqrt{nd^{-1}}$ for a small constant $c_{0}$. From (3.6), it holds that

where the first term above follows from $\|\bm{Z}\|=\sqrt{n}$, (5.7), and $\sigma<\sqrt{n/d}.$ ∎

Next we proceed to prove (5.7)-(5.11). Our analysis will frequently use the following important fact about Gaussian random matrix.

[62, Theorem 4.4.5] For any $\bm{X}\in\hbox{\msbm{R}}^{n\times m}$ random matrix whose entries are i.i.d. standard normal random variables. For any $t>0$, it holds

with probability at least $1-2\exp(-t^{2}).$

The next two lemmas provide the proof of (5.7) and (5.8).

If $\bm{\Phi}$ consists of the top $d$ eigenvectors of $\bm{A}$ with $\bm{\Phi}^{\top}\bm{\Phi}=n\bm{I}_{d}$, then

The same bound applies to $\bm{\Phi}^{(i)}$. Let $\bm{\Phi}^{(i)}$ be the eigenvectors associated to the top $d$ eigenvalues of $\bm{A}^{(i)}$ in (3.7) and then

We directly apply Davis-Kahan theorem by letting $\bm{X}=\bm{Z}\bm{Z}^{\top}$ and $\bm{X}_{\bm{E}}=\bm{A}=\bm{Z}\bm{Z}^{\top}+\bm{\Delta}$ in Theorem 6.2 with $\bm{\Delta}=\sigma\bm{W}$. Note that the $d$th largest eigenvalue of $\bm{A}$ is at least $n-\|\bm{\Delta}\|$ and the $(d+1)$th eigenvalue of $\bm{Z}\bm{Z}^{\top}$ is 0. Thus set $\delta$ in Theorem 6.2 as $n-\|\bm{\Delta}\|$ and it holds

where $\|\bm{W}\bm{Z}\|\leq\|\bm{W}\|\|\bm{Z}\|\lesssim\sqrt{nd}\cdot\sqrt{n}=n\sqrt{d}.$ For $(\bm{\Phi}^{(i)},\bm{Z})$, simply use $\bm{\Delta}^{(i)}=\sigma\bm{W}^{(i)}$ which is defined in (3.7) where $\bm{W}^{(i)}$ equals $\bm{W}$ except that its $i$th block row and column are zero. Using Theorem 6.2 again gives

Combining Lemma 6.3 together with the results above gives

where $\bm{Q}=\mathcal{P}(\bm{Z}^{\top}\bm{\Phi})\in\hbox{\msbm{R}}^{d\times d}$ is an orthogonal matrix. For the singular values of $\bm{\Phi}^{\top}\bm{Z}$, we have

Let $\bm{\Phi}$ and $\bm{\Phi}^{(i)}$ be the top $d$ eigenvectors of $\bm{A}$ and $\bm{A}^{(i)}$ in (3.7) with $\bm{\Phi}^{\top}\bm{\Phi}=(\bm{\Phi}^{(i)})^{\top}\bm{\Phi}^{(i)}=n\bm{I}_{d}$ respectively. Then

where $\bm{S}^{(i)}$ is defined in (5.1).

The $d$th largest eigenvalue of $\bm{A}$ is at least $n-\sigma\|\bm{W}\|$ and the $(d+1)$th largest eigenvalue of $\bm{A}^{(i)}$ is at most $\sigma\|\bm{W}\|.$ Thus we have $\delta\geq n-2\sigma\|\bm{W}\|$ and Theorem 6.2 gives

where $\bm{A}-\bm{A}^{(i)}=\sigma(\bm{W}-\bm{W}^{(i)})$ holds. Let $\bm{W}_{i}\in\hbox{\msbm{R}}^{nd\times d}$ be the $i$th block column of $\bm{W}$. We have

Note that $\|\bm{W}_{i}\|\lesssim\sqrt{nd}$ holds for $1\leq i\leq n$ with high probability. Each entry of $\bm{W}_{i}^{\top}\bm{\Phi}^{(i)}$ is an independent $\mathcal{N}(0,n)$ random variable. As a result, $\|\bm{W}_{i}^{\top}\bm{\Phi}^{(i)}\|\lesssim\sqrt{n}(\sqrt{d}+\sqrt{\log n})$ uniformly for all $1\leq i\leq n$ with high probability, following from Theorem 5.4.

Now consider the $i$th block of $\bm{\Phi}-\bm{\Phi}^{(i)}\bm{S}^{(i)}$ and we have

This implies that if $\sigma n^{-\frac{1}{2}}(\sqrt{d}+\sqrt{\log n})<c_{0}$ with a sufficiently small constant $c_{0}$, then

where $\max_{1\leq i\leq n}\|\bm{\Phi}_{i}\|\geq 1.$ In other words, it holds

For $\bm{Q},\bm{Q}^{(i)}$ and $\bm{S}^{(i)}$ defined in (5.1), we have

since $\mathcal{P}(\bm{Z}^{\top}\bm{\Phi}^{(i)}\bm{S}^{(i)})=\bm{Q}^{(i)}\bm{S}^{(i)}.$ Applying Lemma 5.2 gives

where $\sigma_{\min}^{-1}(\bm{Z}^{\top}\bm{\Phi})\lesssim n^{-1}$ follows from Lemma 5.5 and 5.6. ∎

Suppose a sequence of $n$ matrices $\{\bm{M}_{i}\}_{i=1}^{n}$ and $\bm{M}^{\top}=[\bm{M}_{1}^{\top},\cdots,\bm{M}_{n}^{\top}]$ which is independent of $\bm{\Delta}_{i}=\sigma\bm{W}_{i}$. Then

with probability at least $1-O(n^{-2}).$ In particular, the following inequalities hold with probability at least $1-O(n^{-1})$

Denote the SVD of $\bm{M}\in\hbox{\msbm{R}}^{nd\times d}$ as $\bm{M}=\bm{U}\bm{\Sigma}\bm{V}^{\top}$, where $\bm{U}\in\hbox{\msbm{R}}^{nd\times d}$ with $\bm{U}^{\top}\bm{U}=\bm{I}_{d}$, $\bm{\Sigma}\in\hbox{\msbm{R}}^{d\times d}$, and $\bm{V}\in\hbox{\msbm{R}}^{d\times d}$. Note that $\bm{\Delta}_{i}=\sigma\bm{W}_{i}$ where $\bm{W}_{i}$ is an $nd\times d$ Gaussian random matrix.

where $\|\bm{\Sigma}\|=\|\bm{M}\|.$ Since $\bm{U}^{\top}\bm{U}=\bm{I}_{d}$, then $\bm{W}_{i}^{\top}\bm{U}$ is an asymmetric $d\times d$ Gaussian random matrix. Theorem 5.4 guarantees that $\|\bm{W}_{i}^{\top}\bm{U}\|$ is bounded by $\|\bm{W}_{i}^{\top}\bm{U}\|\lesssim\sqrt{d}+\sqrt{\log n}$ with probability at least $1-O(n^{-2})$. As a result, we have $\|\bm{\Delta}_{i}^{\top}\bm{M}\|\lesssim\sigma\|\bm{M}\|(\sqrt{d}+\sqrt{\log n}).$

Now by letting $\bm{M}=\bm{\Phi}^{(i)}-\bm{Z}\bm{Q}^{(i)}$ or $\bm{Z}$ which is independent of $\bm{\Delta}_{i}$, we have

hold uniformly for all $1\leq i\leq n$ with probability at least $1-O(n^{-1}).$ ∎

3 Object matching under uniform random corruption

Before proceeding to the official proof, we first show that the top $d$ eigenvectors of $\bm{A}_{\bm{G}}=[\bm{G}_{ij}]_{1\leq i,j\leq n}$ are equal to those of $\bm{A}$ in (3.2) and (3.3). Note that for the spiked matrix model (3.3) for (PM), the noise matrix $\bm{\Delta}$ is not mean zero, i.e., $\operatorname{\hbox{\msbm{E}}}\bm{\Delta}=-p^{-1}(1-p)d^{-1}\bm{I}_{n}\otimes\bm{J}_{d}$ is a block-diagonal matrix and $\|\operatorname{\hbox{\msbm{E}}}\bm{\Delta}\|\leq p^{-1}(1-p).$

The matrices $[\bm{G}_{ij}]_{1\leq i,j\leq n}$ and $\bm{A}$ share the same top $d$ eigenvectors.

Without loss of generality, we assume $\bm{G}_{i}=\bm{I}_{d}$. Note that $\bm{A}=\bm{Z}\bm{Z}^{\top}+\bm{\Delta}$.

where $X_{ij}\sim$Bernoulli($p$) and $\bm{P}_{ij}$ is a random permutation matrix.

Note that $\bm{A}_{G}$ is a nonnegative matrix with its leading eigenvector $\bm{1}_{nd}$. Also $\bm{A}=\bm{Z}\bm{Z}^{\top}+\bm{\Delta}$ has $\bm{1}_{nd}$ as its leading eigenvector if $\|\bm{\Delta}\|<n.$ As a result, $\bm{J}_{nd}$ and $\bm{I}_{nd}$ do not change the top $d$ eigenvectors of $\bm{A}$ and $\bm{A}_{G}$. ∎

We follow a similar route of proof as presented in Section 5.2. Under assumption of Theorem 5.1, $p>C_{0}n^{-\frac{1}{2}}\sqrt{\log(nd)}$ for some large constant $C_{0}$, the following inequalities hold with probability at least $1-O(n^{-1})$,

Roadmap: Here (5.12) is given in Lemma 5.10 and (5.13) directly follows from Weyl’s inequality and $\|\bm{\Delta}\|\leq n/2$ for $p>C_{0}n^{-\frac{1}{2}}\sqrt{\log(nd)}$; Lemma 5.11 and 5.13 give (5.14) and (5.15) respectively; The estimation of (5.16) is provided in Lemma 5.14; and Corollary 5.15 implies both (5.17) and (5.18).

where $p^{-1}n^{-\frac{1}{2}}\sqrt{\log(nd))}<1.$ The estimations in (3.6) and (5.2) are bounded by

where the first term is bounded by using (5.14). ∎

The estimation of $\|\bm{\Delta}\|$ uses the matrix Bernstein inequality, see Theorem 6.4 in the appendix.

The operator norm of $\bm{\Delta}$ is bounded by

with probability least $1-O(n^{-1}d^{-1}).$ In particular, if $p\geq C_{0}n^{-\frac{1}{2}}\sqrt{\log(nd)}$ for some large constant $C_{0}$,

where $\|\operatorname{\hbox{\msbm{E}}}\bm{\Delta}\|=p^{-1}(1-p).$

Let $\bm{Z}_{ij}$ be an $nd\times nd$ matrix whose $(i,j)$- and $(j,i)$-block equal $\bm{\Delta}_{ij}$ and $\bm{\Delta}_{ji}$ respectively and all the other blocks are 0, i.e.,

is a symmetric matrix where $\{\bm{e}_{i}\}_{i=1}^{n}$ are the canonical basis in $\hbox{\msbm{R}}^{n}.$ Here

Let’s first compute its variance: for $i<j$, we have

By using the independence between $X_{ij}$ and $\bm{P}_{ij}$, the expectations of $\bm{\Delta}_{ij}\bm{\Delta}_{ij}^{\top}$ and $\bm{\Delta}_{ij}^{\top}\bm{\Delta}_{ij}$ are

where $\operatorname{\hbox{\msbm{E}}}\left[(\bm{P}_{ij}-d^{-1}\bm{J}_{d})(\bm{P}_{ij}-d^{-1}\bm{J}_{d})^{\top}\right]=\operatorname{\hbox{\msbm{E}}}\left[(\bm{P}_{ij}-d^{-1}\bm{J}_{d})^{\top}(\bm{P}_{ij}-d^{-1}\bm{J}_{d})\right]=\bm{I}_{d}-d^{-1}\bm{J}_{d}.$ As a result, it holds that

which implies $\|\sum_{i<j}\operatorname{\hbox{\msbm{E}}}\bm{Z}_{ij}\bm{Z}_{ij}^{\top}\|\leq p^{-2}(1-p^{2})n$.

Applying Bernstein’s inequality results in

with probability at least $1-O(n^{-1}d^{-1}).$ ∎

If $\bm{\Phi}$ consists of the top $d$ eigenvectors of $\bm{A}$ with $\bm{\Phi}^{\top}\bm{\Phi}=n\bm{I}_{d}$, then

holds under (5.12). The same bound applies to $\bm{\Phi}^{(i)}$. Let $\bm{\Phi}^{(i)}$ be the eigenvectors associated to the top $d$ eigenvalues of $\bm{A}^{(i)}$, and then

We directly apply Davis-Kahan theorem by letting $\bm{X}=\bm{Z}\bm{Z}^{\top}$ and $\bm{X}_{\bm{E}}=\bm{A}$ in Theorem 6.2. First we specify the spectral gap:

where $\|\bm{\Delta}\|\leq p^{-1}\sqrt{n\log(nd)}$ and $\|\bm{\Delta}\bm{Z}\|\leq\|\bm{\Delta}\|\|\bm{Z}\|=\sqrt{n}\|\bm{\Delta}\|.$ For $(\bm{\Phi}^{(i)},\bm{Z})$, simply using Theorem 6.2 again with $\bm{A}^{(i)}=\bm{Z}\bm{Z}^{\top}+\bm{\Delta}^{(i)}$ leads to

where $\|\bm{\Delta}^{(i)}\|\leq\|\bm{\Delta}\|\leq p^{-1}\sqrt{n\log(nd)}<n/4.$

Combining Lemma 6.3 with the bounds above gives

This provides a lower bound for the smallest singular value of $\bm{Z}^{\top}\bm{\Phi}$ and $\bm{Z}^{\top}\bm{\Phi}^{(i)}$ which follows from

Now we proceed to prove (5.15)-(5.18) which rely on the following lemma.

Let $\bm{M}\in\hbox{\msbm{R}}^{nd\times d}$ be a matrix with its $j$th block $\bm{M}_{j}$, independent of $\bm{\Delta}_{i}$. For each fixed $1\leq i\leq n$, it holds that

with probability at least $1-O(n^{-2}d^{-2}).$

with $\operatorname{\hbox{\msbm{E}}}\bm{\Delta}_{ij}=0$ for $i\neq j$ and $\bm{\Delta}_{ii}=-p^{-1}(1-p)d^{-1}\bm{J}_{d}.$ It holds that

Now we apply the Bernstein inequality (Theorem 6.4) to estimate the first term above which is a sum of mean zero independent random matrices.

We first compute $\operatorname{\hbox{\msbm{E}}}\sum_{j\neq i}(\bm{\Delta}_{ij}\bm{M}_{j})^{\top}(\bm{\Delta}_{ij}\bm{M}_{j})$. For each $j$, we have

because the cross terms of $\bm{\Delta}_{ij}^{\top}\bm{\Delta}_{ij}$ are of mean zero and $\operatorname{\hbox{\msbm{E}}}\bm{P}_{ij}=d^{-1}\bm{J}_{d}.$ Therefore,

For $\operatorname{\hbox{\msbm{E}}}\sum_{j\neq i}(\bm{\Delta}_{ij}\bm{M}_{j})(\bm{\Delta}_{ij}\bm{M}_{j})^{\top}$, we first compute $\operatorname{\hbox{\msbm{E}}}\bm{\Delta}_{ij}\bm{M}_{j}(\bm{\Delta}_{ij}\bm{M}_{j})^{\top}$:

and the variance of $\sum_{j\neq i}\bm{\Delta}_{ij}\bm{M}_{j}$ is bounded by

Each term $\bm{\Delta}_{ij}\bm{M}_{j}$ is bounded by

where $\|\bm{\Delta}_{ij}\|\leq p^{-1}\|(X_{ij}-p)\bm{I}_{d}+(1-X_{ij})\bm{P}_{ij}\|\leq 2p^{-1}$. Now applying Bernstein inequality gives

Let $\bm{\Phi}$ and $\bm{\Phi}^{(i)}$ be the top $d$ eigenvectors of $\bm{A}$ and $\bm{A}^{(i)}$ with $\bm{\Phi}^{\top}\bm{\Phi}=(\bm{\Phi}^{(i)})^{\top}\bm{\Phi}^{(i)}=n\bm{I}_{d}$ respectively. Then

The $d$th largest eigenvalue of $\bm{A}$ is at least $n-\|\bm{\Delta}\|$ and the $(d+1)$th largest eigenvalue of $\bm{A}^{(i)}$ is at most $\|\bm{\Delta}\|.$ Thus we have $\delta\geq n-2\|\bm{\Delta}\|$ and

where $\bm{A}-\bm{A}^{(i)}=\bm{\Delta}-\bm{\Delta}^{(i)}$ holds.

Then it holds with probability at least $1-O(n^{-1}d^{-2})$ that

where $\|\bm{\Delta}_{i}\|\lesssim p^{-1}\sqrt{n\log(nd)}$ uses Lemma 5.10 and $\|\bm{\Delta}_{i}^{\top}\bm{\Phi}^{(i)}\|\leq p^{-1}\sqrt{n\log(nd)}\max_{1\leq j\leq n}\|\bm{\Phi}_{j}^{(i)}\|$ follows from the independence between $\bm{\Delta}_{i}$ and $\bm{\Phi}^{(i)}$ and Lemma 5.12.

holds uniformly for all $1\leq i\leq n$ with probability at least $1-O(n^{-1}d^{-2})$. Next we will show that $\max_{1\leq j\leq n}\|\bm{\Phi}_{j}^{(i)}\|\lesssim\max_{1\leq i\leq n}\|\bm{\Phi}_{i}\|.$ Let $j^{\prime}$ be the index such that $\|\bm{\Phi}_{j^{\prime}}^{(i)}\|=\max_{1\leq j\leq n}\|\bm{\Phi}_{j}^{(i)}\|$, then applying triangle inequality gives

This implies that if $p>C_{0}n^{-\frac{1}{2}}\sqrt{\log(nd)}$ with a sufficiently large constant $C_{0}$, then

The lower bound for the smallest singular value of $\bm{\Phi}^{\top}\bm{\Phi}^{(i)}$ directly follows from