Random Graph Matching in Geometric Models: the Case of Complete Graphs

Given two point clouds $\{X_{1},\ldots,X_{n}\}$ and $\{Y_{1},\ldots,Y_{n}\}$ in ${\mathbb{R}}^{d}$, we construct two weighted graphs on the vertex set $[n]$ with weighted adjacency matrices $A$ and $B$ as follows. For each $i,j$, let $A_{ij}\overset{\rm ind}{\sim}W(\cdot|X_{i},X_{j})$ and $B_{ij}\overset{\rm ind}{\sim}W(\cdot|Y_{i},Y_{j})$, for some probability transition kernel $W$. The coordinates are correlated through a latent matching as follows: Consider a Gaussian model

where $X_{i},Z_{i}$’s are iid ${\mathcal{N}}(0,I_{d})$ vectors and $\pi^{*}$ is uniform on $S_{n}$, the set of all permutations on $[n]$. In matrix form, we have

where $X,Y,Z\in\mathbb{R}^{n\times d}$ are matrices whose rows are $X_{i}$’s, $Y_{i}$’s and $Z_{i}$’s respectively, $\Pi^{*}\in\mathfrak{S}_{n}$ denotes the permutation matrix corresponding to $\pi^{*}$, and $\mathfrak{S}_{n}$ is the collection of all permutation matrices. Given the observation $A$ and $B$, the goal is to recover the latent correspondence $\pi^{*}$.

Of particular interest are the following special cases:

• •

  Dot-product model: The observations are complete graphs with pairwise inner products as edge weights, namely, $A_{ij}=\langle X_{i},X_{j}\rangle$ and $B_{ij}=\langle Y_{i},Y_{j}\rangle$. As such, the weighted adjacency matrices are $A=XX^{\top}$ and $B=YY^{\top}$, both Wishart matrices. It is clear that from $A$ and $B$ one can reconstruct $X$ and $Y$ respectively, each up to a global orthogonal transformation on the rows. In this light, the model is also equivalent to the so-called Procrustes Matching problem [MDK⁺16, DL17, GJB19], where $Y$ in (1) undergoes a further random orthogonal transformation – see Appendix A for a detailed discussion.

• •

  Distance model: The edge weights are pairwise squared distances $A_{ij}=\|X_{i}-X_{j}\|^{2}$ and $B_{ij}=\|X_{i}-X_{j}\|^{2}$. This setting corresponds to the classical problem of multi-dimensional scaling (MDS), where the goal is to reconstruct the coordinates (up to global shift and orthogonal transformation) from the distance data (cf. [BG05]).

• •

  Random Dot Product Graph (RDPG): In this model, the observed data are two graphs with adjacency matrices $A$ and $B$, where $A_{ij}\overset{\rm ind}{\sim}{\rm Bern}\left(\kappa(\left\langle X_{i},X_{j}\right\rangle)\right)$ and $B_{ij}\overset{\rm ind}{\sim}{\rm Bern}\left(\kappa(\left\langle X_{i},X_{j}\right\rangle)\right)$ conditioned on $X$ and $Y$, and $\kappa:{\mathbb{R}}\to[0,1]$ is some link function, e.g. $\kappa(t)=e^{-t^{2}/2}$. In this way, we observe two instances of RDPG that are correlated through the underlying points and the latent matching. See [AFT⁺17] for a recent survey on RDPG.

• •

  Random Geometric Graph (RGG): Similar to RDPG, $A_{ij}\overset{\rm ind}{\sim}{\rm Bern}(\kappa(\|X_{i}-X_{j}\|))$ conditioned on $X_{1},\ldots,X_{n}$ for some link function $\kappa:{\mathbb{R}}_{+}\to[0,1]$ applied to the pairwise distances. The second RGG instance $B$ is constructed in the same way using $Y_{1},\ldots,Y_{n}$. A simple example is $\kappa(t)={\mathbf{1}_{\left\{{t\leq r}\right\}}}$ for some threshold $r>0$, where each pair of points within distance $r$ is connected [Gil61]; see the monograph [Pen03] for a comprehensive discussion on RGG.

Let us mention that the model where the two point clouds are directly observed has been recently studied by [DCK19, DCK20] in the context of feature matching and independently by [KNW22] as a geometric model for the planted matching problem, extending the previous work in [CKK⁺10, MMX21, DWXY21] with iid weights to a geometric (low-rank) setting. In this model, $X$ and $Y$ in (2) are observed and the maximum likelihood estimator (MLE) of $\pi^{*}$ amounts to solving

