SIGMA: A Structural Inconsistency Reducing Graph Matching Algorithm

Graph matching finds the correspondence of nodes across two correlated graphs and lies at the core of many applications. When graph side information is not available, the node correspondence is estimated on the sole basis of network topologies. In this paper, we propose a novel criterion to measure the graph matching accuracy, structural inconsistency (SI), which is defined based on the network topological structure. Specifically, SI incorporates the heat diffusion wavelet to accommodate the multi-hop structure of the graphs. Based on SI, we propose a Structural Inconsistency reducing Graph Matching Algorithm (SIGMA), which improves the alignment scores of node pairs that have low SI values in each iteration. Under suitable assumptions, SIGMA can reduce SI values of true counterparts. Furthermore, we demonstrate that SIGMA can be derived by using a mirror descent method to solve the Gromov-Wasserstein distance with a novel K-hop-structure-based matching costs. Extensive experiments show that our method outperforms state-of-the-art methods.


Introduction
Graph matching (network alignment) aims to identify the correspondence of nodes across two related graphs, by minimizing the sum of pair-wise matching costs of corresponding nodes/edges. It has many applications, such as linking user accounts in social platforms [Shu et al., 2017, Wang et al., 2019b, aligning entities in knowledge graphs [Sun et al., 2020, Pei et al., 2020, and matching keypoints across two images [Sarlin et al., 2020, Wang et al., 2020a, to name a few.
The choice of matching costs is critical to the performance of graph matching. When node/edge side information, e.g. attributes, is available, one can rely on expert knowledge to design handcrafted node embeddings for the construction of matching costs Tong, 2016, Meng et al., 2016]. Recently, to alleviate the requirement of expert knowledge, end-to-end deep learning methods are proposed to learn the node embedding automatically based on the anchor links (labelled node correspondences) [Zanfir and Sminchisescu, 2018, Fey et al., 2019, Wang et al., 2019a, Sarlin et al., 2020, Wang et al., 2020b. However, in many tasks, side information may not be available Krevl, 2014, Heimann et al., 2021]. Under such circumstance, the matching costs can merely be estimated on the basis of the network topologies of the two graphs [Patro and Kingsford, 2012, Saraph and Milenković, 2014, Sun et al., 2015, Mémoli, 2011, Maretic et al., 2019, Xu et al., 2019b, Heimann et al., 2021.
Among these network-topology based methods, one line of works align graphs by borrowing tools from optimal transport (OT) and achieve the state-of-the-art performance [Mémoli, 2011, Maretic et al., 2019, Xu et al., 2019b, Titouan et al., 2019, Chowdhury and Mémoli, 2019, Barbe et al., 2020. Basically, these methods utilize OT to exploit the geometric properties of the metric spaces that underlie the graphs, and estimate the node correspondence by calculating the transport plan of the Gromov-Wasserstein distance. Typically, the matching costs in these methods are calculated depending on the one-hop neighborhood information of each node.
Recently, Heimann et al. [2021] reveals the connection between the graph matching accuracy and the matched neighborhood consistency (MNC), i.e., the Jaccard simliarity between the one-hop neighborhoods of matched node pairs. In particular, they show that accurate matching essentially entails high MNC. However, it remains unknown how the OT-based methods are related to MNC theoretically. Instead, they propose a method, RefiNA, by maximizing the numerators of MNC, which can be regarded as an opposite measure of the matching cost. Though the authors intend to find matching with high MNC, RefiNA actually diverges from their purpose, as the denominator of MNC is ignored and a high numerator does not ensure a high MNC. Moreover, both OT-based methods and RefiNA may suffer from one-hop structural indistinguishability and tend to misalign nodes that have low degrees, as they discard the abundant network topological information and estimate the graph matching merely depending on the one-hop neighborhoods structure.
In this paper, we propose a novel criterion to measure the graph matching accuracy, structural inconsistency (SI), which utilizes the abundant network topological information. Specifically, SI incorporates the heat diffusion wavelet [Donnat et al., 2018] to accommodate the multi-hop structure of the graphs. If a pair of nodes incurs zero SI, they are K-hop indistinguishable. On the other hand, SI values of true counterparts should be zero for two isomorphic graphs. In practice, the SI between a true corresponding pairs may only close to 0 due to the structural noise.
Based on SI, we propose a Structural Inconsistency reducing Graph Matching Algorithm (SIGMA), which improves the alignment scores of node pairs that have low SI values in each iteration. Under suitable assumptions, SIGMA can reduce SI values of true counterparts. As SI considers the multi-hop topological information, SIGMA can avoid the misguidance of one-hop structural indistinguishability and find graph matching of higher accuracy. Besides, we also show that our method increases a soft version of MNC values, which is positively related to MNC and should also be high/low for accurate/inaccurate matching. Furthermore, we demonstrate that SIGMA can be derived by using a mirror descent method to solve the Gromov-Wasserstein distance with a novel K-hop-structure-based matching costs. As a result, SIGMA can be regarded as the first OT-based method that are guaranteed to increase MNC in each iteration in the graph matching literature.
Extensive experimental results demonstrate that the proposed method outperforms state-of-the-art methods. The rest of the paper is organized as follows. In Sec. 2, we review the background and discuss the relationship between SI and matching accuracy. Based on SI, a new graph matching method SIGMA is proposed in Sec. 4. We give theoretical analysis in Sec. 5. Due to the limited space, all proof are given in the appendix. Experimental results are demonstrated in Sec. 6 and the appendix.

Preliminary
In this section, we give basic preliminary to facilitate the later discussion and review the main background.
Here, we first give the notations used in the main text. We use bold lowercase symbols, bold uppercase letters, and uppercase calligraphic fonts to denote vectors (e.g. x), matrices (e.g. A), and sets (e.g. A), respectively. (·) is the transpose of a vector or a matrix. A[i, :] and A [:, j] are the i-th row and the j-th column of matrix A respectively. For two matrices A and B that are of the same size, A, B = trace(A B) is the Frobenius dot-product. We denote the cardinality of set A by |A|.

Graph Matching
We consider undirected graphs without self-loops G = (V, E, W), where V is the set of vertices, E is the set of edges, and W = [W ij ] ∈ R |V|×|V| is a (binary or weighted) adjacency matrix. If there is an edge (i, j) connecting vertices i and j, i.e., (i, j) ∈ E, the entry W ij represents the weight of the edge; Otherwise, W ij = 0. The degree of each vertex i, written as d i , is defined as the sum of the weights of all the edges incident to it, i.e., d i = |V| j=1 W ij . The degree matrix D = [D ij ] has diagonal elements equal to the degrees D ii = d i , and zeros elsewhere. The Laplacian matrix of G is defined as L = D − W. The k-hop neighborhood of node i, denoted by N k (i), is the subgraph of G induced by all nodes that are k or fewer hops away from node i.
Given the source graph G s = (V s , E s , W s ) and the target graph G t = (V t , E t , W t ), assuming |V s |≤ |V t | without loss of generality graph matching aims to find an injective mapping π : V s → V t that matches nodes in G s to nodes in G t . In this paper, we focus on equal-size graph matching, i.e., V = |V s |= |V t |. In such case, π is restricted to be a bijective. For graphs with different number of nodes, one can add dummy isolated nodes into the smaller graph and make them equal-sized [Gold and Rangarajan, 1996, Zaslavskiy et al., 2008, Jiang et al., 2017. π is also commonly represented as a matching matrix [Fey et al., 2019, Heimann et al., 2021), where

Optimal Transport
Given two discrete measures α = m i=1 p i δ xi (x) and β = n j=1 q j δ yj (y), Optimal transport (OT) addresses the problem of optimally transporting Villani, 2008]. The p-Wasserstein distance between p and q is defined as where K ij is the distance between x i and y j , K p = [K p ij ] ∈ R m×n + and the feasible domain of transport plan T is given by the set of coupling measures Π(p, q) = {T ∈ R m×n + |T1 n = p, T 1 m = q}. Gromov-Wasserstein (GW) distance is a generalization of Wasserstein distance [Mémoli, 2011]. Let X and Y be two sample spaces. Endowing the spaces X and Y with metrics (distances) d X and d Y , the GW distance is defined as . , x m ∈ X and y 1 , . . . , y n ∈ Y.
Optimal transport can be applied to graph matching. The source graph G s and the target graph G t are modelled as two probability distributions µ s = [µ s i ] and µ t = [µ t i ] respectively, where A dissimilarity matrix B s can be assigned to G s to encode the structural information. Each entry B s ij models the distance between nodes i and j in G s , such as the edge weight-based distance [Xu et al., 2019b], the shortest path distance, or the harmonic distance [Verma and Zhang, 2017]. B t is defined similarly. Basically, the node correspondence is estimated by calculating the transport plan of the Gromov-Wasserstein distance that is defined as Node i ∈ V s is matched toî ∈ V t , whereî = argmax i T * ii and T * attains the minimum of Eq. (2).
Note that to better approximate matching matrices (1), one can relax the The matching matrix can then be interpreted as a transport plan from 1 to 1.

Heat Diffusion Wavelet
For each graph z = s, t, the heat diffusion on a graph is the solution to the discrete heat equation known as the Laplacian exponential diffusion kernel [Hammond et al., 2013]. As t approaches zero, high-order terms vanish, meaning that the kernel depicts local connectivity [Tsitsulin et al., 2018, Donnat et al., 2018. The K-th order approximation of exp(−tL z ) is The i-th row, Ψ[i, :], is the K-th order HDW for node i [Donnat et al., 2018] and measures the network connectivity between node i and each node j. B z = [B z ij ] that characterizes node dissimilarities within each graph is defined as where Ψ z ij is defined as Eq. (3), z = s, t and constantΨ > max i,j,z Ψ z ij .

Structural Inconsistency
Based on the local neighborhoods of node i and node i , they can be embedded into a V -dimensional space, which yields embedding vectors z i and z i respectively. The j-th dimension of z i (resp. j -th dimension of z i ) models the distance between i and j (resp. i and j ). In this paper, we instantiate the embeddings based on HDW, i.e., Given transport plan T ∈ Π(1, 1), we model the structural inconsistency (SI) of a pair of nodes (i, i ) ∈ V s × V t as a transport distance   Figure 1a and thus structural inconsistency occurs.
We illustrate SI(A, a) when K = 1 in Figure 1. Here the transport plans are depicted by the thin dashed lines and we omit them in SI. In Figure 1a, node B (resp. C) is wrongly matched to node c (resp. b). Nodes A and B (a and b) are connected by an edge while no edge connects a and c(A and C), which is depicted by the two red(yellow) circles. Therefore, structural inconsistency incurs. In Figure  1b, node A and its one-hop neighbors are correctly matched, which leads to SI(A, a) = 0.
One can derive a SI matrix. The SI for all node pairs (i, i ) given the matching matrix T, can be written where h denotes the element-wise square and ⊗ is the outer product of two vectors. We further define the structural inconsistency of true counterparts (CSI) as CSI(i; T) = SI(i, π * (i); T), where π * is the true matching. SI is closely related to alignment accuracy. On the one hand, if node pair (i,î) have zero SI values,î either is the true counterpart π * (i) or have the same local topology as π * (i), which is formalized in the following theorem.
Theorem 1 Assume G s and G t are isomorphic and there exists T that yields SI(i,î; T) = 0. Then, ifî = π * (i), N t k (î) and N t k (π * (i)) are isomorphic for all k ≤ K.
An example of Theorem 1 is matching two star graphs, each consisting of one center node connected to n − 1 peripheral nodes [Heimann et al., 2021]. Matching peripheral nodes to each other in any order leads to zero SI's, whatever the true correspondence of the peripheral nodes is. This inspires us to use larger neighborhoods to improve the performance, since N t k+1 (î) and N t k+1 (π * (i)) may not be isomorphic when N t k (î) and N t k (π * (i)) are. On the other hand, if G s is isomorphic to G t up to a few noisy or missing edges (the vetex sets remain the same), CSI (i; T * ) values are low where T * = [T * ii ] is matrix form of true matching. This is stated in the following theorem.
Theorem 2 Let ∆ k ∈ R V ×V denote the perturbation to the edge structure, which is given by ∆ k = (L s ) k − (L t ) k whereL t = T * L t T * is the Laplacian matrix of the registered target graph. Assume that for all i, the perturbation to the local neighborhood is bounded as k = max i j∈N s Then CSI (i; T * ) can be bounded from above as follows: where t is the propagation time.
Note that the propagation time t is often set to a small value in order to capture the local topology [Tsitsulin et al., 2018, Donnat et al., 2018. Thus, for a fixed K, when the perturbation k is small, CSI(i; T * ) has a small upper bound and even is close enough to 0. It is reasonable to conjecture that keeping decreasing CSI will lead to an optimal result, which is embodied in our method design.

Methodology
We update the matching matrix based on the SI matrix. The high level intuition is that each iteration decreases alignment scores for node pairs that have relatively high SI values. The recurrence of SIGMA reads where τ = 0, 1, . . . , is the element-wise multiplication and the exponentiation is also element-wise. The projection can be solved via Sinkhorn-Knopp algorithm with linear convergence rate [Sinkhorn and Knopp, 1967]. Actually, iteration (7) is the mirror descent update for calculating the GW distance. To solve problem the most common algorithm is gradient descent that iteratively updates T as follows where the gradient reads Mirror descent takes into account the geometry of the feasible domain [Bubeck et al., 2015], and is better-suited to Problem (8). Mirror descent is given by the recurrence where D ψ (T, T (τ ) ) stands for the Bregman divergence. When one selects the generalized KL-divergence KL(A|B) = ij A ij log Aij Bij − A ij + B ij as the Bregman divergence, (11) becomes iteration (7).
Algorithm 1 SIGMA 1: Input: Total rounds T , graphs G s and G t , initial transport plan T (0) . 2: Output: Correspondence set C. Update T (τ +1) by iteration 7. 6: end for 7:T = T (T ) 8: Initialize correspondence set C = ∅ 9: for i ∈ V s do 10: SIGMA is summarized in Algorithm 1. The matching cost is based on heat diffusion wavelets. It updates the matching matrix using (7). By choosing the largest index fromT[i, :] for each i, we find the correspondence.

Theoretical Analysis
In this section, we first prove that iteration (7) can reduce CSI(i; T) for all i under suitable assumptions. Then the convergence of SIGMA is proved. We also discuss its relationship to state-of-the-art methods.
With appropriate η (τ ) , iteration 7 can reduce alignment score T ii if SI(i, i ; T (τ ) ) is high, which is formalized in the following lemma.
We now state our main theorem.

Combining Lemma 3 and Theorem 4, if we have
) for all i and i =î, we can reduce CSI(i; T) at iteration τ .
By establishing the relationship between SIGMA and optimal transport based methods, we can borrow the convergence result of mirror descent and prove that SIGMA converges as follows.
Theorem 5 SIGMA converges to an -stationary point with the number of mirror descent iterations O( 1 2 ). The main steps of the proof follow Zhang and He [2018]. We further remark that the true matching matrix T * is a fixed point of update (7).
Computational complexity. The complexity for calculaing the heat diffusion wavlets is O(KV 3 ). Due to the two matrix multiplication steps of the term B s TB t , calculating the gradient involves computational complexity O(V 3 ). Note that the matrix multiplication is highly parallelizable and well-suited to modern computing architectures like GPUs. O(V 2 ) cost is induced by the matrix-vector multiplication in the projection operation. Therefore, if we run projected mirror descent for N iterations, each of which involves T matrix-vector multiplications in the projection, O(N T V 2 + (N + K)V 3 ) computational cost is required in learning the matching matrix. Such complexity is of the same order as state-of-the-art methods like Xu et al. [2019b], Konar

Connection to GWL Xu et al. [2019b]
The gradient that arises from Problem (2) [Xu et al., 2019b,a, Titouan et al., 2019, Barbe et al., 2020 and has the following form: where h denotes the element-wise square and ⊗ is the outer product of two vectors [Peyré et al., 2016]. Note that if we set µ z i = 1 for all i and z = s, t, (12) is identical to (10).
GWL uses proximal point method to solve Problem (2) [Xu et al., 2019b]. If one single mirror descent step is adopted to solve the sub-problem and the marginals are set as µ z i = 1, the proximal point method is identical to our mirror descent method.

Connection to RefiNA Heimann et al. [2021]
When K = 1, the pairwise SI has a more intuitive interpretation.
Proposition 6 When one only uses one-hop neighborhood information (i.e., K = 1), is the set of nodes onto which π maps i's neighborsÑ s k (i; π) = {j ∈ V t |∃j ∈ N s k (i) s.t. π(j) = j }, and the bijective is given by π(i) = argmax i [T ii ].

Chen et al. 2020 and Heimann et al. 2021 define a similar quantity.
Specifically, the matched neighborhood consistency (MNC) of node i in G s and node i in G t is the Jaccard similarity of the setsÑ s 1 (i; π) and N t 1 (i ), i.e., MNC(i, i ; π) = |Ñ s 1 (i;π)∩N t 1 (i )| |Ñ s 1 (i;π)∪N t 1 (i )| . SI(i, i ; T) and MNC(i, i ; π) can be related as follows.
Proposition 7 When propagation time t = 1 andÑ s 1 (i; π)∩N t 1 (i ) is not empty, According to Proposition 7, Lemma 3 and Theorem 4, SIGMA with K = 1 and t = 1 can also be viewed as increasing a soft version of MNC values, which is positively related to MNC and should also be high/low for accurate/inaccurate matching. Combined with the description in Sec. 4, SIGMA can be regarded as the first OT-based method that are guaranteed to increase MNC in each iteration in the graph matching literature. In addition, by characterizing the structure of a larger neighborhood via HDW, we alleviate the problem of near structural indistinguishability and improves the matching performance.

Experiments
In this section, we compare SIGMA and state-of-the-art methods. Due to the limited space, additional experimental results are provided in the appendix.

Experimental Setup
Computing infrastructure. All codes are implemented in Python 3.6 and the package dependencies are listed in the requirements.txt file in the code submission. The experiments are conducted on a Ubuntu server with two CPUs (Intel Xeon (R) CPU E5-2680 v4 @ 2.40GHz), 4 Nvidia 2080 Ti graphics cards, and 378 GB of RAM.
Baselines. Our baselines are unsupervised methods including (1) RefiNA [Heimann et al., 2021] that also explicitly considers the neighborhood consistency; (2) REGAL  that obtains the node embeddings based on matrix factorization; (3) GRAMPA [Fan et al., 2020] that is proposed to match correlated Erdős-Rényi graphs and involves computing all eigenvalues and the associated eigenvectors of the adjacency matrices; (4) MM [Konar and Sidiropoulos, 2020] that formulates graph matching as maximizing a monotone supermodular set function subject to matroid intersection constraints; (5) GWL [Xu et al., 2019b] that matches graphs based on OT; (6) GDD [Scott and Mjolsness, 2021] that involves solving a three-level nested optimization problem. As is shown in their papers, RefiNA and MM are sensitive to the noise and thus take the output of CONE-Align [Chen et al., 2020] as initialization. For fair comparison, GWL and SIGMA also use the output of CONE-Align to initialize the transport plan. We implement REGAL, GWL, and GDD based on the corresponding open-sourced code.
Evaluation Metric. We report node correctness (NC) that is defined as NC = |C∩C real | |C real | , where C and C real are the learned and ground truth node correspondences respectively.

Matching Permutated Networks
Datasets. Following the literature (eg. , Konar and Sidiropoulos, 2020, Heimann et al., 2021), we match two randomly permutated graphs, that is, the target graph is obtained by permutating the nodes of the original graph (source graph). SIGMA is tested against baselines on standard benchmark datasets, including PPI Yeast [Breitkreutz et al., 2007], Arxiv [Leskovec et al., 2007], and LastFM ASIA [Rozemberczki and Sarkar, 2020]. These datasets are listed in Table 1. Similar to , Konar and Sidiropoulos [2020], Heimann et al. [2021], we also add structural noise to the target graph by adding q|E s | edges.
Results. We report the performance of SIGMA and baselines in Figure 2, in which K is set to 3. The similarity matrix in GRAMPA incurs quadruple computational complexity and takes thousands of hours to obtain on LastFM ASIA and Arxiv. GRAMPA thus does not scale and its performance is not reported on these two datasets. OT-based methods and RefiNA take into account the local topology and have good performance on all datasets. SIGMA outperforms state-of-the-art methods in all cases. The improvement becomes more significant with the noise increasing.
Runtime. In Table 2, we compare the average runtimes of all different methods across noise levels. Time used to calculate the heat diffusion wavelets in SIGMA is included. The two fastest methods per dataset are in bold. SIGMA is faster than its closest competitors in accuracy (Figure 2).

Related Work
Matching attributed graphs. Deep learning-based methods take as input the node/edge attributes of the graphs and learn node embeddings suitable for graph matching [Zanfir and Sminchisescu, 2018, Fey et al., 2019, Wang et al., 2019a, Sarlin et al., 2020, Wang et al., 2020b, Yan et al., 2020. These supervised methods, however, require a large amount of anchor links that are often not available in practice. Optimal transport-based methods exploit structural information and match graphs in an unsupervised manner [Maretic et al., 2019, Xu et al., 2019b, Chowdhury and Mémoli, 2019. To further incorporate node attributes, Fused Gromov-Wasserstein distance-based methods are proposed to exploit both node attributes and structure information [Titouan et al., 2019, Barbe et al., 2020]. Node or edge attributes can be incorporated into the design of SI.
Graph distance. Comparison among graphs is ubiquitous in analyzing graphstructured data. Spectral distances usually do not take into account all the structural information, focusing only on the Laplacian matrix eigenvectors and ignoring a large portion of the structure encoded in eigenvectors [Jovanović andStanić, 2012, Gera et al., 2018]. The cut distance [Lovász, 2012] and the graph edit distance [Bougleux et al., 2017, Raveaux, 2021 require solving difficult discrete optimization problems. Another family of graph distances that is closely related to this paper is the graph diffusion distance [Hammond et al., 2013, Tsitsulin et al., 2018, Scott and Mjolsness, 2021. Hammond et al. 2013 assume the nodes of the two graphs are already matched and calculate the graph diffusion distance. Scott and Mjolsness 2021 extend Hammond et al. [2013]'s work and address the graph matching problem. Their method involves solving a three-level optimization problem and estimates the node correspondence in the innermost problem by matching eigenvalues of the heat diffusion kernel matrix. Note that our Eq. (8) also defines a distance between two graphs.

Conclusion
In this paper, we consider the problem of matching unattributed graphs without anchor links. A novel quantity, SI, is defined to measure the correspondence quality of a pair of nodes. We reveal that SI has a close relationship to matching quality. On the one hand, if the SI value of a pair of nodes is zero, they are either true counterparts or have the same local topology. On the other hand, true matching entails low SI values of matched node pairs. Based on SI, we propose a graph matching method named SIGMA which can reduce low SI values under suitable assumptions. Empirical results demonstrate the good performance of our method.

A.2 Missing Proofs of Sec. 5
Proof of Lemma 3. The mirror descent recurrence is T (τ +1) = Proj Π(1,1) T (τ ) The projection is achieved by Sinkhorn scaling [?], that is, where Y = T (τ ) G(T (τ ) ), and u, v ∈ R V + . The recurrence implies that, The feasible domain requires that where the first inequality is due to the fact that u i ≥ 0 and Y ii ≥ 0. Combining Eq. (14), (13), and 1 Proof of Theorem 4. By the definition of the structural inconsistency, Where the second equality is due to the fact that N s ii for all i and i = i * , The right hand side of the inequality is exactly SI(i, i * ; T (τ ) ).

Proof of Theorem 5
We can rewrite the transport plan T as a vector variable x ∈ R V 2 + and the Gromov-Wasserstein distance Eq. (8) as objective f (x) = x Ax, respectively.
It is easy to verify that 1. The feasible domain is a closed convex set.
2. f (x) is a weakly convex function.
3. The norm of the gradient ∇f (x) for all feasible x is bounded.
4. The optimal objective value, denoted as f min , exists and f min > −∞.
By Corollary 3.1 of Zhang and He [2018], SIGMA converges to an -stationary point with the number of mirror descent iterations O( 1 2 ).

B.1 Additional Experiment Details
We first include more experimental details about parameter choice in this subsection. For experiments in the main paper, we set the propagation time t = 1×10 −3 and the number of hops K = 3.
The influence of K. We empirically study the influence of K. On the one hand, increasing K can increase the distinguishability of nodes and thus boost the matching accuracy. On the other hand, SIGMA with a larger K is easier to be affected by structural noise, as is implied by the upper bound in Theorem 2. Figure 3 demonstrates that K = 3 can achieve a balance. With the amount of noise increasing, SIGMA with a larger K deteriorates faster. SIGMA with K = 1 outperforms RefiNA, which demonstrates that SIGMA is more related to MNC theoretically.  Sensitivity to propagation time t. The range of propagation time is determined by line search, which is demonstrated in Figure 4. On all datasets, SIGMA is not sensitive to the choice of propagation time t and achieves good performance when 1 × 10 −6 ≤ t ≤ 1 × 10 −3 . Note that when t is too small, Ψ I − tL where I is the identity matrix and SIGMA is reduced to only considering one-hop information. Therefore, the accuracy slightly drops when t ≤ 1 × 10 −5 .  We further demonstrate the performance of SIGMA by matching real-world correlated graphs, i.e., the structural difference between the two graphs is induced by real noise instead of simulated noise. The experiments are conducted on graphs of Autonomous Systems peering information inferred from Oregon Route-Views (RV) 1 between March 31 st , 2001 and April 28 th , 2001. The University's RV project 2 was a tool for Internet operators to obtain real-time border gateway protocol information about the global routing system. We match RV on March 31 st to RVs between March 31 st and April 28 th . The dataset statistics are listed in the Table 3. We still set t = 1 × 10 −3 and K = 3 for SIGMA. For fair comparison, RefiNA, MM, GWL, and SIGMA use the ouput of CONE-Align to initialize the matching matrix.

B.2 Matching Real-world Related networks
We report the results on Table 4. As RV changes over time, the performance of all methods decreases except the one on April 28 th . This is possibly due to monthly patterns in RVs. SIGMA outperforms all baselines. GDD does not take as input the result of CONE-Align and is a deterministic algorithm. The standard deviation is thus 0. Other 0.000 deviation values are induced by rounding to three decimal places.