A polynomial-time iterative algorithm for random graph matching with non-vanishing correlation

The random graph matching problem is motivated by questions from various applied fields such as social network analysis [40, 41], computer vision [8, 4], computational biology [48, 49, 19] and natural language processing [30]. In biology, an important problem is to identify proteins with similar structures/functions across different species. Toward this goal, directly comparing amino acid sequences that constitute proteins is often complicated, since genetic mutations of the species can result in significant variations of such a sequence. However, despite these variations, proteins typically maintain similar functions within each species’ metabolism. In light of this, biologists employ graph-based representations, such as Protein-Protein Interaction (PPI) graphs, for each species. Under the assumption that the topological structures of PPI graphs are similar across species, researchers can then effectively match proteins with similar functions by taking advantage of such similarity (and by possibly taking advantage of some other domain knowledge). This approach turns out to be successful and offers a nuanced understanding of phenomena such as the evolution of protein complexes. We refer the reader to [19] for more information on the topic. In the domain of social networks, data privacy is a fundamental and complicated problem. One complication arises from the fact that a user may have accounts in multiple social network platforms, where the user shares similar content or engages in comparable activities. Thus, from such similarities in different platforms, it is possible to infer their identities by aligning the respective graphs representing user interactions and content consumption. That is to say, it is possible to use graph matching techniques to deanonymize private social networks using information gleaned from public social platforms. In this scope, well-known examples include the deanonymization of Netflix using data from IMDb [40] and the deanonymization of Twitter using Flickr [41]. Viewing this problem from the opposite perspective, a deeper understanding on the random graph matching problem may offer insights on better mechanism to protect data privacy. With aforementioned applications in mind, the graph matching problem has been extensively studied recently from a theoretical point of view.

It is natural that different applications have their own distinct features, which mathematically boils down to a careful choice of underlying random graph model suitable for the desired application. Similar to most previous works on the random graph matching problem, in this paper we consider the correlated Erdős–Rényi random graph model, which is possibly an over idealization for any realistic network but nevertheless offers a good playground to develop insights and methods for this problem in general. By the collective efforts as in [10, 9, 31, 53, 52, 27, 13, 14], it is fair to say that we have a fairly complete understanding on the information thresholds for the problem of correlation detection and vertex matching. In contrast, the understanding on the computational aspect is far from being complete, and in what follows we briefly review the progress on this front.

A huge amount of efforts have been devoted to developing efficient algorithms for graph matching [42, 54, 34, 33, 24, 47, 1, 18, 21, 22, 5, 11, 12, 39, 26, 35, 36, 37, 28, 29, 38]. Previous to this work, arguably the best result is the recent work [38] (see also [28, 29] for a remarkable result on partial recovery of similar flavor when the average degree is $O(1)$), where the authors substantially improved a groundbreaking work [36] and obtained a polynomial time algorithm which succeeds as long as the correlation is above the square root of the Otter’s constant (the Otter’s constant is around ${0.338}$). In terms of methods, the present work seems drastically different from these existing works; instead, this work is closely related to our previous work [17] on a polynomial-time iterative algorithm that matches two correlated Gaussian Wigner matrices. One is encouraged to see [38, Section 1.1] and [17, Section 1.1] for a more elaborated review on previous algorithms. In addition, one is encouraged to see [17, Section 1.2] for discussions on novel features of this iterative algorithm, especially in comparison with the message-passing algorithm [29, 43] and a recent greedy algorithm for aligning two independent graphs [15].

While the present work can be viewed as an extension of [17], we do feel that we have conquered substantial obstacles and have made substantial contributions to this topic, as we describe below.

• •

  While in [38] (see also [29]) polynomial-time matching algorithms were obtained when the correlation is above the threshold from the Otter’s constant, our work establishes a polynomial-time algorithm as long as the average degree grows polynomially and the correlation is non-vanishing. In addition, the power in the running time only tends to $\infty$ as the correlation tends to $0$ (for each fixed $\alpha<1$), and we are under the feeling that this is the best possible; such feeling is supported by a recent work on the complexity of low-degree polynomials for graph matching [16].

• •

  From a conceptual point of view, our work demonstrates the robustness of the iterative matching algorithm proposed in [17]. This type of “algorithmic universality” is closely related to the universality phenomenon in random matrix theory, which roughly speaking postulates that the particular distribution for entries of random matrices is often irrelevant for spectral properties under investigation. Our work also encourages future study on even more ambitious perspectives for robustness, for instance algorithms that are robust with assumptions on underlying random graph models. This is of major interest since realistic networks are usually captured better by more structured graph models, such as the random geometric graph model [50], the random growing graph model [44] and the stochastic block model [45].

• •

  In terms of techniques, our work employs a method which argues that Gaussian approximation is valid typically. There are a couple of major challenges : (1) Gaussian approximation is valid only in a rather weak sense that the Radon-Nikodym derivative is not too large; (2) we can not simply ignore the atypical cases when Gaussian approximation is invalid since we have to analyze the conditional behavior given all the information the algorithm has exploited so far. The latter raises a major obstacle in our proof for theoretical guarantee; see Section 3.1 for more detailed discussions on how this obstacle is addressed.

In addition, it is natural to suspect that our work brings us one step closer to understanding computational phase transitions for random graph matching problems as well as understanding algorithmic perspectives for other matching problems (see, e.g., [6]). We refer an interested reader to [17, Section 1.3] and omit further discussions here.

We record in this subsection some notation conventions, and we point out that a list of commonly used notations is included at the end of the paper for better reference.

Denote the identity matrix by $\mathrm{I}$ and the zero matrix by $\mathrm{O}$. For a $d\!*\!m$ matrix $\mathrm{A}$, we use $\mathrm{A}^{*}$ to denote its transpose, and let $\|\mathrm{A}\|_{\mathrm{HS}}$ denote the Hilbert-Schmidt norm of $\mathrm{A}$. For $1\leq s\leq\infty$, define the $s$-norm of $\mathrm{A}$ by $\|\mathrm{A}\|_{s}=\sup\{\|\mathrm{A}x^{*}\|_{s}:\|x\|_{s}=1\}$. Note that when $\mathrm{A}$ is a symmetric square matrix, we have for $\tfrac{1}{s}+\tfrac{1}{t}=1$