which is a linear assignment (max-weight matching) problem on the weighted complete bipartite graph with weight matrix $YX^{\top}$. In the sequel we shall refer to this setting as the linear assignment model, which we also study in this paper for the sake of proving impossibility results for the more difficult graph matching problem, as the coordinates are latent and only pairwise information are available.

Fig. 1 elucidates the logical connections between the aforementioned models. Among these, linear assignment model is the most informative, followed by the dot product model and the distance model, whose further stochastically degraded versions are RDPG and RGG, respectively. As a first step towards understanding graph matching on geometric models, in this paper we study the case of weighted complete graphs in the dot product and distance models.

By analyzing the MLE (2) in the stronger linear assignment model (1), [KNW22] identified a critical scaling of dimension $d$ at $\log n$:

• •

  In the low-dimensional regime of $d\ll\log n$, accurate reconstruction requires the noise level $\sigma$ to be vanishingly small. More precisely, with high probability, the MLE (2) recovers the latent $\pi^{*}$ perfectly (resp. with a vanishing fraction of errors) provided that $\sigma=o(n^{-2/d})$ (resp. $\sigma=o(n^{-1/d})$).

• •

  In the high-dimensional regime of $d\gg\log n$, it is possible for $\sigma^{2}$ to be as large as $\frac{d}{(4+o(1))\log n}$. Since the dependency between the edges weakens as the latent dimension increases,¹¹1For the Wishart matrix, it is known [JL15, BG18] that the total variation between the joint law of the off-diagonals and their iid Gaussian counterpart converges to zero provided that $d=\omega(n^{3})$. Analogous results have also been obtained in [BDER16] showing that high-dimensional RGG is approximately Erdős-Rényi. this is consistent with the known results in the correlated Erdős-Rényi and Wigner model. For example, to match two GOE matrices with correlation coefficient $\rho$, the sharp reconstruction threshold is at $\rho^{2}=\frac{(4+o(1))\log n}{n}$ [Gan21b, WXY21].

In this paper we mostly focus on the low-dimensional setting as this is the regime where geometric graph ensembles are structurally distinct from Erdős-Rényi graphs. Our main findings are two-fold:

1. 1.

   The same reconstruction thresholds remain achievable even when the coordinates are latent and only inner-product or distance data are accessible.

2. 2.

   Furthermore, these thresholds cannot be improved even when the coordinates are observed.

To make these results precise, we start with the dot-product model with $A=XX^{\top}$ and $B=YY^{\top}$, and $Y=\Pi^{*}X+\sigma Z$ according to (1). In this case the MLE turns out to be much more complicated than (2) for the linear assignment model. As shown in Appendix B, the MLE takes the form

where the integral is with respect to the Haar measure on the orthogonal group $O(d)$, $A^{1/2}\triangleq U\Lambda^{1/2}\in{\mathbb{R}}^{n\times d}$ based on the SVD $A=U\Lambda U^{\top}$, and similarly for $B^{1/2}$. It is unclear whether the above Haar integral has a closed-form solution,²²2The integral in (3) can be reduced to computing $\int{\rm d}Q\exp(\langle\Lambda,Q\rangle)$ for a diagonal $\Lambda$, which, in principle, can be evaluated by Taylor expansion and applying formulas for the joint moments of $Q$ in [Mat13, Theorem 2.2]. let alone how to optimize it over all permutations. Next, we turn to its approximation.

As we will show later, in the low-dimensional case of $d=o(\log n)$, meaningful reconstruction of the latent matching is information-theoretically impossible unless $\sigma$ vanishes with $n$ at a certain speed. In the regime of small $\sigma$, Laplace’s method suggests that the predominant contribution to the integral in (3) comes from the maximum $\langle B^{1/2},\Pi A^{1/2}Q\rangle$ over $Q\in O(d)$. Using the dual form of the nuclear norm $\|X\|_{*}=\max_{Q\in O(d)}\langle X,Q\rangle$, where $\|X\|_{*}$ denotes the sum of all singular values of $X$, we arrive at the following approximate MLE:

We stress that the above approximation to the MLE (3) is justified for the low-dimensional regime where $\sigma$ is small. In the high-dimensional (high-noise) case, the approximate MLE actually takes on the form of a quadratic assignment problem (QAP), which is the MLE for the well-studied iid model [CK16]; in the special case of the dot-product model, it amounts to replacing the nuclear norm in (4) by the Frobenius norm. We postpone this discussion to Section 4.

To measure the accuracy of a given estimator $\widehat{\pi}$, we define

as the fraction of nodes whose matching is correctly recovered. The following result identifies the threshold at which the approximate MLE achieves perfect or almost perfect recovery.