$$\|\mathrm{A}\|_{s}=\sup\{y\mathrm{A}x^{*}:\|x\|_{s}=1,\|y\|_{t}=1\}=\sup\{\|y\mathrm{A}\|_{t}:\|y\|_{t}=1\}=\|\mathrm{A}^{*}\|_{t}=\|\mathrm{A}\|_{t}\,.$$

We further denote the operator norm of $\mathrm{A}$ by $\|\mathrm{A}\|_{\mathrm{op}}=\|\mathrm{A}\|_{2}$. If $m=d$, denote $\mathrm{det(A)}$ and $\mathrm{tr(A)}$ the determinant and the trace of $\mathrm{A}$, respectively. For $d$-dimensional vectors $x,y$ and a $d\!*\!d$ symmetric matrix $\Sigma$, let $\langle x,y\rangle_{\Sigma}=x\Sigma y^{*}$ be the “inner product” of $x,y$ with respect to $\Sigma$, and we further denote $\|x\|_{\Sigma}^{2}=x\Sigma x^{*}$. For two vectors $\gamma,\mu\in\mathbb{R}^{d}$, we say $\gamma\geq\mu$ (or equivalently $\mu\leq\gamma$) if their entries satisfy $\gamma(i)\geq\mu(i)$ for all $1\leq i\leq d$. The indicator function of a set $A$ is denoted by $\mathbf{1}_{A}$.

Without further specification, all asymptotics are taken with respect to $n\to\infty$. We also use standard asymptotic notations: for two sequences $\{a_{n}\}$ and $\{b_{n}\}$, we write $a_{n}=O(b_{n})$ or $a_{n}\lesssim b_{n}$, if $|a_{n}|\leq C|b_{n}|$ for some absolute constant $C$ and all $n$. We write $a_{n}=\Omega(b_{n})$ or $a_{n}\gtrsim b_{n}$, if $b_{n}=O(a_{n})$; we write $a_{n}=\Theta(b_{n})$ or $a_{n}\asymp b_{n}$, if $a_{n}=O(b_{n})$ and $a_{n}=\Omega(b_{n})$; we write $a_{n}=o(b_{n})$ or $b_{n}=\omega(a_{n})$, if $\frac{a_{n}}{b_{n}}\to 0$ as $n\to\infty$. We write $a_{n}\sim b_{n}$ if $\frac{a_{n}}{b_{n}}\to 1$.

We denote by $\mathrm{Ber}(p)$ the Bernoulli distribution with parameter $p$, denote by $\mathrm{Bin}(n,p)$ the Binomial distribution with $n$ trials and success probability $p$, denote by $\mathcal{N}(\mu,\sigma^{2})$ the normal distribution with mean $\mu$ and variance $\sigma^{2}$, and denote by $\mathcal{N}(\mu,\Sigma)$ the multi-variate normal distribution with mean $\mu$ and covariance matrix $\Sigma$. We say $(X,Y)$ is a pair of correlated binomial random variables, denoted as $\mathrm{CorBin}(N,M,p;m,\rho)$ for $m\leq\min\{N,M\}$ and $\rho\in[0,1]$, if $(X,Y)\sim\big{(}\sum_{k=1}^{N}b_{k},\sum_{k=1}^{M}b^{\prime}_{k}\big{)}$ with $b_{k},b^{\prime}_{l}\sim\mathrm{Ber}(p)$ such that $\{b_{k},b^{\prime}_{k}\}$ are independent with $\{b_{1},\ldots,b_{N},b^{\prime}_{1},\ldots,b^{\prime}_{M}\}\setminus\{b_{k},b^{\prime}_{k}\}$ and the covariance between $b_{k}$ and $b^{\prime}_{k}$ is $\rho$ when $k\leq m$, and that $b_{k}$ is independent with $\{b_{1},\ldots,b_{N},b^{\prime}_{1},\ldots,b^{\prime}_{M}\}\setminus\{b_{k}\}$ and $b^{\prime}_{l}$ is independent with $\{b_{1},\ldots,b_{N},b^{\prime}_{1},\ldots,b^{\prime}_{M}\}\setminus\{b^{\prime}_{l}\}$ when $k,l>m$. We say $X$ is a sub-Gaussian variable, if there exists a positive constant $C$ such that $\mathbb{P}(|X|\geq t)\leq 2e^{-{t^{2}}/{C^{2}}}$, and we use $\|X\|_{\psi_{2}}=\inf\big{\{}C>0:\mathbb{E}[\exp\{\frac{X^{2}}{C^{2}}\}]\leq 2\big{\}}$ to denote its sub-Gaussian norm.

Similar to [17], we make appropriate preprocessing on random graphs such that we only need to consider graphs with directed edges. We first make a technical assumption that we only need to consider the case when $\rho$ is a sufficiently small constant, which can be easily achieved by keeping each edge independently with a sufficiently small constant probability.

Now, we define $\overrightarrow{G}$ from $G$. For any $u\neq v\in V$, we do the following:

• •

  if $(u,v)\in E(G)$, then independently among all such $(u,v)$:

  	$\displaystyle\mbox{with probability }\frac{1}{2}-\frac{q}{4},\mbox{ set }\overrightarrow{(u,v)}\in\overrightarrow{G},\overrightarrow{(v,u)}\not\in\overrightarrow{G}\,,$
  	$\displaystyle\mbox{with probability }\frac{1}{2}-\frac{q}{4},\mbox{ set }\overrightarrow{(u,v)}\not\in\overrightarrow{G},\overrightarrow{(v,u)}\in\overrightarrow{G}\,,$
  	$\displaystyle\mbox{with probability }\frac{q}{4},\mbox{ set }\overrightarrow{(u,v)}\in\overrightarrow{G},\overrightarrow{(v,u)}\in\overrightarrow{G}\,,$
  	$\displaystyle\mbox{with probability }\frac{q}{4},\mbox{ set }\overrightarrow{(u,v)}\not\in\overrightarrow{G},\overrightarrow{(v,u)}\not\in\overrightarrow{G}\,;$

• •

  if $(u,v)\not\in E(G)$, then set $\overrightarrow{(u,v)}\not\in\overrightarrow{G},\overrightarrow{(v,u)}\not\in\overrightarrow{G}$.

We define $\overrightarrow{\mathsf{G}}$ from $\mathsf{G}$ in the same manner such that $\overrightarrow{G}$ and $\overrightarrow{\mathsf{G}}$ are conditionally independent given $(G,\mathsf{G})$. We continue to use the convention that $\overrightarrow{G}_{u,v}=\mathbf{1}_{\{\overrightarrow{(u,v)}\in\overrightarrow{G}\}}$. It is then straightforward to verify that $\{\overrightarrow{G}_{u,v}:u\neq v\}$ and $\{\mathsf{\overrightarrow{\mathsf{G}}_{u,v}:u\neq v}\}$ are two families of i.i.d. Bernoulli random variables with parameter $\frac{q}{2}$. In addition, we have

$\displaystyle\mathbb{E}[\overrightarrow{G}_{u,v}{\overrightarrow{\mathsf{G}}_{\pi(u),\pi(v)}}]=\mathbb{E}[\overrightarrow{G}_{u,v}{\overrightarrow{\mathsf{G}}_{\pi(v),\pi(u)}}]=\frac{q(q+\rho(1-q))}{4}\,.$

Thus, $\overrightarrow{G},\overrightarrow{\mathsf{G}}$ are edge-correlated directed graphs, denoted as $\overrightarrow{\mathcal{G}}(n,\hat{q},\hat{\rho})$, such that $\hat{q}=\frac{q}{2}$ and $\hat{\rho}=\frac{1-q}{2-q}\rho$. Also note that $\hat{q}\geq n^{-\alpha+o(1)}$ and $\hat{\rho}\in[\frac{\rho}{3},\frac{\rho}{2})$ since $q\leq 1/2$. From now on we will work on the directed graph $(\overrightarrow{G},\overrightarrow{\mathsf{G}})$.

For a pair of standard bivariate normal variables $(X,Y)$ with correlation $u$, we define $\phi:[-1,1]\mapsto[0,1]$ by (below the number 10 is somewhat arbitrarily chosen)

$\displaystyle\phi(u)=\mathbb{P}(|X|\geq 10,|Y|\geq 10)\,.$

(2.1)

In addition, we define

$${}\iota_{\mathrm{ub}}=\sup_{x\in(0,1]}\Big{\{}\frac{\phi(x)-\phi(0)}{x^{2}}\Big{\}}\mbox{ and }\iota_{\mathrm{lb}}=\inf_{x\in(0,1]}\Big{\{}\frac{\phi(x)-\phi(0)}{x^{2}}\Big{\}}\,.$$

(2.2)

From the definition we know $\phi$ is strictly increasing and by [17, Claims 2.6 and 2.8] we have $\phi^{\prime}(0)=0,\phi^{\prime\prime}(0)>0$, and thus both $\iota_{\mathrm{ub}}$ and $\iota_{\mathrm{lb}}$ are positive and finite. Also we write $\mathfrak{a}=\phi(1)=\mathbb{P}(|X|\geq 10)$. Recalling in Subsection 2.1 it was shown that $\rho$ can be assumed to be a sufficiently small constant, from now on we will assume that

$${}\hat{\rho}\leq\rho\leq\min\big{\{}\mathfrak{a}-\mathfrak{a}^{2},\iota_{\mathrm{ub}}^{-1},\tfrac{1}{10}\big{\}}\,.$$

(2.3)

Let $\kappa=\kappa(\hat{\rho})$ be a sufficiently large constant depending on $\hat{\rho}$ whose exact value will be decided later in (2.10). Set $K_{0}=\kappa$. We then arbitrarily choose a sequence $A=(u_{1},u_{2},\ldots,u_{K_{0}})$ where $u_{i}$’s are distinct vertices in $V$, and we list all the sequences of length $K_{0}$ with distinct elements in $\mathsf{V}$ as $\mathsf{A}_{1},\mathsf{A}_{2},\ldots,\mathsf{A}_{\mathtt{M}}$ where $\mathtt{M}=\mathtt{M}(n,\hat{\rho},\hat{q})=\prod_{i=0}^{K_{0}-1}(n-i)$. As in [17], for each $1\leq\mathtt{m}\leq\mathtt{M}$, we will run a procedure of initialization and iteration and clearly for one of them (although a priori we do not know which one it is) we are running the algorithm as if we have $K_{0}$ true pairs as seeds. For convenience, when describing the initialization and the iteration we will drop $\mathtt{m}$ from notations, but we emphasize that this procedure will be applied to each $\mathsf{A}_{\mathtt{m}}$. Having this clarified, we take a fixed $\mathtt{m}$ and denote $\mathsf{A}_{\mathtt{m}}=(\mathsf{u}_{1},\ldots,\mathsf{u}_{K_{0}})$. In what follows, we abuse the notation and write $V\setminus A$ when regarding $A$ as a set (similarly for $\mathsf{A}_{\mathtt{m}}$).

We next describe our initialization procedure. As discussed earlier, to the contrary of the case for Wigner matrices, we have to investigate the neighborhood around a seed up to a certain depth in order to get information for a large number of vertices. To this end, we choose an integer $\chi\leq\frac{1}{1-\alpha}$ as the depth such that

$$(n\hat{q})^{\chi}=o\big{(}ne^{-(\log\log n)^{100}}\big{)}\mbox{ and }(n\hat{q})^{\chi+1}=\Omega\big{(}ne^{-(\log\log n)^{100}}\big{)}$$

(2.4)