A few remarks are in order:

• •

  In fact we will show the following nonasymptotic estimate that implies (6): For all sufficiently small $\varepsilon$, if $\sigma^{-d}>16n2^{2/\varepsilon}$, then $\mathsf{overlap}(\widehat{\pi}_{\mathrm{AML}},\pi^{*})\geq 1-\varepsilon$ with probability tending to one.

• •

  The estimator (4) has previously appeared in the literature of Procrustes matching [GJB19], albeit not as an approximation to the MLE in a generative model. See Appendix A for a detailed discussion.

• •

  Unlike linear assignment, it is unclear how to solve the optimization in (4) over permutations efficiently. Nevertheless, for constant $d$ we show that it is possible to find an approximate solution in time that is polynomial in $n$ that achieves the same statistical guarantee as in Theorem 1. Indeed, note that (3) is equivalent to the double maximization

  $$\widehat{\Pi}_{\mathrm{AML}}=\arg\max_{\Pi\in\mathfrak{S}_{n}}\max_{Q\in O(d)}\langle B^{1/2},\Pi A^{1/2}Q\rangle.$$

  (7)

  Approximating the inner maximum over a suitable discretization of $O(d)$, each maximization over $\Pi$ for fixed $Q$ is a linear assignment problem, which can be solved in $O(n^{3})$ time. In Section 3, we provide a heuristic that shows (7) can be further approximated by the classical spectral algorithm of Umeyama [Ume88] which is much faster in practice and achieves good empirical performance. For $d$ that grows with $n$, it is an open question to find a polynomial-time algorithm that attains the (optimal, as we show next) threshold in Theorem 1.

Next, we proceed to the more difficult distance model, where $A_{ij}=\|X_{i}-X_{j}\|^{2}$ and $B_{ij}=\|Y_{i}-Y_{j}\|^{2}$. Deriving the exact MLE in this model appears to be challenging; instead, we apply the estimator (4) to an appropriately centered version of the data matrices. Let $\mathbf{1}\in{\mathbb{R}}^{n}$ denotes the all-one vector and define $\mathbf{F}=\frac{1}{n}\mathbf{1}\mathbf{1}^{\top}$. Then $A=-2XX^{\top}+a\mathbf{1}^{\top}+\mathbf{1}a^{\top}$ and $B=-2YY^{\top}+b\mathbf{1}^{\top}+\mathbf{1}b^{\top}$, where $a=(\|X_{i}\|^{2})$ and $b=(\|Y_{i}\|^{2})$. Strictly speaking, the vectors $a$ and $b$ are correlated with the ground truth $\pi^{*}$, since $b$ can be viewed as a noisy version of $\Pi^{*}a$; however, we expect them to inform very little about $\pi^{*}$ because such scalar-valued observations are highly sensitive to noise (analogous to degree matching in correlated Erdős-Rényi graphs [DMWX21, Section 1.3]). As such, we ignore $a$ and $b$ by projecting $A$ and $B$ to the orthogonal complement of the vector $\mathbf{1}$. Specifically, we compute, as commonly done in the MDS literature (see e.g. [SRZF03, OMK10]),

It is easy to verify that $\widetilde{A}=\widetilde{X}\widetilde{X}^{\top}$ and $\widetilde{B}=\widetilde{Y}\widetilde{Y}^{\top}$, where $\widetilde{X}=(I-\mathbf{F})X$ and $\widetilde{Y}=(I-\mathbf{F})Y$ consist of centered coordinates $\widetilde{X}_{i}=X_{i}-\bar{X}$ and $\widetilde{Y}_{i}=Y_{i}-\bar{Y}$ respectively, with $\bar{X}=\frac{1}{n}\sum_{i=1}^{n}X_{i}$ and $\bar{Y}=\frac{1}{n}\sum_{i=1}^{n}Y_{i}$. Overall, we have reduced the distance model to a dot product model where the latent coordinates are now centered.

One can show that the MLE of $\Pi^{*}$ given the reduced data $(\widetilde{A},\widetilde{B})$ is of the same Haar-integral form (3). Using again the small-$\sigma$ approximation, we arrive at the following estimator by applying (4) to the centered data $\widetilde{A}$ and $\widetilde{B}$:

Finally, we state an impossibility result for the linear assignment model, proving that the perfect and almost perfect recovery threshold of $\sigma=o(n^{-2/d})$ and $\sigma=o(n^{-1/d})$ obtained by analyzing the MLE in [KNW22] are in fact information-theoretically necessary. Complementing Theorem 1 and Theorem 2, this result also establishes the optimality of the estimator (4) and (9) for their respective model.