(this is possible since $\alpha<1$). We choose such $\chi$ since on the one hand we wish to see a large number of vertices near the seed and on the other hand we want to reveal only a vanishing fraction of the edges. Now for $1\leq k\leq K_{0}$, define the seeds

$\displaystyle\aleph^{(0)}_{k}=\{u_{k}\},\Upsilon^{(0)}_{k}=\{\mathsf{u}_{k}\}\,,\mbox{ and }\vartheta_{0}=\varsigma_{0}=1/n\,.$

(2.5)

Then for $1\leq a\leq\chi$, we iteratively define the $a$-neighborhood of each seed by

		$\displaystyle\aleph^{(a)}_{k}=\Big{\{}v\in V\setminus\big{(}\cup_{1\leq k\leq K_{0},0\leq j\leq a-1}\aleph^{(j)}_{k}\big{)}:\overrightarrow{G}_{v,u}=1\mbox{ for some }u\in\aleph^{(a-1)}_{k}\Big{\}}\,,$		(2.6)
		$\displaystyle\Upsilon^{(a)}_{k}=\Big{\{}\mathsf{v}\in\mathsf{V}\setminus\big{(}\cup_{1\leq k\leq K_{0},0\leq j\leq a-1}\Upsilon^{(j)}_{k}\big{)}:\overrightarrow{\mathsf{G}}_{\mathsf{v},\mathsf{u}}=1\mbox{ for some }\mathsf{u}\in\Upsilon^{(a-1)}_{k}\Big{\}}\,.$

Also, for $1\leq a\leq\chi$ we iteratively define

		$\displaystyle\vartheta_{a}=\mathbb{P}(X\geq 1)\mbox{ where }X\sim\mathrm{Bin}(\vartheta_{a-1}n,\hat{q})\,,$		(2.7)
		$\displaystyle\varsigma_{a}=\mathbb{P}(X\geq 1,Y\geq 1)\mbox{ where }(X,Y)\sim\mathrm{CorBin}(\vartheta_{a-1}n,\vartheta_{a-1}n,\hat{q};\varsigma_{a-1}n,\hat{\rho})\,.$

We will show in Subsection 3.3 that actually we have

$$|\aleph^{(a)}_{k}|/n,|\Upsilon^{(a)}_{k}|/n\approx\vartheta_{a}\mbox{ and }|\aleph^{(a)}_{k}|\cap\Upsilon^{(a)}_{k}|/n\approx\varsigma_{a}\,.$$

Let $\mathtt{d}_{\chi}=\mathtt{d}_{\chi}(n,\hat{q})$ be the minimal integer such that $\mathbb{P}(\mathrm{Bin}(n\vartheta_{\chi},\hat{q})\geq\mathtt{d}_{\chi})<\frac{1}{2}$, and set

		$\displaystyle\Gamma^{(0)}_{k}=\aleph^{(\chi+1)}_{k}=\Big{\{}v\in V\setminus\big{(}\cup_{1\leq k\leq K_{0},0\leq j\leq\chi}\aleph^{(j)}_{k}\big{)}:\sum_{u\in\aleph^{(\chi)}_{k}}\overrightarrow{G}_{v,u}\geq\mathtt{d}_{\chi}\Big{\}}\,,$		(2.8)
		$\displaystyle\Pi^{(0)}_{k}=\Upsilon^{(\chi+1)}_{k}=\Big{\{}\mathsf{v}\in\mathsf{V}\setminus\big{(}\cup_{1\leq k\leq K_{0},0\leq j\leq\chi}\Upsilon^{(j)}_{k}\big{)}:\sum_{\mathsf{u}\in\Upsilon^{(\chi)}_{k}}\overrightarrow{\mathsf{G}}_{\mathsf{v},\mathsf{u}}\geq\mathtt{d}_{\chi}\Big{\}}\,.$

And we further define

		$\displaystyle\vartheta=\vartheta_{\chi+1}=\mathbb{P}(X\geq\mathtt{d}_{\chi})\mbox{ where }X\sim\mathrm{Bin}(\vartheta_{\chi}n,\hat{q})\,,$		(2.9)
		$\displaystyle\varsigma=\varsigma_{\chi+1}=\mathbb{P}(X,Y\geq\mathtt{d}_{\chi})\mbox{ where }X,Y\sim\mathrm{CorBin}(\vartheta_{\chi}n,\vartheta_{\chi}n,\hat{q};\varsigma_{\chi}n,\hat{\rho})\,.$

We may then choose $K_{0}=\kappa$ sufficiently large such that $K_{0}\geq\frac{10^{34}\iota_{\mathrm{ub}}^{2}\hat{\rho}^{-20}}{\iota_{\mathrm{lb}}^{2}(\mathfrak{a}-\mathfrak{a}^{2})^{2}}$ and

$$\frac{\log(K_{0}\iota_{\mathrm{lb}}^{2}(\mathfrak{a}-\mathfrak{a}^{2})^{2}\hat{\rho}^{20}/10^{30}\iota_{\mathrm{ub}}^{2})}{\log({K_{0}\iota_{\mathrm{lb}}^{4}\hat{\rho}^{24}\varepsilon_{0}^{2}}/16\cdot 10^{30}\iota_{\mathrm{ub}}^{2})}\leq 1.01\mbox{ where }\varepsilon_{0}=\frac{\varsigma-\vartheta^{2}}{2(\vartheta-\vartheta^{2})}\,.$$

(2.10)

In addition, we define $\Phi^{(0)},\Psi^{(0)}$ to be $K_{0}\!*\!K_{0}$ matrices by

$\displaystyle\Phi^{(0)}=\mathrm{I}\mbox{ and }\Psi^{(0)}=\frac{\varsigma-\vartheta^{2}}{2(\vartheta-\vartheta^{2})}\mathrm{I}\,,$

(2.11)

and in the iterative steps we will also construct $K_{t},\varepsilon_{t},\Gamma^{(t)}_{k},\Pi^{(t)}_{k}$ and $\Phi^{(t)},\Psi^{(t)}$ for $t\geq 1$. Similarly as in [17], the matrices $\Phi^{(t)}$ and $\Psi^{(t)}$ are supposed to approximate the cardinalities of the sets $\big{\{}\Gamma^{(t)}_{k},\Pi^{(t)}_{k}\big{\}}$ in the following sense. Write