Theorem 3(i) slightly improves the necessary condition for perfect recovery in [KNW22] from $\sigma=O(n^{-2/d})$ to $\sigma=o(n^{-2/d})$. For almost perfect recovery, the negative result in [KNW22] is limited to MLE, while Theorem 3 holds for all algorithms. Moreover, the necessary condition in Theorem 3(ii) was conjectured in [KNW22, Conjecture 1.4, item 1], which we now resolve in the positive. Finally, while our focus is in the low-dimensional case of $d=o(\log n)$, we also provide necessary conditions that hold for general $d$. (See Appendix E for details).

In view of Fig. 1, since the negative results in Theorem 3 are proved for the strongest model and the positive results in Theorem 2 are for the weakest model, we conclude that for all three models, namely, linear assignment, dot-product, and distance model, the thresholds for exact and almost perfect reconstruction is given by $n^{-2/d}$ and $n^{-1/d}$, respectively.

The positive results of Theorem 1 and Theorem 2 are proved in Appendix C and Appendix D. Here we briefly describe the proof strategy in the dot product model. Suppose we want to bound the probability that the approximate MLE $\widehat{\Pi}_{\mathrm{AML}}$ in (4) makes more than $t$ number of errors. Denote by ${\rm d}(\pi_{1},\pi_{2})\triangleq\sum_{i=1}^{n}{\mathbf{1}_{\left\{{\pi_{1}(i)\neq\pi_{2}(i)}\right\}}}$ the Hamming distance between two permutations $\pi_{1},\pi_{2}\in S_{n}$. Without loss of generality, we will assume that $\pi^{*}=\mathrm{Id}$. By the orthogonal invariance of $\|\cdot\|_{*}$, we can assume, for the sake of analysis, that $A^{1/2}=X$ and $B^{1/2}=Y$. Applying (7),

For each fixed $\Pi$ and $Q$, averaging over the noise yields, for some absolute constant $c_{0}$,

In the remaining argument, there are three places where the structure of the orthogonal group $O(d)$ plays a crucial role:

1. 1.

   The quantity in (11) turns out to depend on $\Pi$ through its cycle type and on $Q$ through its eigenvalues. Crucially, the eigenvalues of an orthogonal matrix $Q$ lie on the unit circle, denoted by $(e^{\mathrm{i}\theta_{1}},\ldots,e^{\mathrm{i}\theta_{d}})$, with $|\theta_{\ell}|\leq\pi$. We then show that the error probability in (11) can be further bounded by, for some absolute constant $C_{0}$,

   $$(C_{0}\sigma)^{d(n-{\mathfrak{c}})}\left(\prod_{\ell=1}^{d}\frac{C_{0}\sigma}{\sigma+|\theta_{\ell}|}\right)^{n_{1}},$$

   (12)

   where $n_{1}$ is the number of fixed points in $\pi$ and ${\mathfrak{c}}$ is the total number of cycles.

2. 2.

   In order to bound (10), we take a union bound over $\pi$ and another union bound over an appropriate discretization of $O(d)$. This turns out to be much subtler than the usual $\delta$-net-based argument, as one needs to implement a localized covering and take into account the local geometry of the orthogonal group. Specifically, note that the error probability in (11) becomes larger when $\pi$ is near $\mathrm{Id}$ and when $Q$ is near $I_{d}$ (i.e. the phases $|\theta_{\ell}|$’s are small); fortunately, the entropy (namely, the number of such $\pi$ and such $Q$ within a certain resolution) also becomes smaller, balancing out the deterioration in the probability bound. This is the second place where the structure of $O(d)$ is used crucially, as the local metric entropy of $O(d)$ in the vicinity of $I_{d}$ is much lower than that elsewhere.

3. 3.

   Controlling the approximation error of the nuclear norm is another key step. Note that for any matrix norm of the dual form $\|A\|=\sup_{\|Q\|^{\prime}\leq 1}\langle A,Q\rangle$, where $\|\cdot\|^{\prime}$ is the dual norm of $\|\cdot\|$, the standard $\delta$-net argument (cf. [Ver18, Lemma 4.4.1]) yields a multiplicative approximation $\max_{Q\in N}\left\langle A,Q\right\rangle\geq(1-\delta)\|A\|$, where $N$ is any $\delta$-net of the dual norm ball. In general, this result cannot be improved (e.g. for Frobenius norm); nevertheless, for the special case of nuclear norm, this approximation ratio can be improved from $1-\delta$ to $1-\delta^{2}$, as the following result of independent interest shows. This improvement turns out to be crucial for obtaining the sharp threshold.

   Lemma 1.

   Let $N\subset O(d)$ be a $\delta$-net in operator norm of the orthogonal group $O(d)$. For any $A\in\mathbb{R}^{d\times d}$,

   $$\max_{Q\in N}\left\langle A,Q\right\rangle\geq\left(1-\frac{\delta^{2}}{2}\right)\|A\|_{*}.$$

   (13)