$${}\mathfrak{a}_{t}=\begin{cases}\mathfrak{a},&t\geq 1\,;\\ \vartheta,&t=0\,.\end{cases}$$

(2.12)

Then somewhat informally, we expect that

	$\displaystyle\frac{1}{n}|\Gamma^{(t)}_{i}|,\frac{1}{n}|\Pi^{(t)}_{i}|\approx\mathfrak{a}_{t}\,,$		(2.13)
	$\displaystyle\frac{\frac{1}{n}|\Gamma^{(t)}_{i}\cap\Gamma^{(t)}_{j}|-\frac{1}{n}|\Gamma^{(t)}_{i}|-\frac{1}{n}|\Gamma^{(t)}_{j}|+\mathfrak{a}_{t}^{2}}{\mathfrak{a}_{t}-\mathfrak{a}_{t}^{2}}\approx\frac{\frac{1}{n}|\Gamma^{(t)}_{i}\cap\Gamma^{(t)}_{j}|-\mathfrak{a}_{t}^{2}}{\mathfrak{a}_{t}-\mathfrak{a}_{t}^{2}}\approx\Phi^{(t)}(i,j)\,,$		(2.14)
	$\displaystyle\frac{\frac{1}{n}|\Pi^{(t)}_{i}\cap\Pi^{(t)}_{j}|-\frac{1}{n}|\Pi^{(t)}_{i}|-\frac{1}{n}|\Pi^{(t)}_{j}|+\mathfrak{a}_{t}^{2}}{\mathfrak{a}_{t}-\mathfrak{a}_{t}^{2}}\approx\frac{\frac{1}{n}|\Pi^{(t)}_{i}\cap\Pi^{(t)}_{j}|-\mathfrak{a}_{t}^{2}}{\mathfrak{a}_{t}-\mathfrak{a}_{t}^{2}}\approx\Phi^{(t)}(i,j)\,,$		(2.15)
	$\displaystyle\frac{\frac{1}{n}|\Gamma^{(t)}_{i}\cap\Pi^{(t)}_{j}|-\frac{1}{n}|\Gamma^{(t)}_{i}|-\frac{1}{n}|\Pi^{(t)}_{j}|+\mathfrak{a}_{t}^{2}}{\mathfrak{a}_{t}-\mathfrak{a}_{t}^{2}}\approx\frac{\frac{1}{n}|\Gamma^{(t)}_{i}\cap\Pi^{(t)}_{j}|-\mathfrak{a}_{t}^{2}}{\mathfrak{a}_{t}-\mathfrak{a}_{t}^{2}}\approx\Psi^{(t)}(i,j)\,.$		(2.16)

As in [17, Lemma 2.1], in order to facilitate our analysis later we will also need an important property on the eigenvalues of $\Phi^{(t)}$ and $\Psi^{(t)}$:

		$\displaystyle\Phi^{(t)}\mbox{ has at least }\frac{3}{4}K_{t}\mbox{ eigenvalues between $0.9$ and $1.1$}\,,$		(2.17)
	and	$\displaystyle\Psi^{(t)}\mbox{ has at least }\frac{3}{4}K_{t}\mbox{ eigenvalues between $0.9\varepsilon_{t}$ and $1.1\varepsilon_{t}$}\,.$		(2.18)

We will show in Subsection 3.3 that (2.13)–(2.18) are satisfied for $t=0$. The main challenge is to construct $(\Gamma_{k}^{(t+1)},\Pi_{k}^{(t+1)})$ and $\Phi^{(t+1)},\Psi^{(t+1)}$ such that (2.13)–(2.18) hold for $t+1$, under the inductive assumption that (2.13)–(2.18) hold for $t$. We conclude this subsection with some bounds on $(\vartheta_{k},\varsigma_{k})$.

We reiterate that in this subsection we are describing the iteration for a fixed $1\leq\mathtt{m}\leq\mathtt{M}$ and eventually this iterative procedure will be applied to each $\mathtt{m}$. Define

$\displaystyle K_{t+1}=\frac{1}{\varkappa}K_{t}^{2}\mbox{ where }\varkappa=\varkappa(\hat{\rho})=\frac{10^{30}\iota_{\mathrm{ub}}^{2}\hat{\rho}^{-20}}{\iota_{\mathrm{lb}}^{2}(\mathfrak{a}-\mathfrak{a}^{2})^{2}}$

(2.19)

for $t\geq 0$. Since we have assumed $K_{0}\geq 10^{4}\varkappa$, we can then prove by induction that

$${}10^{30}\hat{\rho}^{20}(\mathfrak{a}-\mathfrak{a}^{2})^{2}K_{t}^{2}\geq K_{t+1}\geq 10^{4}K_{t}\,.$$

(2.20)

We now suppose that $(\Gamma^{(s)}_{k},\Pi^{(s)}_{k})_{1\leq k\leq K_{s}}$ and $\Phi^{(s)},\Psi^{(s)}$ have been constructed for $s\leq t$ (which will be implemented inductively via (2.29) as we describe next). Recall that we are working under the assumption that (2.13)–(2.18) hold for $s\leq t$. For $v\in V$ and $\mathsf{v}\in\mathsf{V}$, define $D^{(t)}_{v},\mathsf{D}^{(t)}_{\mathsf{v}}\in\mathbb{R}^{K_{t}}$ to be the “normalized degrees” of $v$ in $\Gamma^{(t)}_{k}$ and of $\mathsf{v}$ in $\Pi^{(t)}_{k}$ as follows:

	$\displaystyle D^{(t)}_{v}(k)$	$\displaystyle=\frac{1}{\sqrt{(\mathfrak{a}_{t}-\mathfrak{a}_{t}^{2})n\hat{q}(1-\hat{q})}}\sum_{u\in V}(\mathbf{1}_{u\in\Gamma^{(t)}_{k}}-\mathfrak{a}_{t})(\overrightarrow{G}_{v,u}-\hat{q})\,,$		(2.21)
	$\displaystyle\mathsf{D}^{(t)}_{\mathsf{v}}(k)$	$\displaystyle=\frac{1}{\sqrt{(\mathfrak{a}_{t}-\mathfrak{a}_{t}^{2})n\hat{q}(1-\hat{q})}}\sum_{\mathsf{u}\in\mathsf{V}}(\mathbf{1}_{\mathsf{u}\in\Pi^{(t)}_{k}}-\mathfrak{a}_{t})(\overrightarrow{\mathsf{G}}_{\mathsf{v,u}}-\hat{q})\,.$

Recalling (2.12), we note that there is a difference between the definition (2.21) for $t=0$ and $t\geq 1$; this is because $\Gamma^{(0)}_{k}$ and $\Pi^{(0)}_{k}$ may only contain a vanishing fraction of vertices. We also point out that similar to [17], in the above definition we used the “centered” version of $\mathbf{1}_{u\in\Gamma^{(t)}_{k}}$ and $\mathbf{1}_{\mathsf{u}\in\Pi^{(t)}_{k}}$ since (2.13) suggests that intuitively each vertex $u$ (respectively, $\mathsf{u}$) has probability approximately $\mathfrak{a}_{t}$ to belong to $\Gamma^{(t)}_{k}$ (respectively, $\Pi^{(t)}_{k}$); such centering will be useful for our proof later as it leads to additional cancellation.

Assuming Lemma 2.2, we can then write $\Phi^{(t)}$ and $\Psi^{(t)}$ as their spectral decompositions:

$$\Phi^{(t)}=\sum^{K_{t}}_{i=1}\lambda^{(t)}_{i}\left({\nu^{(t)}_{i}}\right)^{*}\left(\nu^{(t)}_{i}\right)\mbox{ and }\Psi^{(t)}=\sum_{i=1}^{K_{t}}\mu^{(t)}_{i}\left({\xi^{(t)}_{i}}\right)^{*}\left(\xi^{(t)}_{i}\right)$$

(2.22)

where

$$\lambda^{(t)}_{i}\in(0.9,1.1),\mu^{(t)}_{i}\in(0.9\varepsilon_{t},1.1\varepsilon_{t})\mbox{ for }1\leq i\leq\frac{3K_{t}}{4}$$

(2.23)

and $\nu_{i}^{(t)},\xi_{i}^{(t)}$ are the unit eigenvectors with respect to $\lambda_{i}^{(t)},\mu_{i}^{(t)}$ respectively. Next, for $s,t$ we define $\mathrm{M}_{\Gamma}^{(t,s)},\mathrm{M}_{\Pi}^{(t,s)},\mathrm{P}_{\Gamma,\Pi}^{(t,s)}$ to be $K_{t}\!*\!K_{s}$ matrices as follows:

	$\displaystyle\mathrm{M}_{\Gamma}^{(t,s)}(i,j)$	$\displaystyle=\frac{|\Gamma^{(t)}_{i}\cap\Gamma^{(s)}_{j}|-\mathfrak{a}_{s}|\Gamma^{(t)}_{i}|-\mathfrak{a}_{t}|\Gamma^{(s)}_{j}|+\mathfrak{a}_{s}\mathfrak{a}_{t}n}{\sqrt{(\mathfrak{a}_{s}-\mathfrak{a}_{s}^{2})(\mathfrak{a}_{t}-\mathfrak{a}_{t}^{2})}n}\,,$		(2.24)
	$\displaystyle\mathrm{M}_{\Pi}^{(t,s)}(i,j)$	$\displaystyle=\frac{|\Pi^{(t)}_{i}\cap\Pi^{(s)}_{j}|-\mathfrak{a}_{s}|\Pi^{(t)}_{i}|-\mathfrak{a}_{t}|\Pi^{(s)}_{j}|+\mathfrak{a}_{s}\mathfrak{a}_{t}n}{\sqrt{(\mathfrak{a}_{s}-\mathfrak{a}_{s}^{2})(\mathfrak{a}_{t}-\mathfrak{a}_{t}^{2})}n}\,,$
	$\displaystyle\mathrm{P}_{\Gamma,\Pi}^{(t,s)}(i,j)$	$\displaystyle=\frac{|\pi(\Gamma^{(t)}_{i})\cap\Pi^{(s)}_{j}|-\mathfrak{a}_{s}|\Gamma^{(t)}_{i}|-\mathfrak{a}_{t}|\Pi^{(s)}_{j}|+\mathfrak{a}_{s}\mathfrak{a}_{t}n}{\sqrt{(\mathfrak{a}_{s}-\mathfrak{a}_{s}^{2})(\mathfrak{a}_{t}-\mathfrak{a}_{t}^{2})}n}\,.$

These matrices actually represent the covariance matrices for random vectors of the form $D^{(t)}_{v}$ and $\mathsf{D}^{(s)}_{\pi(v)}$. To get a rough intuition of this, we (formally incorrectly) regard $D_{v}^{(s)}$ as a linear combination of $\{G_{u,v}\}$ with deterministic coefficients (and the same applies to $\mathsf{D}^{(s)}_{\mathsf{v}}$). Then we can see that for all $v\in V$, the “correlation” between $D^{(t)}_{v}(i)$ and $D^{(s)}_{v}(j)$ equals

$$\frac{1}{\sqrt{(\mathfrak{a}_{t}-\mathfrak{a}_{t}^{2})(\mathfrak{a}_{s}-\mathfrak{a}_{s}^{2})}n}\sum_{u\in V\setminus A}(\mathbf{1}_{u\in\Gamma^{(t)}_{i}}-\mathfrak{a}_{t})(\mathbf{1}_{u\in\Gamma^{(s)}_{j}}-\mathfrak{a}_{s})=\mathrm{M}_{\Gamma}^{(t,s)}(i,j)\,.$$