The proof of Theorem 1 is completed by combining (12) with a union bound over a specific discretization of $O(d)$, whose cardinality satisfies the desired eigenvalue-based local entropy estimate, followed by a union bound over $\pi$ which can be controlled using moment generating function of the number of cycles in a random derangement.

The information-theoretic lower bounds in Theorem 3 for the linear assignment model are proved in Appendix E. Here we sketch the main ideas. We first derive a necessary condition for almost perfect recovery that holds for any $d$ via a simple mutual information argument [HWX17]: On one hand, the mutual information $I(\pi^{*};X,Y)$ can be upper bounded by the Gaussian channel capacity as $\frac{nd}{2}\log(1+\sigma^{-2})$. On the other hand, to achieve almost perfect recovery, $I(\pi^{*};X,Y)$ needs be asymptotically equal to the full entropy $H(\pi^{*})$ which is $(1-o(1))\log n$. These two assertions together immediately imply that $\frac{nd}{2}\log(1+\sigma^{-2})\geq((1-o(1))\log n$, which further simplifies to $\sigma=n^{-(1-o(1))/d}$ when $d=o(\log n)$. However, for constant $d$, this necessary condition turns out to be loose and the main bulk of our proof is to improve it to the optimal condition $\sigma=o(n^{-1/d})$. To this end, we follow the program recently developed in [DWXY21] in the context of the planted matching model by analyzing the posterior measure of the latent $\pi^{*}$ given the data $(X,Y)$.

To start, a simple yet crucial observation in [DWXY21] is that to prove the impossibility of almost perfect recovery, it suffices to show a random permutation sampled from the posterior distribution is at Hamming distance $\Omega(n)$ away from the ground truth with constant probability. As such, it suffices to show there is more posterior mass over the bad permutations (those far away from the ground truth) than that over the good permutations (those near the ground truth) in the posterior distribution. To proceed, we first bound from above the total posterior mass of good permutations by a truncated first moment calculation applying the large deviation analysis developed in the proof of the positive results. To bound from below the posterior mass of bad permutations, we aim to construct exponentially many bad permutations $\pi$ whose log likelihood $L(\pi)$ is no smaller than $L(\pi^{*})$. A key observation is that $L(\pi)-L(\pi^{*})$ can be decomposed according to the orbit decomposition of $(\pi^{*})^{-1}\circ\pi$:

where ${\mathcal{O}}$ denotes the set of orbits in $(\pi^{*})^{-1}\circ\pi$ and for any orbit $O=(i_{1},i_{2},\ldots,i_{t})$,

Thus, the goal is to find a collection of vertex-disjoint orbits $O$ whose total lengths add up to $\Omega(n)$ and each of which is augmenting in the sense that $\Delta(O)\geq 0$. Here, a key difference to [DWXY21] is that in the planted matching model with independent edge weights studied there, short augmenting orbits are insufficient to meet the $\Omega(n)$ total length requirement; instead, [DWXY21] resorts to a sophisticated two-stage process that first finds many augmenting paths then connects then into long cycles. Fortunately, for the linear assignment model in low dimensions of $d=\Theta(1)$, as also observed in [KNW22] in their analysis of the MLE, it suffices to look for augmenting $2$-orbits and take their disjoint unions. More precisely, we show that there are $\Omega(n)$ many vertex-disjoint augmenting $2$-orbits. This has already been done in [KNW22] using a second-moment method enhanced by an additional concentration inequality. It turns out that the correlation among the augmenting $2$-orbits is mild enough so that a much simpler argument via a basic second-moment calculation followed by an application of Turán’s theorem suffices to extract a large vertex-disjoint subcollection. Finally, these vertex-disjoint augmenting $2$-orbits give rise to exponentially many permutations that differ from the ground truth by $\Omega(n)$.

Finally, we briefly remark on perfect recovery, for which it suffices to focus on the MLE (2) which minimizes the error probability for uniform $\pi^{*}$. In view of the likelihood decomposition given in (14), it further suffices to prove the existence of an augmenting $2$-orbit. This can be easily done using the second-moment method. A similar strategy was adopted in [DCK19], but our first-moment and second-moment estimates are tighter and hence yield nearly optimal conditions.