This justifies our definition of $\mathrm{M}_{\Gamma}^{(t,s)}$ which aims to record the correlation between $D^{(t)}_{v}$ and $D^{(s)}_{v}$. Similarly under the same simplification we have $\mathrm{M}_{\Pi}^{(t,s)}$ (respectively, $\hat{\rho}\mathrm{P}_{\Gamma,\Pi}^{(t,s)}$) is the correlation matrix between $\mathsf{D}^{(t)}_{\mathsf{v}}$ and $\mathsf{D}^{(s)}_{\mathsf{v}}$ (respectively, between $D^{(t)}_{v}$ and $\mathsf{D}^{(s)}_{\pi(v)}$). In addition, from (2.13)–(2.16) we expect that $\mathrm{M}_{\Gamma}^{(t,t)},\mathrm{M}_{\Pi}^{(t,t)}\approx\Phi^{(t)}$ and $\mathrm{P}_{\Gamma,\Pi}^{(t,t)}\approx\Psi^{(t)}$. Note that $\mathrm{M}_{\Gamma},\mathrm{M}_{\Pi}$ are accessible by the algorithm but $\mathrm{P}_{\Gamma,\Pi}$ is not (since it relies on the latent matching). We further define two linear subspaces as follows:

	$\displaystyle\mathrm{W}^{(t)}$	$\displaystyle\overset{\triangle}{=}\big{\{}x\in\mathbb{R}^{K_{t}}:x\mathrm{M}_{\Gamma}^{(t,s)}=0,x\mathrm{M}_{\Pi}^{(t,s)}=0,\mbox{ for all }s<t\big{\}}\,,$		(2.25)
	$\displaystyle\mathrm{V}^{(t)}$	$\displaystyle\overset{\triangle}{=}\mathrm{span}\big{\{}\nu^{(t)}_{1},\nu^{(t)}_{2},\ldots,\nu^{(t)}_{\frac{3}{4}K_{t}}\big{\}}\cap\mathrm{span}\big{\{}\xi^{(t)}_{1},\xi^{(t)}_{2},\ldots,\xi^{(t)}_{\frac{3}{4}K_{t}}\big{\}}\cap\mathrm{W}^{(t)}\,.$

We refer to [17, Remark 3.3] for underlying reasons of the definition above. Note that the number of linear constraints posed on $\mathrm{W}^{(t)}$ are at most $2\sum_{i=1}^{t}K_{i-1}$. So

$$\dim(\mathrm{V}^{(t)})\geq\frac{3}{4}K_{t}+\frac{3}{4}K_{t}+\dim(\mathrm{W}_{t})-2K_{t}\geq\frac{1}{2}K_{t}-2\sum_{i=1}^{t}K_{i-1}\overset{\eqref{eq-increasing-dimension}}{\geq}0.49K_{t}\,.$$

As proved in [17, (2.10) and (2.11)], we can choose $\eta^{(t)}_{1},\eta^{(t)}_{2},\ldots,\eta^{(t)}_{\frac{1}{12}K_{t}}$ from $\mathrm{V}^{(t)}$ such that

	$\displaystyle\eta^{(t)}_{i}\mathrm{M}_{\Gamma}^{(t,t)}\big{(}\eta^{(t)}_{j}\big{)}^{*}=\eta^{(t)}_{i}\mathrm{M}_{\Pi}^{(t,t)}\big{(}\eta^{(t)}_{j}\big{)}^{*}=\eta^{(t)}_{i}\Psi^{(t)}\big{(}\eta^{(t)}_{j}\big{)}^{*}=0\,,$		(2.26)
	$\displaystyle\eta^{(t)}_{i}\Phi^{(t)}\big{(}\eta^{(t)}_{i}\big{)}^{*}=1,\quad 2\varepsilon_{t}\geq\eta^{(t)}_{i}\Psi^{(t)}\big{(}\eta^{(t)}_{i}\big{)}^{*}\geq 0.5\varepsilon_{t}\,.$		(2.27)

Furthermore, we must have $\big{\|}\eta^{(t)}_{i}\big{\|}^{2}\in(\frac{1}{2},2)$. As in [17], we will project the degrees $D^{(t)}_{v},\mathsf{D}^{(t)}_{\mathsf{v}}$ to a set of carefully chosen directions in the space spanned by all $\eta_{i}$’s. These directions are defined as follows: we sample $\beta^{(t)}_{k}(j)$ as i.i.d. uniform variables on $\{-1,1\}$. By [17, Proposition 2.4], these $\beta^{(t)}_{k}(j)$’s satisfy [17, (2.21)–(2.24)] with probability at least 0.5. As in [17], we will keep resampling until these requirements are satisfied. Define

$\displaystyle\sigma_{k}^{(t)}=\sqrt{\frac{12}{K_{t}}}\sum_{j=1}^{\frac{1}{12}K_{t}}\beta^{(t)}_{k}(j)\eta_{j}^{(t)}\mbox{ for }k=1,2,\ldots,K_{t+1}\,.$

(2.28)

We sample i.i.d. standard normal variables $\{W^{(t)}_{v}(i),\mathsf{W}^{(t)}_{\mathsf{v}}(i):1\leq i\leq\frac{K_{t}}{12}\}$ and complete our iteration by setting

		$\displaystyle\Gamma^{(t+1)}_{k}=\Big{\{}v\in V:\frac{1}{\sqrt{2}}\Big{|}\sqrt{\frac{12}{K_{t}}}\langle\beta^{(t)}_{k},W^{(t)}_{v}\rangle+\langle\sigma^{(t)}_{k},D^{(t)}_{v}\rangle\Big{|}\geq 10\Big{\}}\,,$		(2.29)
		$\displaystyle\Pi^{(t+1)}_{k}=\Big{\{}\mathsf{v}\in\mathsf{V}:\frac{1}{\sqrt{2}}\Big{|}\sqrt{\frac{12}{K_{t}}}\langle\beta^{(t)}_{k},\mathsf{W}^{(t)}_{v}\rangle+\langle\sigma^{(t)}_{k},\mathsf{D}^{(t)}_{\mathsf{v}}\rangle\Big{|}\geq 10\Big{\}}\,.$

In the above, we introduced a Gaussian smoothing $\{W^{(t)}_{v}(i),\mathsf{W}^{(t)}_{\mathsf{v}}(i):1\leq i\leq K_{t}\}$. We believe this is not essential but provides technical convenience: on the one hand it probably simply reduces the efficiency of the algorithm since it weakens the signal, but on the other hand it facilitates the analysis since it brings the distribution closer to Gaussian. In addition, we have used the absolute value of a random variable instead of a random variable itself, with the purpose of introducing more symmetry as in [17] (e.g., to bound (3.95) below). Recall that we have explained that $\mathrm{M}_{\Gamma}^{(t,t)}$ records the covariance matrix of $D^{(t)}_{v}$ for all $v\in V$. Thus, we expect that the correlation between $\langle\sigma^{(t)}_{k},D^{(t)}_{v}\rangle$ and $\langle\sigma^{(t)}_{l},D^{(t)}_{v}\rangle$ is approximately

$$\frac{12}{K_{t}}\sum_{i,j=1}^{\frac{1}{12}K_{t}}\beta^{(t)}_{k}(i)\beta^{(t)}_{l}(j)\eta^{(t)}_{i}\mathrm{M}_{\Gamma}^{(t,t)}\big{(}\eta^{(t)}_{j}\big{)}^{*}\overset{\eqref{equ-vector-orthogonal},\eqref{equ-vector-unit}}{=}\frac{12}{K_{t}}\sum_{i=1}^{\frac{1}{12}K_{t}}\beta^{(t)}_{k}(i)\beta^{(t)}_{l}(i)=\frac{12}{K_{t}}\big{\langle}\beta^{(t)}_{k},\beta^{(t)}_{l}\big{\rangle}\,.$$

In particular, the variance of each $\langle\sigma^{(t)}_{i},D^{(t)}_{v}\rangle$ is approximately $1$. Similarly, we can show the correlation between $\langle\sigma^{(t)}_{i},\mathsf{D}^{(t)}_{\mathsf{v}}\rangle$ and $\langle\sigma^{(t)}_{j},\mathsf{D}^{(t)}_{\mathsf{v}}\rangle$ is approximately $\frac{12}{K_{t}}\big{\langle}\beta^{(t)}_{k},\beta^{(t)}_{l}\big{\rangle}$, and the correlation between $\langle\sigma^{(t)}_{i},D^{(t)}_{v}\rangle$ and $\langle\sigma^{(t)}_{j},\mathsf{D}^{(t)}_{\pi(v)}\rangle$ is approximately $\hat{\rho}\cdot\frac{12}{K_{t}}\big{\langle}\hat{\beta}^{(t)}_{i},\hat{\beta}^{(t)}_{j}\big{\rangle}$, where

$$\hat{\beta}^{(t)}_{k}(j)=\Big{(}\eta^{(t)}_{j}\Psi^{(t)}\big{(}\eta^{(t)}_{j}\big{)}^{*}\Big{)}^{1/2}\cdot\beta^{(t)}_{k}(j)\,.$$

(2.30)

(Here we also used that (2.16) implies that $\mathrm{P}_{\Gamma,\Pi}^{(t,t)}\approx\Psi^{(t)}$). Recall our desire for (2.13)–(2.16) to hold for $t+1$. Thus the signal contained in each pair at time $t+1$ is approximately

$\displaystyle\varepsilon_{t+1}=\frac{1}{(\mathfrak{a}-\mathfrak{a}^{2})}\Big{(}\phi\Big{(}\frac{\hat{\rho}}{2}\frac{12}{K_{t}}\sum_{j=1}^{\frac{K_{t}}{12}}\eta^{(t)}_{j}\Psi^{(t)}\big{(}\eta^{(t)}_{j}\big{)}^{*}\Big{)}-\phi(0)\Big{)}\,.$

(2.31)

By (2.27), we have that

$$\varepsilon_{t+1}\in\big{[}\frac{\iota_{\mathrm{lb}}\hat{\rho}^{2}}{4(\mathfrak{a}-\mathfrak{a}^{2})}(0.5\varepsilon_{t})^{2},\frac{\iota_{\mathrm{ub}}\hat{\rho}^{2}}{4(\mathfrak{a}-\mathfrak{a}^{2})}(2\varepsilon_{t})^{2}\big{]}\,.$$

(2.32)

Recalling (2.3), we have $\varepsilon_{t+1}\leq\varepsilon_{t}^{2}$, and thus (recall from (2.10) that $\varepsilon_{0}<\tfrac{1}{2}$)

$${}\varepsilon_{t+1}\leq\varepsilon_{t}\leq\ldots\leq\varepsilon_{0}\leq\tfrac{1}{2}\,,$$

(2.33)

which verifies our statement that the signal $\varepsilon_{t}$ in each pair is decreasing. We can then finish the iteration by defining $\Phi^{(t+1)},\Psi^{(t+1)}$ to be $K_{t+1}\!*\!K_{t+1}$ matrices such that

		$\displaystyle\Phi^{(t+1)}(i,j)=(\mathfrak{a}-\mathfrak{a}^{2})^{-1}\Big{\{}\phi\Big{(}\frac{12}{K_{t}}\langle{\beta}^{(t)}_{i},{\beta}^{(t)}_{j}\rangle\Big{)}-\mathfrak{a}^{2}\Big{\}}\,,$		(2.34)
		$\displaystyle\Psi^{(t+1)}(i,j)=(\mathfrak{a}-\mathfrak{a}^{2})^{-1}\Big{\{}\phi\Big{(}\frac{\hat{\rho}}{2}\frac{12}{K_{t}}\langle\hat{\beta}^{(t)}_{i},\hat{\beta}^{(t)}_{j}\rangle\Big{)}-\mathfrak{a}^{2}\Big{\}}\,.$

Next, we state a lemma which then inductively justifies (2.17) and (2.18).