In dpmcsuss/iGraphMatch: Tools for Graph Matching
Abstract
iGraphMatch is an R package developed for matching the corresponding vertices between two edge-correlated graphs.
The package covers three categories of prevalent graph matching algorithms including relaxation-based, percolation-based, and spectral-based, which are applicable to matching graphs under the most general settings: weighted directed graphs of different order and graphs of multiple layers with auxiliary graph matching techniques.
We provide versatile options to incorporate prior information in the form of seeds with or without noise and similarity scores.
We also implement an S4 class that overloads popular operations for matrices in R for efficient computation of sparse plus low-rank matrices, which is a common structure we can decompose the matrices into in the process of matching graphs.
In addition, iGraphMatch provides functions to summarize the graph matching results in terms of several evaluation measures and visualize the matching performance.
Finally, the package also enables the users to sample correlated random graph pairs from classic random graph models to generate data for simulations.
This paper illustrates the practical applications of the package to the analysis of graph matching by detailed examples using real data from social networks and bioinformatics.
Introduction {#sec:intro}
The graph matching (GM) problem seeks to find an alignment between the vertex sets of graphs that best preserves common structure across graphs.
This is often posed as minimizing edge disagreements of two graphs over all alignments.
Formally, given $A$ and $B$, two adjacency matrices corresponding to two graphs $G_1=(V_1, E_1)$ and $G_2=(V_2, E_2)$, the goal is to find
\begin{align}
\argmin_{P\in\Pi}\lVert A-PBP^T \rVert_F^2
\end{align}
where $\Pi$ is the set of all permutation matrices.
GM has wide applications in diverse fields, such as pattern recognition (@PatternRec1, @PatternRec2, @PatternRec3), machine learning (@ML1, @ML2), bioinformatics (@bio3, @bio2), neuroscience (@neuro), social network analysis (@SocialNetwork), and knowledge graph queries (@Hu2018-hd).
More generally, the problem of discovering some true latent alignment between two networks can often be posed as variations on the above problem by adjusting the objective function for the setting.
The well-known graph isomorphism problem is a special case of GM problem when there exists a bijection between the nodes of two graphs which exactly preserves the edge structure.
In terms of computational complexity, GM is equivalent to the NP-hard quadratic assignment problem, which is considered a very challenging problem where few theoretical guarantees exist, even in special cases (@QAP).
For certain problems where the graphs are nearly isomorphic, polynomial-time algorithms do exist (@friendly, @Umeyama) but these methods frequently break down for more challenging instances.
This paper presents the detailed functionality of the iGraphMatch R package which serves as a practical tool for the use of prevalent graph matching methodologies.
These algorithms utilize either the spectral embedding of vertices (@Umeyama), or relaxations of the objective function (@PATH, @FAQ), or apply ideas from percolation theory (@Percolation, @ExpandWhenStuck).
The iGraphMatch package provides versatile options of working with graphs in the form of matrices, igraph objects or lists of either, and matching graphs under a generalized setting: weighted, directed, graphs of a different order, and multilayer graphs.
In addition, the iGraphMatch package incorporates prior information: seeds and similarities for all the implemented algorithms.
Seeds, or anchors, refer to partial knowledge of the alignment of two graphs.
In practice, seeds can be users with the same name and location across different social networks or pairs of genes with the same DNA sequences.
Some algorithms like the percolation algorithm (@Percolation, @ExpandWhenStuck) which matches two graphs by propagating matching information to neighboring pairs require seeds to kick off.
All algorithms improve substantially by incorporating seeds and can achieve accurate matching in polynomial time (@SGM).
Similarity scores are another commonly used prior which measures the similarity between pairs of nodes across the graphs.
In the bioinformatics area, BLAST similarity score is an example of similarity scores that plays an important role in aligning two PPI networks (@IsoRank).
Similarity scores are usually generated from nodal covariates that are observed in both networks (@BLAST, @SimilarityScore).
While under many scenarios the availability of exact partial matches, or hard seeding, is not realistic and expensive, the package also enables utilizing noisy prior information.
Similarity scores incorporate uncertainty by assigning the pair of nodes with higher similarity scores a bigger chance to match.
Seeds with uncertainty and even error can still be handled by self-correcting graph matching algorithms like the Frank-Wolfe algorithm initialized at the noisy partial matching, called soft seeding.
@soft_seeding showed that the Frank-Wolfe algorithm with soft seeding scheme converges quickly to the true alignment under the correlated Erd\H{o}s-R\'enyi model with high probability.
Thus, the original intractable problem is reduced to be solvable in polynomial time.
Although there exist some open source software and packages containing graph matching functionality, iGraphMatch package provides a centralized repository for common graph matching methodologies with flexibility, tools for developing graph matching problem methodology, as well as metrics for evaluating and tools for visualizing matching performance.
Among the alternative GM packages, the most relevant ones include
the igraph (@igraph) package which focuses on descriptive network analysis and graph visualization based on igraph objects and provides a single graph matching algorithm,
the GraphM (@PATH) package which implements several GM algorithms proposed between 1999 and 2009 in C, and the
Corbi (@Corbi) R package which is particularly designed for studies in bioinformatics and SpecMatch (@SpecMatch) which only involves implementations of spectral embedding based GM algorithms and written in C/C++.
None of these packages provide the breadth of tools, flexibility, and ease-of-use provided by the iGraphMatch package.
The rest of this paper is organized as follows.
Section \@ref(sec:background) describes the theoretical representations of the implemented GM algorithms, correlated random graph models and evaluation metrics.
Section \@ref(sec:usage) discusses the functionality and usage of R functions in the package, illustrated on the synthetic correlated graph pairs.
Section \@ref(sec:example) presents more complex examples on real data with several functions involved in the analysis and section \@ref(sec:conclusion) gives guidelines on using different GM algorithms under different circumstances and concludes the paper.
Graph matching background {#sec:background}
In this section, we give background on graph matching and related problems followed by descriptions of the principal algorithms implemented in iGraphMatch.
For simplicity, we state all the algorithms in the context of matching undirected, unweighted graphs with the same cardinality.
All algorithms can also be directly applied to directed and weighted graphs.
In the second subsection, we discuss the techniques for matching graphs with a different number of vertices along with other extensions.
To conclude the section, we introduce the statistical models for correlated networks and discuss measures for the goodness of matching.
For the remainder of this paper we use the following notation.
Let $G_1=(V_1,E_1)$ and $G_2=(V_2,E_2)$ denote two graphs with $n$ vertices.
Let $A$ and $B$ be their corresponding binary symmetric adjacency matrices.
In the setting of seeded graph matching, suppose without loss of generality, the first $s$ pairs of nodes are seeds for simplicity.
In iGraphMatch, much more flexible seed specifications are possible, which will be illustrated in examples for usage of the package in section \@ref(sec:usage).
Accordingly, let $A$ and $B$ be partitioned as:
\begin{equation} \label{eq:seed_blocks}
A=
\begin{bmatrix}
A_{11} & A_{21}^T \
A_{21} & A_{22}
\end{bmatrix}\text{ and }
B=
\begin{bmatrix}
B_{11} & B_{21}^T \
B_{21} & B_{22}
\end{bmatrix}
\end{equation}
where $A_{11}, B_{11}\in{0, 1}^{s\times s}$ denote seed-to-seed adjacencies, $A_{21}, B_{21}\in{0, 1}^{(n-s)\times s}$ denote nonseed-to-seed adjacencies and $A_{22}, B_{22}\in{0, 1}^{(n-s)\times (n-s)}$ denote nonseed-to-nonseed adjacencies.
Let $S$ be an $n$-by-$n$ real-valued matrix of similarity scores.
Let $\Pi$ be the set of all permutation matrices and $\mathcal{D}$ be the set of all doubly stochastic matrices.
Assignment problems
Matching or assignment problems are core problems in combinatorial optimization and appear in numerous fields (@AP).
As we illustrate in Eq. \eqref{eq:ob_func}, a general version of the graph matching problem is equivalent to the quadratic assignment problem (QAP).
Similarly, QAP is related to the linear assignment problem (LAP) which also plays a role in GM.
The LAP asks how to assign $n$ items (eg. workers or nodes in $G_1$) to $n$ other items (eg. tasks or nodes in $G_2$) with minimum cost.
Let $C$ denote an ${n\times n}$ cost matrix, where $C_{ij}$ denotes the cost of matching $i$ to $j$, then the LAP is to find
\begin{equation} \label{eq:LAP}
\begin{aligned}
& \argmin_{P\in \Pi}
& \mathrm{trace}(C^TP)
\end{aligned}
\end{equation}
LAP is solvable in $O(n^3)$ time and there are numerous exact and approximate methods for both general (@lap_solver, @Hungarian) and special cases, such as sparse cost matrices (@Volgenant1996-o).
The statement of QAP resembles LAP, except that the cost function is expressed as a quadratic function.
Given two $n\text{-by-}n$ matrices $A$ and $B$ which can represent flows between facilities and the distance between locations respectively, or the adjacency matrices of two unaligned graphs, the objective function for QAP is:
\begin{equation}\label{eq:qap}
\argmin_{P\in \Pi} \mathrm{trace}(APBP^T).
\end{equation}
This problem is NP-hard (@QAP) leading to a core challenge for any graph matching approach.
As will be illustrated in the rest of the section, some matching algorithms reduce the graph matching problem to solving a LAP.
For these algorithms, we include similarity scores $S$ by adding an additional term $\mathrm{trace}(S^TP)$ to the reduced objective function.
Graph matching algorithms {#subsec:gm}
In the iGraphMatch package, we implement three types of prevalent GM algorithms.
The first group uses relaxations of the objective function, including convex, concave, and indefinite relaxations.
The second group consists of algorithms that apply ideas from percolation theory, where matching information is spread from an initial set of matched nodes.
The last group is based on the spectral embedding of vertices.
Relaxation-based algorithms
These approaches relax the constraint that $P$ is a permutation matrix to require only that $P$ is doubly stochastic, optimizing over $\mathcal{D}$, the convex hull of $\Pi$.
When $P$ is a permutation matrix
\begin{equation}
\label{eq:ob_func}
\lVert A-PBP^T \rVert_F^2 = \lVert AP-PB \rVert_F^2 = \lVert A \rVert_F^2 + \lVert B \rVert_F^2 - 2\cdot \mathrm{trace}APBP^T.
\end{equation}
However, these equalities do not hold for all $P\in \mathcal{D}$, leading to different relaxations.
The second term of Eq.\eqref{eq:ob_func} is a convex function and optimizing it over $P\in\mathcal{D}$ gives the convex relaxation, where the gradient at $P$ to the convex relaxed objective function is $-4 APB + 2A^TAP + 2PBB^T$.
The last equality in Eq. \eqref{eq:ob_func} shows that minimizing edge disagreements is equivalent to maximizing the number of edge agreements, $\mathrm{trace} APBP^T$, a QAP.
Optimizing the indefinite function over $\mathcal{D}$ gives the indefinite relaxation with gradient $-2APB$ (@FAQ).
| Relaxation | Objective Function | Domain | GM Algorithm | Optimization Guarantee | Optimum Form |
| ---------- | ---------------- | ------------- | ------------ | ---------------- | -------------- |
| None | $\lVert A-PBP^T \rVert_F^2$ | $\Pi$ | NA | | $\Pi$ |
| Indefinite | $\mathrm{tr}ADBD^T$ | $\mathcal{D}$ | FW | Local | $\mathcal{D}$ (often $\Pi$) |
| Convex | $\lVert AD-DB \rVert_F^2$ | $\mathcal{D}$ | FW, PATH | Global | $\mathcal{D}$ |
| Concave | $-\mathrm{tr}(\Delta D)-2\mathrm{tr}(L_1^T D L_2 D^T)$ | $\mathcal{D}$ | PATH | Local | $\Pi$ |
Table: (#tab:relaxation) Summary of relaxation methods for graph matching problem.
Generally, the convex relaxation leads to a solution that is not guaranteed to be near the solution to the original GM.
However, @friendly introduced the class of "friendly" graphs based on the spectral properties of the adjacency matrices to characterize the applicability of the convex relaxation.
Matching two friendly graphs by using the convex relaxation is guaranteed to find the exact solution to the GM problem.
Unfortunately, this class is quite limiting and does not hold for most statistical models or real-world examples.
Another relaxation is the concave relaxation used in the PATH algorithm (@PATH).
The concave relaxation uses the Laplacian matrix defined as $L=D-A$, where $D$ is the diagonal degree matrix with diagonal entries $D_{ii}=\sum_{i=1}^N A_{ij}$.
Assume $L_i$ and $D_i$, $i=1,2$, are the Laplacian matrices and degree matrices for $G_1$ and $G_2$ respectively, then we can rewrite the objective function as
\begin{equation} \label{eq:concave}
\begin{split}
\lVert A-PBP^T\rVert_F^2 & = \lVert AP-PB \rVert_F^2\
& = \lVert (D_1P-PD_2)-(L_1P-PL_2)\rVert_F^2\
& = -\mathrm{trace}(\Delta P)+\mathrm{trace}(L_1^2)+\mathrm{trace}(L_2^2)-2\mathrm{trace}(L_1^T P L_2 P^T),
\end{split}
\end{equation}
where the matrix $\Delta_{ij}=(D_{2_{jj}}-D_{1_{ii}})^2$.
Dropping the terms not dependent on $P$ in equation \ref{eq:concave}, we obtain the concave function $-\mathrm{trace}(\Delta P)-2\mathrm{trace}(L_1^T P L_2 P^T)$ on $\mathcal{D}$.
A summary of the different relaxations is provided in Table \ref{tab:relaxation}.
Relaxing the discrete problem to a continuous problem breaks the equivalence to the original formulation of the edge disagreement and enables employing algorithms based on gradient descent.
Frank Wolfe methodology
@FAQ introduced an algorithm for the relaxed graph matching problem, with each iteration computable in polynomial time, that can find local optima for the relaxations above.
The \Frank-Wolfe (FW) (@FW) methodology is an iterative gradient ascent approach composed of two steps.
The first step finds an ascent direction that maximizes the gradient ascent.
In this case the ascent direction is a permutation matrix which is a vertex of the polytope of doubly stochastic matrices.
For the convex, indefinite, and concave relaxations, this corresponds to a LAP with the gradient as the cost function.
The second step performs a line search along the ascent direction to optimize the relaxed objective function.
As the objectives are all quadratic, this line search is simply requires optimizing a single-variable quadratic function along a line segment.
After the iterative algorithm converges, the final step of the procedure is to project the doubly stochastic matrix back to the set of permutation matrices, which is also a LAP.
The various relaxed forms can all serve as the objective function $f(\cdot)$ in the FW Methodology, but in all cases a matrix $D^0\in \mathcal{D}$ must be chosen to initialize the procedure.
For the convex relaxation, the FW methodology is guaranteed to converge to the global optimum regardless of the $D^0$.
On the other hand, the FW algorithm for the indefinite relaxation is not guaranteed to find a global optimum so the initialization is critical.
In many instances, the optimal solution to the convex relaxation lies in the interior of $\mathcal{D}$.
This can lead to inaccurate solutions after the last projection step.
The local optima for the indefinite relaxation are often at extreme points of $\mathcal{D}$, meaning the final projection often does nothing.
The default initialization for the indefinite problem is at the barycenter matrix, $D^0 = \frac{1}{n}11^T$, but many other initialization procedures can be used.
These include randomized initializations, initializations based on similarity matrices, and initializing the indefinite relaxation at the interior point solution of the convex relaxation (@relax_paper).
When prior information regarding a partial correspondence is known to be noisy, rather than incorporating this information as seeds, one can incorporate it as "soft" seeds which are used to generate the initialization (@soft_seeding).
\begin{algorithm}
\caption{Frank-Wolfe Methodology}
\label{alg:FAQ}
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{$A,B$, doubly stochastic matrix $D^0$, tolerance $\epsilon$}
\Output{permutation matrix $P$}
Set $i=0$;\par
\While{$\Vert D^i-D^{i-1} \Vert^2_F\ge \epsilon$}{
$P^i=$ $\argmin_{P\in\Pi}\mathrm{trace}\nabla f(D^i)^TP$;\par
$D^{i+1}=$ $\argmin_{D\in\mathcal{D}}f(D)$ over line segment from $D^i$ to $P^i$;\par
$i=i+1$;\par
}
Project $D^{i}$ to the nearest $P$ by maximizing $\mathrm{trace}(P^TD)$;\par
\end{algorithm}
When prior information is available in the form of seeds,
the seeded graph matching problem (@SGM) works on the objective function \ref{eq:ob_func} with the permutation matrix $P^{n\times n}$ substituted by $I_s\oplus P^{(n-s)\times (n-s)}$, the direct sum of an $s\times s$ identity matrix and an $(n-s)\times (n-s)$ permutation matrix.
Employing the indefinite relaxed objective function incorporating seeds, we formulate the problem as finding
\begin{align}
\hat{P} % &= \argmin_{P\in \Pi} \mathrm{trace} A(I_s\oplus P)B(I_s\oplus P)^T \
&= \argmax_{P\in\mathcal{D}} 2\cdot\mathrm{trace}P^TA_{21}B_{21}^T+\mathrm{trace}A_{22}PB_{22}P^T
\end{align}
where the gradient to the objective function is
\begin{equation} \label{eq:sgm_gradient}
\nabla f(P)=2\cdot A_{21}B_{21}^T+2\cdot A_{22}PB_{22}.
\end{equation}
In total, this uses the information between seeded nodes and nonseeded nodes and the nonseed-to-nonseed information.
Applying seeded graph matching to the convex relaxation and concave relaxation closely resembles the case of indefinite relaxation.
PATH algorithm
@PATH introduced a convex-concave programming approach to approximately solve the graph matching problem.
The concave relaxation has the same solution as the original graph matching problem.
The PATH algorithm finds a local optimum to the concave relaxation by considering convex combinations of the convex relaxation $F_0(P)$ and the concave relaxation $F_1(P)$ denoted by $F_{\lambda} =(1 - \lambda) F_0 + \lambda F_1$.
Starting from the solution to the convex relaxation ($\lambda=0$) the algorithm iteratively performs gradient ascent using the FW methodology at $F_\lambda$, increasing $\lambda$ after each iteration, until $\lambda = 1$.
Percolation-based algorithms
Under the FW methodology, all the nodes admit a correspondence but the (relaxed) matching correspondence evolves through iterations.
On the other hand, percolation approaches start with a set of seeds, adding one new match at each iteration.
The new matches are fixed and hence not updated in future iterations.
Each iteration expands the set of matched nodes by propagating the current matching information to neighbors.
The guiding intuition is that more matched neighbors are an indicator of a more plausible match, an intuition analogous to the gradient ascent approaches above.
We will present two algorithms in this category where the ExpandWhenStuck algorithm is an extension to the Percolation algorithm.
There are some distinctions about the inputs and outputs of percolation methods compared to the above relaxation methods.
Percolation Algorithm
@Percolation provides a simple and fast approach to solve the graph matching problem by starting with a handful of seeds and propagating to the rest of the graphs.
At each iteration, the matching information up to the current iteration is encoded in a subpermutation matrix $P$ where $P_{ij}=1$ if $i$ is matched to $j$, and $0$ otherwise.
The Percolation algorithm searches for the most promising new match among the unmatched pairs through the mark matrix, $M=APB$,
which is the gradient of the indefinite relaxation when extended to sub-doubly stochastic matrices.
When similarity scores are available, they are added to the mark matrix to combine topological structure and similarity scores.
Adopting analogous partitions on the adjacency matrices as in equation \eqref{eq:seed_blocks}, we let $A_{21}, B_{21}$ denote sub-matrix corresponding to potential adjacencies between unmatched and matched nodes.
Since all the candidates of matched pairs are permanently removed from consideration, we need only consider $M'=A_{21}B_{21}^T$, the sub-matrix of $M$ corresponding to the unmatched nodes in both graphs.
As a result, the Percolation algorithm only uses matched-to-unmatched information to generate new matches.
Moreover, the mark matrix $M$ can also be interpreted as encoding the number of matched neighboring pairs for each pair of nodes $i\in V_1$, $j\in V_2$.
Suppose $u, u'\in V_1$, $v,v'\in V_2$, $[u,u']\in E_1$ and $[v,v']\in E_2$, then $(u',v')$ is a neighboring pair of $(u,v)$.
In each iteration, while there remain unmatched nodes with more than $r$ matched neighboring pairs, the percolation algorithm matches the pair of nodes with the highest score $M_{uv}$,
and adds one mark to all the neighboring pairs of $(u,v)$.
Note that the algorithm may stop before all nodes are matched, leading to the return of a partial match.
There is only one tuning parameter in the Percolation algorithm, the threshold $r$ which controls a tradeoff between quantity of matches and quality of matches.
With a small threshold, the algorithm has a larger chance of matching wrong pairs.
If $r$ is larger, then the algorithm might stop before matching many pairs (@ExpandWhenStuck).
The Percolation algorithm can be generalized to matching weighted graphs by making an adjustment to how we measure the matching information from the neighbors.
Since we prefer to match edges with smaller weight differences and higher absolute weights, we propose to adopt the following update formula for the score associated with each pair of nodes $(i,j)$:
$$M_{ij}=M_{ij} + \sum_{u\in N(i)}\sum_{v\in N(j)}1-\frac{|w_{iu}-w_{jv}|}{\max(|w_{iu}|, |w_{jv}|)}.$$
Thus, the score contributed by each neighboring pair of $(i,j)$ is a number in $[0,1]$.
\begin{algorithm}[H]
\SetAlgoLined
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{$A$, $B$, $s$ pairs of seeds, threshold $r$}
\Output{(sub-)permutation matrix $P$}
Initialize the sub-permutation matrix $P$ incorporating seeds;\par
Calculate the mark matrix $M=APB$;\par
Set the rows and columns of $M$ corresponding to seeds to minus infinity;\par
\While{$max(M)\ge r$}{
$P_{ij}\leftarrow 1$, where $[i,j]$ is the index of $max(M)$;\par
$M\leftarrow APB$;\par
Set $i^{th}$ row and $j^{th}$ column of $M$ to minus infinity;\par
}
\caption{Percolation Algorithm}
\label{alg:perco}
\end{algorithm}
ExpandWhenStuck Algorithm
@ExpandWhenStuck extends the Percolation algorithm to a version that can operate with a smaller number of seeds.
Without enough seeds, when there are no more unmatched pairs with a score higher or equal to the threshold $r$, the Percolation algorithm would stop even if there are still unmatched pairs.
ExpandWhenStuck uses all pairs of nodes with at least one matched neighboring pair, $M_{ij}\geq 1$, as new seeds to restart the matching process by adding one mark to all of the new seeds' neighboring pairs, without updating the matched set.
If the updated mark matrix consists of new pairs with marks greater or equal to $r$, then the percolation algorithm continues, leading to larger matched sets.
Spectral-based algorithm
Another class of graph matching algorithms uses the spectral properties of adjacency matrices.
IsoRank Algorithm
@IsoRank propose the IsoRank algorithm that uses neighboring topology and similarity scores and exploits spectral properties of the solution.
The IsoRank algorithm is also based on the relaxation-based algorithms by encoding the topological structure of two graphs in $ADB$, which is again proportional to the gradient of the indefinite relaxation.
However, the representations of each term of $ADB$ are slightly different.
$A$ and $B$ are the column-wise normalized adjacency matrices and $D$ is not necessarily a doubly stochastic matrix yet $D_{ij}$ still indicates how promising it is to match $i\in V_1$ to $j\in V_2$.
Similar to the idea of Percolation algorithm, the intuition is that the impact of a pair of matched nodes is evenly distributed to all of their neighbors to propagate plausible matches.
This is achieved by solving the eigenvalue problem
\begin{align}\label{eq:Iso1}
\mathrm{vec}(D)=(A\otimes B) \mathrm{vec}(D),
\end{align}
where $\mathrm{vec}(D)$ denotes the vectorization of matrix $D$, and the right hand side is equivalent to $ADB$.
To combine network-topological structure and similarity scores in the objective function, the normalized similarity score $E$ is added to the right hand side of Eq. \ref{eq:Iso1}, where $E=S/\|S\|_1$, and $|\cdot|_1$ denotes the L1 norm.
Note that when similarity score is not available as prior information, we can also construct a doubly stochastic similarity score matrix from seeds by taking $I_{s\times s}\oplus \frac{1}{n-s}11_{(n-s)\times (n-s)}$.
To solve the eigenvalue problem \ref{eq:Iso1}, we resort to the power method.
Finally, the global alignment is generated by a greedy algorithm or using the algorithms for solving the linear assignment problem (LAP).
\begin{algorithm}
\caption{IsoRank Algorithm}
\label{alg:IsoRank}
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{$A,B$, similarity scores $S$}
\Output{permutation matrix $P$}
Column-wise normalize the adjacency matrices $A$ and $B$;\par
Normalize similarity scores: $E=S/|S|_1$;\par
Initialize $D^0=E$ and tolerance $\epsilon$;\par
\While{$| D^i-D^{i-1}|\ge \epsilon$}{
Calculate $D^{i+1}=AD^iB+E$;\par
Normalize $D^{i+1}=D^{i+1}/|D^{i+1}|$;\par
$i=i+1$;\par
}
Extract node mapping from $\hat{D}$ using a greedy method or by solving a LAP;
\end{algorithm}
Umeyama algorithm
@Umeyama is a spectral approach to find approximate solutions to the graph matching problem.
Assuming eigendecompositions of adjacency matrices $A$ and $B$ as $A=U_A\Lambda_AU_A^T$ and $B=U_B\Lambda_BU_B^T$, let $|U_A|$ and $|U_B|$ be matrices which takes absolute values of each element of $U_A$ and $U_B$.
Such modification to the eigenvector matrices guarantees the uniqueness of eigenvector selection.
The global mapping is obtained by minimizing the differences between matched rows of $U_A$ and $U_B$:
\begin{align}
\hat{P}=\mathop{\argmin}{P\in\Pi}\lVert |U_A|-P|U_B|\rVert_F=\mathop{\argmax}{P\in\Pi}\mathrm{trace}(|U_B||U_A|^TP)
\end{align}
The Umeyama algorithm can be generalized to matching directed graphs by eigendecomposing the Hermitian matrices $E_A$ and $E_B$ derived from the asymmetric adjacency matrices of the directed graphs.
The Hermitian matrix for the adjacency matrix $A$ is defined as $E_A=A_S+iA_N$, where $A_S=(A+A^T)/2$ is a symmetric matrix, $A_N=(A-A^T)/2$ is a skew-symmetric matrix and $i$ is the imaginary unit.
Similarly, we can define the Hermitian matrix for $B$.
Assume the eigendecompositions of $E_A$ and $E_B$ as follows:
\begin{align}
E_A=W_A\Gamma_AW_A^, \quad E_B=W_B\Gamma_BW_B^
\end{align}
and we aim at searching for:
\begin{align}
\hat{P}=\mathop{\argmax}_{P\in\Pi}\mathrm{trace}(|W_B||W_A|^TP)
\end{align}
Note that the Umeyama algorithm works on the condition that two graphs are isomorphic or nearly isomorphic.
\begin{algorithm}
\caption{Umeyama Algorithm}
\label{alg:Umeyama}
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\Input{$A,B$}
\Output{permutation matrix $P$}
Compute the eigendecompositions of $A$ and $B$: $A=U_A\Lambda_AU_A^T$, $B=U_B\Lambda_BU_B^T$ ;\par
Solve the LAP: $P=\mathop{\argmax}_{P\in\Pi}\mathrm{trace}(|U_B||U_A|^TP)$;\par
\end{algorithm}
Auxiliary graph matching tools
Centering technique
Instead of encoding the non-adjacencies by zeros in the adjacency matrices, the centering technique (@centering) assigns negative values to such edges.
The first approach is encoding non-adjacent node-pairs as $-1$ with centered adjacency matrices $\tilde{A}=2A-\textbf{J}$ and $\tilde{B}=2B-\textbf{J}$, where $\textbf{J}$ is a matrix of all ones.
An alternative approach relies on modeling assumptions where the pair of graphs are correlated but do not share a global structure.
We match $\tilde{A} = A-\Lambda_A$ and $\tilde{B} = B-\Lambda_B$, where $\Lambda$ is an $n\text{-by-}n$ matrix with $ij$-th entry denoting an estimated marginal probability of an edge.
In general, $\Lambda$ is unknown but there are methods in the literature to estimate $\Lambda$.
Matching centered graphs changes the rewards for matching edges, non-edges, and the penalties for mismatches.
Adapting the centering technique for the problem at hand can be used to find specific types of correspondences.
This can also be combined with constructing multilayer networks out of single layer networks to match according to multiple criteria (@Lin, @GRAMPA).
The centering technique can be applied to any of the implemented graph matching algorithm.
It is especially useful when padding graphs with differing numbers of vertices to distinguish isolated vertices from padded vertices.
Padding graphs of different orders
Until this section, we have been considering matching two graphs whose vertex sets are of the same cardinality.
However, matching graphs with different orders are commonly seen in real-world problems.
Suppose $A\in{0,1}^{n\times n}$ and $B\in{0,1}^{n_c\times n_c}$ with $n_c<n$.
One can then pad the smaller graph with extra vertices to match the order of the larger graph, $\tilde{B}=B\oplus \textbf{0}_{n-n_c}$ and match $A$ and $\tilde{B}$.
Every implemented graph matching algorithm in the iGraphMatch package automatically handles input graphs with a different number of vertices using sparse padding with minimal memory impact.
Since the isolated vertices and the padded vertices share the same topological structure, it can be useful to center the original graphs first then pad the smaller graph in the same manner.
This approach serves to differentiate between isolated vertices the padded ones.
It's theoretically verified that in the correlated \ErdosRenyi graph model, the centered padding scheme is guaranteed to find the true correspondence between the nodes of $G_1$ and the induced subgraph of $G_2$ under mild conditions even if $|V_1|\ll|V_2|$, but the true alignment is not guaranteed without centering (@centering).
Exploiting sparse and low-rank structure
Many real-world graphs, especially large graphs, are often very sparse with $o(n^2)$ and often $\theta(n)$ edges.
This can increase the difficulty of the graph matching problem due to the fact that there are fewer potential edges to match, but sparse graphs also come with computational advantages.
We rely on igraph and Matrix for efficient storage of these structures as well as the efficient implementation of various matrix operations.
We also use the LAPMOD algorithm for sparse LAP problems (@Volgenant1996-o) (see below).
Similarly, a low-rank structure appears in some of the procedures including starting at the rank-1 barycenter matrix and the different centering schemes.
Since low-rank matrices are generally not sparse and visa-versa we implemented the splr S4 class, standing for sparse plus low-rank matrices.
This class inherits from the Matrix class and includes slots for an $n\times n$ sparse matrix x and $n\times d$ dense matrices a and b, to represent matrices of the form x + tcrossprod(a, b).
This class implements efficient methods for matrix multiplication and other operations that exploit the sparse and low-rank structure of the matrices.
Specifically, these methods often require only $O(\|x\|_0) + O(nd)$ storage as opposed to $O(n^2)$ required for densely stored matrices, and enjoy analogous computational advantages.
While users can also use these matrices explicitly, most use of them is automatic within functions such as init_start and center_graph and the matrices can largely be used interchangeably with other matrices.
LAP methods {#sssec:lap_methods}
Multiple graph matching methods include solving an LAP and so we have included multiple methods for solving LAPs into the package.
Specifically we implement the Jonker-Volgenant algorithm (@lap_solver) for dense cost matrices and the LAPMOD algorithm (@Volgenant1996-o) for sparse cost matrices.
Both algorithms are implemented in C to provide improved performance.
The LAPMOD approach is typically advantageous when the number of non-zero entries is less than 50%.
We also depend on the clue package for the solve_LSAP function which implements the Hungarian algorithm (@Papadimitriou1998-lz) for solving an LAP.
Each of these methods can be used independently of a specific graph matching method using the do_lap function.
Multi-layered graph matching
Frequently, networks edges may have categorical attributes and from these categories, we can construct multilayer graphs (@multilayer), where each layer in the networks contains edges from specific categories.
For matching two multilayer graphs, the standard graph matching problem can be extended as $\sum_{l=1}^{m}\|A^{(l)} - PB^{(l)}P^T\|_F^2$ where $m$ denotes the number of categories and $A^{(l)}, B^{(l)}$ are the adjacency matrices for the $l$th layers in each graph.
Note, we assume that the layers are aligned, so that layer $l$ corresponds to the same edge-types in both multi-layer networks.
For an igraph object, the function split_igraph can be used to convert a single object with categorical edge attributes into a list with each element only containing the edges with a specific attribute value.
The implemented algorithms can seamlessly match multi-layer graphs, which are encoded as a list of either igraph objects or matrix-like objects.
We also implemented a matrixlist S4 class that implements many standard matrix operations so that algorithms can be easily extended to work with multilayer networks.
Correlated random graph models {#graph-models}
The correlated \ErdosRenyi model (@FAQ) is essential in the theoretical study of graph matching algorithms.
In a single \ErdosRenyi graph, each edge is present in the graph independently with probability $p$.
The correlated \ErdosRenyi model provides a joint distribution for a pair of graphs, where each graph is marginally distributed as an \ErdosRenyi graph and corresponding edge-pairs are correlated.
\begin{defn}[Correlated \ErdosRenyi]
Suppose $p \in (0,1)$ and $\rho\in [\max{\frac{-p}{1-p}, \frac{p-1}{p}},1]$, a pair of adjacency matrices $(A, B)\sim \text{CorrER}(pJ, \rho J)$ if:
For each $1\le u \le v \le n, A_{uv}$ are independent with $A_{uv} \sim Bernoulli(p)$.\;
For each $1\le u \le v \le n, B_{uv}$ are independent with $B_{uv} \sim Bernoulli(p)$.\;
$A_{uv}$ and $B_{u'v'}$ are independent unless $u=u', v=v'$ in which case the Pearson correlation $\mathrm{corr}(A_{uv},B_{uv})=\rho$.
\end{defn}
To sample a pair of correlated \ErdosRenyi graphs with edge probability $p$, and Pearson correlation $\rho$, we first sample three independent \ErdosRenyi graphs $G_1$, $Z_0$ and $Z_1$ with edge probabilities $p$, $p(1-\rho)$ and $p+\rho(1-p)$ respectively.
Let $G_2 = (Z_1 \cap G_1) \bigcup (Z_0\cap G_1^c)$.
@Percolation provide an alternative formulation for the correlated \ErdosRenyi model.
First, one samples a single random \ErdosRenyi graph $G$ with edge probability $p'$.
Conditioned on $G$, each edge in $G$ is present independently in $G_1,G_2$ with probability $s'$.
These two parameterizations are related to each other by the relationship $s'=p+\rho(1-p)$ and $p'=p/(p+\rho(1-p))$.
The original parameterization is slightly more general because it allows for the possibility of negative correlation.
In addition to homogeneous correlated \ErdosRenyi random graphs, we also implement heterogeneous generalizations of this model.
The stochastic block model (@SBM) and the random dot product graphs (RDPG) model (@rdpg) can both be regarded as extensions of the \ErdosRenyi model.
The stochastic block model is useful to represent the community structure of graphs by dividing the graph into $K$ groups.
Each node is assigned to a group and the probability of edges is determined by the group memberships of the vertex pair.
For the RDPG model, each vertex is assigned a latent position in $\mathbb{R}^d$ and edge probabilities are given by the inner product between the latent positions of the vertex pair.
For both of these models, we can consider correlated graph-pairs where marginally they arise from one of these models and again corresponding edge pairs are correlated.
Measures for goodness of matching {#sec:measure}
The ability to assess the quality of the match when ground truth is unavailable is critical for the usage of the matching approaches.
There are various topological criteria that can be applied to measure the quality of matching results.
At the graph level, the most frequently used structural measures include matching pairs (MP), edge correctness (EC), and the size of the largest common connected subgraph (LCCS) (@measure).
MP counts the number of correctly matched pairs of nodes, thus can only be used when the true alignment is available.
Global counts of common edges (CE) and common non-edges (CNE) can be defined as
\begin{align}
CE=\frac{1}{2}\sum_{i,j}1{A_{ij}=PBP^T_{ij} =1}\quad CNE=\frac{1}{2}\sum_{i,j}1{A_{ij}=PBP^T_{ij} = 0},
\end{align}
along with error counts such extra edges (EE) and missing edges (ME),
\begin{align}
EE=\frac{1}{2}\sum_{i,j}1{A_{ij}=0=1 - PBP^T_{ij}}\quad ME=\frac{1}{2}\sum_{i,j}1{A_{ij}= 1 = 1 - PBP^T_{ij}},
\end{align}
EC measures the percentage of correctly aligned edges, that is the fraction CE$/|E_1|$.
The LCCS denotes the largest subset of aligned vertices such that the corresponding induced subgraphs of each graph are connected.
Matches with a larger LCCS are often preferable to those with many isolated components.
Another group of criteria measures the goodness of matching at the vertex level.
Informally, we aim at testing the hypotheses
$$H_0^{(v)}: \text{the vertex } v \text{ is not matched correctly by } P^, $$
$$H_a^{(v)}: \text{the vertex } v \text{ is matched correctly by } P^$$
for each vertex $v$.
The goal is to address if the permutation matrix found by graph matching algorithm is significantly different from the one sampled from uniformly distributed permutation matrices (@row_perm).
Unfortunately, vertex-level matching criteria have only received limited attention in the literature, however, we include two test statistics to measure fit.
The row difference statistic is the L$1$-norm of the difference between $A$ and $P^B{P^}^T$, namely
\begin{equation}
T_d(v,P^):=\Vert A{v\cdot}-(P^B{P^}^T){v\cdot}\Vert_1.
\end{equation}
Intuitively, a correctly matched vertex $v$ should induce a smaller $T_d(v,P^)$,
which for unweighted graphs corresponds to the number of edge disagreements induced by matching $v$.
Alternatively, the row correlation statistic is defined as
\begin{equation}
T_c(v,P^):=1 - corr(A{v\cdot},(P^B{P^}^T)_{v\cdot}).
\end{equation}
We expect the empirical correlation between the neighborhoods of $v$ in $A$ and $P^B{P^*}^T$ to be larger for a correctly matched vertex.
We employ permutation testing ideas to the raw statistics as a normalization across vertices.
Let us take the row difference statistic for example.
The guiding intuition is that if $v$ is correctly matched, the number of errors induced by $P^$ across the neighborhood of $v$ in $A$ and $B$ (i.e., $T_d(v, P^)$) should be significantly smaller than the number of errors induced by a randomly chosen permutation $P$ (i.e., $T_d(v, P)$).
With this in mind, let $\mathbb{E}_P$ and $\mathrm{Var}_P$ denote the conditional expectation and variance of the raw statistic with $P$ uniformly sampled over all permutation matrices.
The normalization is then given by
\begin{equation}
T_p(v,P^):=\frac{T(v,P^)-\mathbb{E}_PT(v,P)}{\sqrt{Var_PT(v,P)}}
\end{equation}
where $T(v, P)$ can be either of the two test statistics we introduced earlier.
In addition to measuring match quality, these vertex-wise statistics can also serve as a tool to find which vertices have no valid match in another network, i.e. the vertex entity is present in one network but not the other.
R functions and usage {#sec:usage}
The R package iGraphMatch offers versatile options for graph matching and subsequent analysis.
Here we describe the usage of the package including sampling random correlated graph pairs, graph matching, and evaluating matching results.
Sampling correlated random graph pairs
We first illustrate the usage of functions for sampling correlated random graph pairs.
The usage of graph matching will be demonstrated on the graph-pairs sampled using these methods.
Functions of the form sample_correlated_*_pair for sampling random graph pairs have the common syntax:
sample_correlated_*_pair(n, ***model parameters***,
permutation = 1:n, directed = FALSE, loops = FALSE)
The argument n specifies the number of nodes in each graph, and the default options are to sample a pair of undirected graphs without self-loop whose true alignment is the identity.
The permutation argument can be used to permute the vertex labels of the second graph.
The model parameters arguments vary according to different random graph models and typically consist of parameters for marginal graph distributions and for correlations between the corresponding edges.
The functions each return a named list of two igraph objects.
For the homogeneous correlated \ErdosRenyi graph model, the model parameters are p, the global edge probability, and corr, the Pearson correlation between aligned vertex-pairs.
For example, to sample a pair of graphs with 5 nodes from $\mathrm{CorrER}(0.5J, 0.7J)$ we run
knitr::opts_chunk$set(warning=FALSE, message=FALSE)
options(scipen = 1, digits = 2)
library(iGraphMatch)
set.seed(1)
cgnp_pair <- sample_correlated_gnp_pair(n = 5, corr = 0.7, p = 0.5)
(cgnp_g1 <- cgnp_pair$graph1)
cgnp_g1[]
cgnp_g2 <- cgnp_pair$graph2
Since we didn't obscure the vertex correspondence by assigning a value to the permutation argument, the underlying true alignment is the identity.
For the more general heterogeneous correlated \ErdosRenyi graph model, one needs to specify an edge probability matrix and a Pearson correlation matrix.
To sample a pair of graphs from the heterogeneous correlated \ErdosRenyi model again with 5 nodes in each graph, and with random edge probabilities and Pearson correlations:
set.seed(123)
p <- matrix(runif(5^2, .5, .8),5)
c <- matrix(runif(5^2, .5, .8),5)
ieg_pair <- sample_correlated_ieg_pair(n = 5, p_mat = p, c_mat = c)
Since the default is undirected graphs without self-loops, the entries of p and c along and below the diagonal are effectively ignored.
The stochastic block model requires block-to-block edge probabilities stored in the pref.matrix argument and the block.sizes argument indicates the size of each block, along with the Pearson correlation parameter corr.
Next, we sample a pair of graphs from the stochastic block model with two blocks of size 2 nodes and 3 nodes respectively, within-group edge probabilities of .7 and .5, across-group edge probability of .001, and Pearson correlation equal to .5.
pm <- cbind(c(.7, .001), c(.001, .5))
sbm_pair <- sample_correlated_sbm_pair(n = 5, pref.matrix = pm,
block.sizes = c(2,3), corr = 0.5)
The iGraphMatch package also provides functions for sampling a pair of correlated random graphs with junk vertices, i.e vertices that don't have true correspondence in the other graph by specifying the number of overlapping vertices in the argument ncore or overlapping block sizes in the argument core.block.sizes.
The iGraphMatch package offers auxiliary tools for centering graphs to penalize the incorrect matches as well, which is implemented in the center_graph function with syntax:
center_graph(A, scheme = c(-1, 1), use_splr = TRUE)
with the first input being either a matrix-like or igraph object.
The scheme argument specifies the method for centering graphs.
Options include a pair of scalars where the entries of the adjacency matrix are linearly rescaled so that their minimum is min(scheme) and their maximum is max(scheme).
Note, scheme = "center" is the same as scheme = c(-1, 1).
Another option is to pass in a single integer, where the returned value is the adjacency matrix minus its best rank-scheme approximation.
The last argument use_splr is a boolean indicating whether to return a splrMatrix object.
We use use_splr = FALSE here to better display the matrices but use_splr = TRUE will often result in improved performance, especially for large sparse networks.
Here, we center the sampled graph cgnp_g1 using different schemes:
center_graph(cgnp_g1, scheme = "center", use_splr = FALSE)
center_graph(cgnp_g1, scheme = 2, use_splr = FALSE)
Users can then use the centered graphs as inputs to the implemented graph matching algorithms, which serve to alter rewards and penalties for common edges, common non-edges, missing edges, and extra edges.
Graph matching methods
The graph matching methods share the same basic syntax:
gm(A, B, seeds = NULL, similarity = NULL, method = "indefinite",
***algorithm parameters***)
| Parameter | Type | Description | Functions |
| ------- | --------- | ---------------- | --------- |
| start | Matrix or character | Initialization of the start matrix for iterations. | FW, convex |
| lap_method | Character | Method for solving linear assignment problem. | FW, convex, PATH, IsoRank |
| max_iter | Number | Maximum number of iterations. | FW, convex, PATH, IsoRank |
| tol | Number | Tolerance of edge disagreements. | FW, convex, PATH |
| r | Number | Threshold of neighboring pair scores. | percolation |
| ExpandWhenStuck | Boolean | TRUE if performs ExpandWhenStuck algorithm. | percolation |
Table: (#tab:gm-alg) Overview of arguments for different graph matching functions.
The first two arguments for graph matching algorithms represent two networks which can be matrices, igraph objects, or two lists of either form in the case of multi-layer matching.
The seeds argument contains prior information on the known partial correspondence of two graphs.
It can be a vector of logicals or indices if the seed pairs have the same indices in both graphs.
In general, the seeds argument takes a matrix or a data frame as input with two columns indicating the indices of seeds in the two graphs respectively.
The similarity parameter is for a matrix of similarity scores between the two vertex sets, with larger scores indicating higher similarity.
Notably, one should be careful with the different scales of the graph topological structure and the vertex similarity information in order to properly address the relative importance of each part of the information.
The method argument specifies a graph matching algorithm to use, and one can choose from "indefinite" (default), "convex", "PATH", "percolation", "IsoRank", "Umeyama", or a self-defined graph matching function which enables users to test out their own algorithms while remaining compatible with the package.
If method is a function, it should take at least two networks, seeds and similarity scores as arguments.
Users can also include additional arguments if applicable.
The self-defined graph matching function should return an object of the "graphMatch" class with matching correspondence, sizes of two input graphs, and other matching details.
As an illustrative example, graph_match_rand defines a new graph matching function which matches by randomly permuting the vertex label of the second graph using a random seed rand_seed.
We then apply this self-defined GM method to matching the correlated \ErdosRenyi graphs sampled earlier with a specified random seed:
graph_match_rand <- function(A, B, seeds = NULL,
similarity = NULL, rand_seed){
totv1 <- nrow(A[[1]])
totv2 <- nrow(B[[1]])
nv <- max(totv1, totv2)
set.seed(rand_seed)
corr <- data.frame(corr_A = 1:nv,
corr_B = c(1:nv)[sample(nv)])
graphMatch(
corr = corr,
nnodes = c(totv1, totv2),
detail = list(
rand_seed = rand_seed
)
)
}
match_rand <- gm(cgnp_g1, cgnp_g2,
method = graph_match_rand, rand_seed = 123)
Other arguments vary for different graph matching algorithms with an overview given in Table \ref{tab:gm-alg}.
The start argument for the FW methodology with "indefinite" and "convex" relaxations takes any $nns\text{-by-}nns$ matrix or an initialization method including "bari", "rds" or "convex".
These represent initializing the iterations at a specific matrix, the barycenter, a random doubly stochastic matrix, or the doubly stochastic solution from "convex" method on the same graphs, respectively.
Moreover, sometimes we have access to side information on partial correspondence with uncertainty.
If we still treat such prior information as hard seeds and pass them through the seeds argument for "indefinite" and "convex" methods, incorrect information can yield unsatisfactory matching results.
Instead, we provide the option of soft seeding by incorporating the noisy partial correspondence into the initialization of the start matrix.
The core function used for initializing the start matrix with versatile options is the init_start function.
Suppose the first two pairs of nodes are hard seeds and another pair of incorrect seed $(3,4)$ is soft seeds:
hard_seeds <- 1:5 <= 2
soft_seeds <- data.frame(seed_A = 3, seed_B = 4)
We generate a start matrix incorporating soft seeds initialized at the barycenter:
as.matrix(start_bari <- init_start(start = "bari", nns = 3,
ns = 2, soft_seeds = soft_seeds))
An alternative is to generate a start matrix that is a random doubly stochastic matrix incorporating soft seeds as follow
set.seed(1)
as.matrix(start_rds <- init_start(start = "rds", nns = 3,
ns = 2, soft_seeds = soft_seeds))
Then we can initialize the Frank-Wolfe iterations at any of the start matrix by specifying the start parameter.
When there are no soft seeds, we no longer need to initialize the start matrix by using init_start first. Instead we can directly assign an initialization method to the start argument in the gm function:
match_rds <- gm(cgnp_g1, cgnp_g2, seeds = hard_seeds,
method = "indefinite", start = "rds")
Below use solution from the convex relaxation as the initialization for the indefinite relaxation.
set.seed(123)
match_convex <- gm(cgnp_g1, cgnp_g2, seeds = hard_seeds,
method = "indefinite", start = "convex")
Now let's match the sampled pair of graphs \code{sbm_pair} from the stochastic block model by using Percolation algorithm.
Apart from the common arguments for all the graph matching algorithms, Percolation has another argument \code{r} representing the minimum number of matched neighbors required for matching a new qualified vertex pair.
Here we adopt the default value which is 2.
Also, at least one of similarity scores and seeds is required for Percolation algorithm to kick off.
Let's utilize the same set of hard seeds and assume there is no available prior information on similarity scores.
sbm_g1 <- sbm_pair$graph1
sbm_g2 <- sbm_pair$graph2
match_perco <- gm(sbm_g1, sbm_g1, seeds = hard_seeds,
method = "percolation", r = 2)
match_perco
Without enough prior information on partial correspondence, Percolation couldn't find any qualifying matches.
Suppose in addition to the current pair of sampled graphs, the above sampled correlated homogeneous and heterogeneous \ErdosRenyi graphs are different layers of connectivity for the same set of vertices.
We can then match the nonseed vertices based on the topological information in all of these three graph layers.
To be consistent, let's still use the Percolation algorithm with threshold \code{r} equal to 2 and the same set of seeds.
matrix_lA <- list(sbm_g1, ieg_pair$graph1, cgnp_g1)
matrix_lB <- list(sbm_g2, ieg_pair$graph2, cgnp_g2)
match_perco_list <- gm(A = matrix_lA, B = matrix_lB, seeds = hard_seeds,
method = "percolation", r = 2)
match_perco_list
With the same amount of available prior information, we are now able to match all the nodes correctly.
Finally, we will give an example of matching multi-layers of graphs using IsoRank algorithm.
Unlike the other algorithm, similarity scores are required for IsoRank algorithm.
Without further information, we adopt the barycenter as the similarity matrix here.
set.seed(1)
sim <- as.matrix(init_start(start = "bari", nns = 5,
soft_seeds = hard_seeds))
match_IsoRank <- gm(A = matrix_lA, B = matrix_lB,
seeds = hard_seeds, similarity = sim,
method = "IsoRank", lap_method = "LAP")
Graph matching functions return an object of class "graphMatch" which contains the details of the matching results, including a list of the matching correspondence, a call to the graph matching function and dimensions of the original two graphs.
match_convex@corr
match_convex@call
match_convex@nnodes
Additionally, "graphMatch" also returns a list of matching details corresponding to the specified method.
Table \ref{tab:gm-value} provides an overview of returned values for different graph matching methods.
With the seeds information, one can obtain a node mapping for non-seeds accordingly
match_convex[!match_convex$seeds]
| Parameter | Description | Functions |
| :------- | :----------------: | ---------: |
| seeds | A vector of logicals indicating if the corresponding vertex is a seed | All the functions |
| soft | The functional similarity score matrix with which one can extract more than one matching candidates | FW, convex, PATH, IsoRank, Umeyama |
| lap_method | Choice for solving the LAP. | FW, convex, Umeyama, IsoRank |
| iter | Number of iterations until convergence or reaches the max_iter. | FW, convex, PATH |
| max_iter | Maximum number of replacing matches. | FW, convex |
| match_order | The order of vertices getting matched. | percolation, IsoRank |
Table: (#tab:gm-value) Overview of return values for different graph matching functions.
The "graphMatch" class object can also be flexibly used as a matrix.
In addition to the returned list of matching correspondence, one can obtain the corresponding permutation matrix in the sparse form.
match_convex[]
Notably, multiplicity is applicable to the "graphMatch" object directly without converting to the permutation matrix.
This enables obtaining the permuted second graph, that is $PBP^T$ simply by
match_convex %*% cgnp_g2
Evaluation of goodness of matching
Along with the graph matching methodology, iGraphMatch has many capabilities for evaluating and visualizing the matching performance.
After matching two graphs, the function summary can be used to get a summary of the overall matching result in terms of commonly used measures including the number of matches, the number of correct matches, common edges, missing edges, extra edges and the objective function value.
The edge matching information is stored in a data frame named edge_match_info.
Note that summary outputs the number of correct matches only when the true correspondence is known by specifying the true_label argument with a vector indicating the true correspondence in the second graph.
Applying the summary function on the matching result match_convex with true_label = 1:5, indicating the true correspondence is the identity that provides these summaries.
summary(match_convex, cgnp_g1, cgnp_g2, true_label = 1:5)
Applying the summary function to a multi-layer graph matching result returns edge statistics for each layer.
summary(match_IsoRank, matrix_lA, matrix_lB)
In realistic scenarios, the true correspondence is not available.
As introduced in section \ref{sec:background}, the user can use vertex level statistics to evaluate match performance.
The best_matches function evaluates a vertex-level metric and returns a sorted data.frame of the vertex-matches with the metrics.
The arguments are the two networks, a specific measure to use, the number of top-ranked vertex-matches to output, and the matching correspondence in the second graph if applicable.
As an example here, we apply best_matches to rank the matches from above with the true underlying alignment
best_matches(cgnp_g1, cgnp_g2, match = match_convex,
measure = "row_perm_stat", num = 3,
true_label = 1:igraph::vcount(cgnp_g1))
Note, best_matches uses seed information from the match parameter and only outputs non-seed matches.
Without the true correspondence, true_label would take the default value and the output data frame only contains the first three columns.
To visualize the matches of smaller graphs, the function plot displays edge discrepancies of the two matched graphs by an adjacency matrix or a ball-and-stick plot, depending on the input format of two graphs.
plot(cgnp_g1, cgnp_g2, match_convex)
plot(cgnp_g1[], cgnp_g2[], match_convex)
The plots for visualizing matching performance of match_convex are shown in Figure \@ref(fig:visualization).
Grey edges and pixels indicate common edges, red ones indicate edges only in the second graph.
If they were present, blue pixels and edges represent missing edges that only exist in the first graph.
The corresponding linetypes are solid, short dash, and long dash.
Examples {#sec:example}
In this section, we demonstrate graph matching analysis using iGraphMatch via examples on real datasets, including communication networks, neuronal networks, and transportation networks.
Table \@ref(tab:dataset-overview) presents brief overviews of the first two datasets.
Note that the number of edges doesn't consider weights for weighted graphs, and for directed graphs, an edge from node $i$ to node $j$ and another edge from $j$ to $i$ will be counted as two edges.
Tables \@ref(tab:edge-summary) and \@ref(tab:edge-summary-trans) summarize the edge correspondence between two graphs under the true alignment including the number of common edges, missing edges, and extra edges in two graphs.
In the first Enron email network example, we demonstrate the usage of Frank-Wolfe methodology and how to improve matching performance by using the centering technique and incorporating adaptive seeds.
In the second example using C. Elegans synapses networks, we illustrate how to use soft matching for a challenging graph matching task using Frank-Wolfe methodology, PATH algorithm and IsoRank algorithm.
Finally, we include an example of matching two multi-layer graphs with similarity scores on the Britain transportation networks.
summarize_graph <- function(gp, name) {
nn <- nrow(gp[[1]][])
if (nn != nrow(gp[[2]][])) {
nn <- paste(nn, "/", nrow(gp[[2]][]))
}
ne <- paste(igraph::ecount(gp[[1]]), "/", igraph::ecount(gp[[2]]))
data.frame(
`Dataset` = name,
`Nodes` = as.character(nn),
`Edges` = ne,
`Correlation` = cor(as.numeric(gp[[1]][]), as.numeric(gp[[2]][])),
`Weighted` =
ifelse(igraph::is_weighted(gp[[1]]), "Yes", "No")
,
`Directed` =
ifelse(igraph::is_directed(gp[[1]]), "Yes", "No")
,
`Loop` = paste(
ifelse(any(igraph::is.loop(gp[[1]])), "Yes", "No"),
"/",
ifelse(any(igraph::is.loop(gp[[2]])), "Yes", "No")
)
)
}
knitr::kable(
rbind(
summarize_graph(Enron, "Enron"),
summarize_graph(C.Elegans, "C. Elegans")
),
col.names = c("Dataset", "\\# Nodes", "\\# Edges", "Correlation", "Weighted", "Directed", "Loop"),
booktabs = TRUE,
caption = "Overview of the Enron and C. Elegans graphs.",
escape = FALSE
)
summarize_edge <- function(gp1, gp2, name) {
nn <- nrow(gp1[])
corr <- graphMatch(data.frame(corr_A = 1:nn, corr_B = 1:nn), nn)
edge_info <- summary(corr, gp1, gp2)$edge_match_info
data.frame(
`Dataset` = name,
`Common Edges` = edge_info$common_edges,
`Missing Edges` = edge_info$missing_edges,
`Extra Edges` = edge_info$extra_edges
)
}
knitr::kable(
rbind(
summarize_edge(Enron[[1]], Enron[[2]], "Enron"),
summarize_edge(C.Elegans[[1]], C.Elegans[[2]], "C. Elegans")
),
col.names = c("Dataset", "Common edges", "Missing edges", "Extra edges"),
booktabs = TRUE,
caption = "Edge summary under the true alignments of the Enron and C. Elegans graphs. The columns indicate the number of common edges, missing edges in $G_1$, and extra edges in $G_2$. For weighted graphs, we define a pair of corresponding edges as a common edge as long as they both have positive weights. ",
escape = FALSE
)
Example: Enron Email Network Data {#sec:Enron}
The Enron email network data was originally made public by the Federal Energy Commission during the investigation into the Enron Corporation (@Enron).
Each node of Enron network represents an email address and if there is at least one email sent from one address to another address, a directed edge exists between the corresponding nodes.
The iGraphMatch package includes the Enron email network data in the form of a pair of igraph objects derived from the original data where each graph represents one week of emails between 184 email addresses.
The two networks are unweighted and directed with edge densities around
r round((igraph::edge_density(Enron[[1]]) + igraph::edge_density(Enron[[2]]))/2, 3) in each graph and the empirical correlation between two graphs is r round(cor(as.numeric(Enron[[1]][]), as.numeric(Enron[[2]][])), 2).
First, let's load packages required for the following analysis:
library(igraph)
library(iGraphMatch)
library(purrr)
library(dplyr)
Visualization of Enron networks
We visualize the aligned Enron networks using the \code{plot} function with vertices sorted by a community detection algorithm (@community_detection) and degree.
For detailed interpretations to figure \@ref(fig:Enron-graph), please refer to figure \@ref(fig:visualization).
g <- igraph::as.undirected(Enron[[1]])
com <- igraph::membership(igraph::cluster_fast_greedy(g))
deg <- rowSums(as.matrix(g[]))
ord <- order(max(deg)*com+deg)
plot(Enron[[1]][][ord,ord], Enron[[2]][][ord,ord])
Note that r sum(degree(Enron[[1]]) == 0) and r sum(degree(Enron[[2]]) == 0) out of the total r vcount(Enron[[1]]) nodes are isolated from the other nodes in two graphs respectively, indicating the corresponding employees haven't sent or received emails from other employees.
This adds difficulty to matching since it's impossible to distinguish the isolated nodes based on topological structure alone.
We first keep only the largest connected component of each graph.
vid1 <- which(largest_cc(Enron[[1]])$keep)
vid2 <- which(largest_cc(Enron[[2]])$keep)
vinsct <- intersect(vid1, vid2)
v1 <- setdiff(vid1, vid2)
v2 <- setdiff(vid2, vid1)
A <- Enron[[1]][][c(vinsct, v1), c(vinsct, v1)]
B <- Enron[[2]][][c(vinsct, v2), c(vinsct, v2)]
The sizes of largest connect components of two graphs are r dim(A)[1] and r dim(B)[1], which are different.
We reorder two graphs so that the first r length(vinsct) nodes are aligned and common to both graphs.
Matching largest connected components using FW Algorithm
Let's assume the Enron email communication network from the second week is anonymous, and we aim at
finding an alignment between the email addresses from the first network and the second one to de-anonymize the latter.
Additionally, we want to find the email addresses that are active in both months.
Suppose no prior information on partial alignment is available in this example.
We match the two largest connected components using the FW algorithm with indefinite relaxation since seeds and similarity scores are not mandatory for this method.
Without any prior information, \code{seeds} and \code{similarity} arguments take default values which are \code{NULL}.
For the \code{start} argument, we assign equal likelihood to all the possible matches by initializing at the barycenter.
Since two graphs are of different sizes, the \code{gm} function automatically pads the smaller graph with extra 0's.
set.seed(1)
match_FW <- gm(A = A, B = B, start = "bari", max_iter = 200)
head(match_FW)
Then, we check the summary of matching performance in terms of matched nodes, matched edges and the graph matching objective function.
summary(match_FW, A, B)
emi <- summary(match_FW, A, B)$edge_match_info
In this example, we can evaluate the matching result based on statistics on matched edges.
Without any seeds or similarity scores, around r paste(round(100*emi$common_edges / (emi$common_edges + emi$missing_edges), 0), "%", sep="") of edges are correctly matched.
Centering the larger graph
We can try to improve performance by centering B by assigning -1 to non-edges, so that we penalize edges that are missing in B but present in A.
A_center <- center_graph(A = A, scheme = "naive", use_splr = TRUE)
B_center <- center_graph(A = B, scheme = "center", use_splr = TRUE)
set.seed(1)
match_FW_center <- gm(A = A_center, B = B_center,
start = "bari", max_iter = 200)
summary(match_FW_center, A, B)
From the summary tables, we would prefer matching Enron networks with the application of the centering scheme, since we get more matched common edges, as well as fewer missing edges and extra edges.
Matching with adaptive seeds
Supposing we have no access to ground truth, we use the \code{best_matches} function to measure and rank the vertex-wise matching performance.
Below shows the 6 matches that minimize the row permutation statistic.
bm <- best_matches(A = A, B = B, match = match_FW_center,
measure = "row_perm_stat")
head(bm)
Since seeded graph matching enhances the graph matching performance substantially (@SGM), it may be useful to use some of these best matches as seeds to improve matching results.
Here, we use adaptive seeds, taking the $ns$ best matches and using them as seeds in a second run of the matching algorithm.
The table below displays edge statistics and objective function values for different number of adaptive seeds used.
The second column in the table shows the matching precision of the adaptive seeds based on ground truth.
Incorporating adaptive seeds and repeating the FW matching procedure on centered graphs further improve the matching results, compared with the case without any adaptive seeds when $ns=0$.
The first 40 pairs of matched nodes ranked by \code{best_matches} function are all correctly matched, and this is also when matching is improved the most.
match_w_hard_seeds <- function(ns){
seeds_bm <- head(bm, ns)
precision <- mean(seeds_bm$A_best == seeds_bm$B_best)
match_FW_center_seeds <- gm(A = A_center, B = B_center,
seeds = seeds_bm, similarity = NULL,
start = "bari", max_iter = 100)
edge_info <- summary(match_FW_center_seeds, A, B)$edge_match_info
cbind(ns, precision, edge_info)
}
set.seed(12345)
map_dfr(seq(from = 0, to = 80, by = 20), match_w_hard_seeds)
set.seed(12345)
map_dfr(seq(from = 0, to = 80, by = 20), match_w_hard_seeds) %>%
knitr::kable(col.names = c("ns", "precision","common", "missing", "extra", "fnorm"),
booktabs = TRUE, digits = 2)
As the number of adaptive seeds increases, the precision of adaptive seeds decreases.
Note that if they are treated as hard seeds, incorrect matches will remain in the matched set and might cause a cascade of errors.
An alternative way is to treat the top-ranked matches as soft seeds embedded in the start matrix to handle the uncertainty.
In this way, adaptive seeds not only provide prior information but also evolve over iterations.
The table below shows that the soft seeding approach always outperforms or performs as good as the hard seeding approach regardless of the number of adaptive seeds being used.
match_w_soft_seeds <- function(ns){
seeds_bm <- head(bm, ns)
precision <- mean(seeds_bm$A_best == seeds_bm$B_best)
start_soft <- init_start(start = "bari",
nns = max(dim(A)[1], dim(B)[1]),
soft_seeds = seeds_bm)
match_FW_center_soft_seeds <- gm(A = A_center, B = B_center,
start = start_soft, max_iter = 100)
edge_info <- summary(match_FW_center_soft_seeds, A, B)$edge_match_info
cbind(ns, precision, edge_info)
}
set.seed(12345)
map_dfr(seq(from = 0, to = 80, by = 20), match_w_soft_seeds)
set.seed(12345)
map_dfr(seq(from = 0, to = 80, by = 20), match_w_soft_seeds) %>%
knitr::kable(col.names = c("ns", "precision","common", "missing", "extra", "fnorm"),
booktabs = TRUE, digits = 2)
Core vertices detection
The function \code{best_matches} can also be used to detect core vertices.
Suppose the ground truth is known and that the first 145 vertices are core vertices.
The mean precision of detecting core vertices and junk vertices using \code{best_matches} function is displayed in figure \ref{fig:core}.
A lower rank is a stronger indicator of a core vertex and a higher rank is a stronger indicator of a junk vertex.
Let $r^C_i, 1\le i\le n_c$ and $r^J_j, 1\le j\le n_j$ denote the ranks associated with each core vertex and each junk vertex.
The figure shows the precision of identifying core vertices at each low rank $r$, i.e. $\frac{1}{r}\sum_{i = 1}^{n_c}1_{r^C_i\le r}$,
and the precision of identifying junk vertices at each high rank $r$, i.e. $\frac{1}{r}\sum_{j = 1}^{n_j}1_{r^J_j\ge n_c+n_j-r}$, which are separated by the vertical lines.
nc <- length(vinsct)
nj <- max(length(v1), length(v2))
core_precision <- 1:nc %>% map_dbl(~mean(bm$A_best[1:.x]<=nc))
junk_precision <- 1:nj %>% map_dbl(~mean(bm$A_best[(nc+.x):(nc+nj)]>nc))
Core detection performance is substantially better than chance, as represented by the dotted horizontal lines.
The top r max(which(core_precision == 1)) are all core vertices indicating good overall performance for core identification.
For junk identification, the junk vertices are ranked r nc + nj - which(bm$A_best > nc) + 1 according to which have the lowest score, indicating that some junk vertices are difficult to identify.
```rMean precision for identifying core and junk vertices for the Enron networks by using the row permutation test. The vertical lines separate the performance of identifying core vertices with low ranks from junk vertices with high ranks. The horizontal lines indicate the performance of a random classifier."}
library(ggplot2)
map_df_core <- data.frame(rank = 1:nc, which = "core",
mean_precision = core_precision)
map_df_junk <- data.frame(rank = (nc+1):(nc+nj), which = "junk",
mean_precision = junk_precision)
map_df <- rbind(map_df_core, map_df_junk)
map_df <- map_df %>%
mutate(random = ifelse(which == "junk", nj/(nc+nj), nc/(nc+nj))) %>%
mutate(start=ifelse(which=="core",0,nc-5),
end=ifelse(which=="core",nc+5,nc+nj))
map_df %>%
ggplot(aes(x=rank,y=mean_precision)) +
geom_line(alpha=.5,size=.3) +
geom_segment(aes(y=random,yend=random, x=start,xend=end),
data=map_df,color="black",alpha=1,linetype=2,size=.1)+
geom_vline(aes(xintercept=nc),size=.1,linetype=2)+
ylab("mean precision")+
theme_minimal()+
theme(text=element_text(family="Palatino"),
axis.line=element_line(size=.3,color="black"),
panel.grid.major=element_line(color="grey95",size=0),
panel.grid.minor=element_line(color="grey98",size=0))
## Example: C. Elegans Network Data {#sec:CE}
The _C. Elegans_ networks consist of the chemical synapses network and the electrical synapses network of the roundworm, where each of `r vcount(C.Elegans[[1]])` nodes represents a neuron and each edge represents the intensity of synapse connections between two neurons (@neuro).
Matching the chemical synapses network to the electrical synapses network is essential for understanding how the brain functions.
These networks are quite sparse with edge densities of `r round(igraph::edge_density(C.Elegans[[1]]), 3)` and `r round(igraph::edge_density(C.Elegans[[2]]), 3)` in each graph and the empirical correlation between two graphs is `r round(cor(as.numeric(C.Elegans[[1]][]), as.numeric(C.Elegans[[2]][])), 2)`.
### A challenging task
For simplicity, we made the networks unweighted and undirected for the experiments, and we assume the ground truth is known to be the identity.
```r
C1 <- C.Elegans[[1]][] > 0
C2 <- C.Elegans[[2]][] > 0
plot(C1[], C2[])
match <- gm(C1, C2, start = Matrix::Diagonal(nrow(C1)))
plot(C1[], C2[], match)
nv <- nrow(C1)
id_match <- graphMatch(data.frame(corr_A = 1:nv, corr_B = 1:nv), nv)
i_sum <- summary(id_match, C.Elegans[[1]], C.Elegans[[2]])
m_sum <- summary(match, C.Elegans[[1]], C.Elegans[[2]], id_match)
i_emi <- i_sum$edge_match_info
m_emi <- m_sum$edge_match_info
Matching the C. Elegans networks is a challenging task.
Figures \@ref(fig:C-Elegans-edge) depict the edge discrepancies of two networks under the true alignment and the matching correspondence using FW algorithm initialized at the true alignment.
The alignment found using FW is not the identity with r m_sum$n_true_match out of r nv nodes correctly matched
and improves upon the identity in terms of the number of edge discrepancies.
For the true alignment, there are
r i_emi$common_edges
edge errors and
r round(i_emi$extra_edges + i_emi$missing_edges)
common edges while the alignment yielded by FW initialized at the true correspondence has
r m_emi$common_edges
edge errors and
r round(m_emi$extra_edges + m_emi$missing_edges)
common edges.
Hence, this graph matching object does not have a solution at the true alignment.
One can try to use other objective functions to enhance the matching result, however we do not investigate this here.
Overall, while most performance measures are poor, our results illustrate the spectrum of challenges for graph matching.
Soft matching: MAP\@3
Considering matching C. Elegans graphs is quite challenging, let's assume 20 pairs of vertices are known as seeds, which are chosen at random.
Accordingly, we generate a similarity matrix with 1's corresponding to seeds, and the rest being barycenter.
seeds <- sample(nrow(C1), 20)
sim <- init_start(start = "bari", nns = nrow(C1), soft_seeds = seeds)
In addition to one-on-one matching, we will also conduct soft matching, which is to find three most promising matches to each non-seed vertex.
We achieve the goal of soft matching by finding the top 3 largest values in each row of the doubly stochastic matrix from the last iteration of Frank Wolfe methodology with indefinite relaxation and PATH algorithm, as well as the normalized matrix from the last iteration of the power method for IsoRank algorithm.
To evaluate the matching performance, we will look at both matching precision: $precision=\frac{1}{n_m-s}\sum_{i\in V_m\setminus S}P_{ii}$, and Mean Average Precision \@ 3 (MAP\@ 3):$MAP@3 = \frac{1}{n_m-s}\sum_{i\in V_m\setminus S}1_{{i\in T_i}}$, where $T_i$ is the set of 3 most promising matches to node $i$.
set.seed(123)
m_FW <- gm(A = C1, B = C2, seeds = seeds,
similarity = sim, method = "indefinite",
start = "bari", max_iter = 100)
m_PATH <- gm(A = C1, B = C2, seeds = seeds,
similarity = NULL, method = "PATH",
epsilon = 1, tol = 1e-05)
m_Iso <- gm(A = C1, B = C2, seeds = seeds,
similarity = as.matrix(sim), method = "IsoRank",
max_iter = 50, lap_method = "LAP")
match_eval <- function(match){
precision <- mean(match$corr_A == match$corr_B)
order <- apply(match$soft, MARGIN = 1, FUN = order, decreasing = TRUE)
top3 <- t(order[1:3,]) - 1:ncol(order)
MAP3 <- mean(apply(top3, MARGIN = 1, FUN = function(v){0 %in% v}))
round(data.frame(precision, MAP3),4)
}
sapply(list(m_FW, m_PATH, m_Iso), match_eval) %>%
knitr::kable(col.names = c("Frank Wolfe", "PATH", "IsoRank"),
booktabs = TRUE, digits = 2)
MAP\@ 3 is slightly higher than precision for each method.
Soft matching provides an alternative way of matching by generating a set of promising matching candidates.
Example: Britain Transportation Network {#sec:Transp}
library(Matrix)
tm <- Transportation[[1]]
cm <- Transportation[[2]]
candidate <- Transportation[[3]]
tn <- nrow(tm[[1]])
wn <- nrow(cm[[1]])
mc <- candidate %>%
with(Matrix::sparseMatrix(i = tem, j = wor, x = 1,
dims = c(tn,wn)))
c_count <- rowSums(mc)
To demonstrate matching multi-layer networks-layers, we consider two graphs derived from the Britain Transportation network (@transportation).
The network reflects the transportation connections in the UK, with five layers representing ferry, rail, metro, coach, and bus.
A smaller template graph was constructed based on a random walk starting from a randomly chosen hub node, a node that has connections in all the layers.
The template graph has r tn nodes and r sum(sapply(tm, sum)) connections in total and is an induced subgraph of the original graph.
Additionally, based on filter methods from \citet{Moorman2018-ha}, the authors of that paper also provided a list of candidate matches for each template node, where the true correspondence is guaranteed to be among the candidates.
The number of candidates ranges from r min(c_count) to r max(c_count) at most, with an average of r round(mean(c_count)) candidates for each template vertex.
Thus, we made an induced subgraph from the transportation network with only candidates, which gave us the world graph with r wn vertices and r sum(sapply(cm, sum)) connections.
tm <- Transportation[[1]]
cm <- Transportation[[2]]
candidate <- Transportation[[3]]
Figure \ref{Fig:trans_net} visualizes the transportation connections for the induced subgraphs, where means of transportation are represented by different colors.
Note that all edges in the template are common edges shared by two graphs,
where r paste0(round(100*sum(tm[[1]]) / sum(cm[[1]][1:tn, 1:tn]),1), "%"), r paste0(round(100*sum(tm[[2]]) / sum(cm[[2]][1:tn, 1:tn]),1), "%"), r paste0(round(100*sum(tm[[3]]) / sum(cm[[3]][1:tn, 1:tn]),1), "%"), r paste0(round(100*sum(tm[[4]]) / sum(cm[[4]][1:tn, 1:tn]),1), "%") and r paste0(round(100*sum(tm[[5]]) / sum(cm[[5]][1:tn, 1:tn]),1), "%") of edges in the world graph are in template for each layer.
All graphs are unweighted, directed, and do not have self-loops.
Tables \@ref(tab:edge-summary-trans) further displays an overview and edge summary regarding each layer of the Britain Transportation Network.
A true correspondence exists for each template vertex in the world graph, our goal is to locate each template vertex in the Britain Transportation network by matching two multi-layer graphs with different number of vertices.
summarize_graph_trans <- function(gp1, gp2, name) {
nn <- nrow(gp1)
if (nn != nrow(gp2)) {
nn <- paste(nn, "/", nrow(gp2))
}
ne <- paste(sum(gp1), "/", sum(gp2))
data.frame(
`Layer` = name,
`Nodes` = as.character(nn),
`Edges` = ne,
`Correlation` = cor(as.numeric(gp1), as.numeric(gp2[1:nrow(gp1), 1:nrow(gp1)]))
)
}
summary_table <-
rbind(
summarize_graph_trans(tm[[1]], cm[[1]], "Ferry"),
summarize_graph_trans(tm[[2]], cm[[2]], "Rail"),
summarize_graph_trans(tm[[3]], cm[[3]], "Metro"),
summarize_graph_trans(tm[[4]], cm[[4]], "Coach"),
summarize_graph_trans(tm[[5]], cm[[5]], "Bus")
)
a <- tm
b <- cm %>% purrr::map(~.x[][1:tn, 1:tn])
edge_layer <- summary(graphMatch(data.frame(corr_A = 1:tn, corr_B = 1:tn), as.integer(tn)), a, b)$edge_match_info
# edge_layer$layer <- c("Ferry", "Rail", "Metro", "Coach", "Bus")
knitr::kable(
cbind(summary_table, edge_layer[,2:4]),
col.names = c("Layer", "\\# Nodes", "\\# Edges", "Correlation","Common", "Missing", "Extra"),
booktabs = TRUE,
caption = "Overview of the Britain Transportation Network layers. Correlation is calculted using the template graph and the aligned induced subgraph of the world graph. The final three columns indicate the number of common edges, missing edges, and extra edges in the aligned subgraph of the world graph.",
escape = FALSE
)
# Edge summary of the Britain Transportation Network layers.
```r Visualization of the template graph (left) and the world graph (right) with corresponding vertices, both derived from the Britain Transportation network with five layers: ferry, rail, metro, coach, and bus. Edges represent transportation transactions and each color indicates a different means of transportation from a different layer of network."}
l <- 1:length(tm)
tg_sub <- l %>% map(~ tm[[.x]] %>% graph_from_adjacency_matrix)
wg_sub <- l %>% map(~ cm[[.x]][1:tn,1:tn] %>% graph_from_adjacency_matrix)
for(l in 1:5){
tg_sub[[l]] <- set_edge_attr(tg_sub[[l]], "layer", value = l)
wg_sub[[l]] <- set_edge_attr(wg_sub[[l]], "layer", value = l)
}
tg_colored <- igraph::union(tg_sub[[1]], tg_sub[[2]], tg_sub[[3]], tg_sub[[4]], tg_sub[[5]])
wg_colored <- igraph::union(wg_sub[[1]], wg_sub[[2]], wg_sub[[3]], wg_sub[[4]], wg_sub[[5]])
tg_layer <- rep(0, ecount(tg_colored))
wg_layer <- rep(0, ecount(wg_colored))
for(l in 1:5){
tg_layer <- tg_layer + ifelse(is.na(edge_attr(tg_colored)[[l]]) | tg_layer != 0, 0, edge_attr(tg_colored)[[l]])
wg_layer <- wg_layer + ifelse(is.na(edge_attr(wg_colored)[[l]]) | wg_layer != 0, 0, edge_attr(wg_colored)[[l]])
}
layout <- layout_nicely(wg_colored)
plot(tg_colored, edge.color = factor(tg_layer),
vertex.size = 5, vertex.label.cex = .666,
edge.curved = TRUE,
edge.arrow.size = 0.2, layout = layout)
plot(wg_colored, edge.color = factor(wg_layer),
vertex.size = 5, vertex.label.cex = .666,
edge.curved = TRUE,
edge.arrow.size = 0.2, layout = layout)
```r
start <- mc %>%
rbind(Matrix(0, nrow = wn - tn, ncol = wn))
start[1:tn, ] <- diag(1 / rowSums(start[1:tn, ])) %*% start[1:tn, ]
similarity <- start * 1e5
Based on the candidates, we specify a start matrix that is row-stochastic which can be used for the \code{start} argument in the graph matching function for FW methodology.
For each row node, its value is either zero or the inverse of the number of candidates for that node.
To ensure that template nodes only get matched to candidates, we constructed a similarity score matrix by taking the start matrix $\times 10^5$, so that a high similarity score is assigned to all the template-candidate pairs.
Then we match the template graph with the world graph using Percolation algorithm.
The template graph stored in \code{tm} and world graph \code{cm} are lists of 5 matrices of dimensions 53 and 2075 respectively.
Since we have no information on seeds, we assign \code{NULL} to the \code{seeds} argument, the Percolation algorithm will initialize the mark matrix using prior information in the similarity score matrix.
match <- gm(A = tm, B = cm, similarity = similarity,
method = "percolation", r = 4)
summary(match, tm, cm)
The \code{summary} function outputs edge statistics and objective function values for each layer separately.
To further improve matching performance, one can replicate all the analysis in the first example on Enron dataset, such as using the centering scheme and adaptive seeds.
Finally, one can refer to the match report to compare matching performance and pick the best one.
Conclusions {#sec:conclusion}
In this work, we detail the methods and usage of the R package iGraphMatch for finding and assessing an alignment between the vertex sets of two edge-correlated graphs.
The package implements common steps for the analysis of graph matching: seamless matching of generalized graphs, evaluation of matching performance, and visualization.
For each of the graph matching methodologies, we provide versatile options for the form of input graphs and the specification of available prior information.
Through the discussion in section \ref{sec:example}, we demonstrate the broad functionality and flexibility of the package by analyzing diverse graph matching problems on real data step by step.
The package also provides tools for simulating correlated graphs which can be used in the development and enhancement of graph matching methods.
Methods for graph matching are still under active development.
We plan to include other novel methods as the field continues to develop.
In the short term we are looking to introduce a suite of additional matching methods that have recently been proposed in the literature.
One of the biggest challenges for graph matching is evaluating the quality of a match, especially at the vertex level.
This has received minimal attention in the previous literature.
We provide measures of goodness of matching on the vertex level and demonstrate their effectiveness empirically.
These baseline methods implement a permutation testing framework for assessing matches that can be readily extended to other metrics.
\section*{Acknowledgments}
The primary authors for the package are Vince Lyzinski, Zihuan Qiao, and Daniel Sussman.
Joshua Agterberg, Lujia Wang, and Yixin Kong also provided important contributions.
We also want to thank all of our users, especially Youngser Park, for their feedback and patience as we continue to develop the package.
This work was supported in part by grants from DARPA (FA8750-20-2-1001 and FA8750-18-0035) and from MIT Lincoln Labs.
References
dpmcsuss/iGraphMatch documentation built on Feb. 15, 2024, 3:26 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Abstract
iGraphMatch is an R package developed for matching the corresponding vertices between two edge-correlated graphs.
The package covers three categories of prevalent graph matching algorithms including relaxation-based, percolation-based, and spectral-based, which are applicable to matching graphs under the most general settings: weighted directed graphs of different order and graphs of multiple layers with auxiliary graph matching techniques.
We provide versatile options to incorporate prior information in the form of seeds with or without noise and similarity scores.
We also implement an S4 class that overloads popular operations for matrices in R for efficient computation of sparse plus low-rank matrices, which is a common structure we can decompose the matrices into in the process of matching graphs.
In addition, iGraphMatch provides functions to summarize the graph matching results in terms of several evaluation measures and visualize the matching performance.
Finally, the package also enables the users to sample correlated random graph pairs from classic random graph models to generate data for simulations.
This paper illustrates the practical applications of the package to the analysis of graph matching by detailed examples using real data from social networks and bioinformatics.
Introduction {#sec:intro}
The graph matching (GM) problem seeks to find an alignment between the vertex sets of graphs that best preserves common structure across graphs. This is often posed as minimizing edge disagreements of two graphs over all alignments. Formally, given $A$ and $B$, two adjacency matrices corresponding to two graphs $G_1=(V_1, E_1)$ and $G_2=(V_2, E_2)$, the goal is to find \begin{align} \argmin_{P\in\Pi}\lVert A-PBP^T \rVert_F^2 \end{align} where $\Pi$ is the set of all permutation matrices. GM has wide applications in diverse fields, such as pattern recognition (@PatternRec1, @PatternRec2, @PatternRec3), machine learning (@ML1, @ML2), bioinformatics (@bio3, @bio2), neuroscience (@neuro), social network analysis (@SocialNetwork), and knowledge graph queries (@Hu2018-hd). More generally, the problem of discovering some true latent alignment between two networks can often be posed as variations on the above problem by adjusting the objective function for the setting.
The well-known graph isomorphism problem is a special case of GM problem when there exists a bijection between the nodes of two graphs which exactly preserves the edge structure. In terms of computational complexity, GM is equivalent to the NP-hard quadratic assignment problem, which is considered a very challenging problem where few theoretical guarantees exist, even in special cases (@QAP). For certain problems where the graphs are nearly isomorphic, polynomial-time algorithms do exist (@friendly, @Umeyama) but these methods frequently break down for more challenging instances.
This paper presents the detailed functionality of the iGraphMatch R package which serves as a practical tool for the use of prevalent graph matching methodologies.
These algorithms utilize either the spectral embedding of vertices (@Umeyama), or relaxations of the objective function (@PATH, @FAQ), or apply ideas from percolation theory (@Percolation, @ExpandWhenStuck).
The iGraphMatch package provides versatile options of working with graphs in the form of matrices, igraph objects or lists of either, and matching graphs under a generalized setting: weighted, directed, graphs of a different order, and multilayer graphs.
In addition, the iGraphMatch package incorporates prior information: seeds and similarities for all the implemented algorithms.
Seeds, or anchors, refer to partial knowledge of the alignment of two graphs.
In practice, seeds can be users with the same name and location across different social networks or pairs of genes with the same DNA sequences.
Some algorithms like the percolation algorithm (@Percolation, @ExpandWhenStuck) which matches two graphs by propagating matching information to neighboring pairs require seeds to kick off.
All algorithms improve substantially by incorporating seeds and can achieve accurate matching in polynomial time (@SGM).
Similarity scores are another commonly used prior which measures the similarity between pairs of nodes across the graphs.
In the bioinformatics area, BLAST similarity score is an example of similarity scores that plays an important role in aligning two PPI networks (@IsoRank).
Similarity scores are usually generated from nodal covariates that are observed in both networks (@BLAST, @SimilarityScore).
While under many scenarios the availability of exact partial matches, or hard seeding, is not realistic and expensive, the package also enables utilizing noisy prior information. Similarity scores incorporate uncertainty by assigning the pair of nodes with higher similarity scores a bigger chance to match. Seeds with uncertainty and even error can still be handled by self-correcting graph matching algorithms like the Frank-Wolfe algorithm initialized at the noisy partial matching, called soft seeding. @soft_seeding showed that the Frank-Wolfe algorithm with soft seeding scheme converges quickly to the true alignment under the correlated Erd\H{o}s-R\'enyi model with high probability. Thus, the original intractable problem is reduced to be solvable in polynomial time.
Although there exist some open source software and packages containing graph matching functionality, iGraphMatch package provides a centralized repository for common graph matching methodologies with flexibility, tools for developing graph matching problem methodology, as well as metrics for evaluating and tools for visualizing matching performance.
Among the alternative GM packages, the most relevant ones include
the igraph (@igraph) package which focuses on descriptive network analysis and graph visualization based on igraph objects and provides a single graph matching algorithm,
the GraphM (@PATH) package which implements several GM algorithms proposed between 1999 and 2009 in C, and the
Corbi (@Corbi) R package which is particularly designed for studies in bioinformatics and SpecMatch (@SpecMatch) which only involves implementations of spectral embedding based GM algorithms and written in C/C++.
None of these packages provide the breadth of tools, flexibility, and ease-of-use provided by the iGraphMatch package.
The rest of this paper is organized as follows. Section \@ref(sec:background) describes the theoretical representations of the implemented GM algorithms, correlated random graph models and evaluation metrics. Section \@ref(sec:usage) discusses the functionality and usage of R functions in the package, illustrated on the synthetic correlated graph pairs. Section \@ref(sec:example) presents more complex examples on real data with several functions involved in the analysis and section \@ref(sec:conclusion) gives guidelines on using different GM algorithms under different circumstances and concludes the paper.
Graph matching background {#sec:background}
In this section, we give background on graph matching and related problems followed by descriptions of the principal algorithms implemented in iGraphMatch.
For simplicity, we state all the algorithms in the context of matching undirected, unweighted graphs with the same cardinality.
All algorithms can also be directly applied to directed and weighted graphs.
In the second subsection, we discuss the techniques for matching graphs with a different number of vertices along with other extensions.
To conclude the section, we introduce the statistical models for correlated networks and discuss measures for the goodness of matching.
For the remainder of this paper we use the following notation.
Let $G_1=(V_1,E_1)$ and $G_2=(V_2,E_2)$ denote two graphs with $n$ vertices.
Let $A$ and $B$ be their corresponding binary symmetric adjacency matrices.
In the setting of seeded graph matching, suppose without loss of generality, the first $s$ pairs of nodes are seeds for simplicity.
In iGraphMatch, much more flexible seed specifications are possible, which will be illustrated in examples for usage of the package in section \@ref(sec:usage).
Accordingly, let $A$ and $B$ be partitioned as:
\begin{equation} \label{eq:seed_blocks}
A=
\begin{bmatrix}
A_{11} & A_{21}^T \
A_{21} & A_{22}
\end{bmatrix}\text{ and }
B=
\begin{bmatrix}
B_{11} & B_{21}^T \
B_{21} & B_{22}
\end{bmatrix}
\end{equation}
where $A_{11}, B_{11}\in{0, 1}^{s\times s}$ denote seed-to-seed adjacencies, $A_{21}, B_{21}\in{0, 1}^{(n-s)\times s}$ denote nonseed-to-seed adjacencies and $A_{22}, B_{22}\in{0, 1}^{(n-s)\times (n-s)}$ denote nonseed-to-nonseed adjacencies.
Let $S$ be an $n$-by-$n$ real-valued matrix of similarity scores.
Let $\Pi$ be the set of all permutation matrices and $\mathcal{D}$ be the set of all doubly stochastic matrices.
Assignment problems
Matching or assignment problems are core problems in combinatorial optimization and appear in numerous fields (@AP). As we illustrate in Eq. \eqref{eq:ob_func}, a general version of the graph matching problem is equivalent to the quadratic assignment problem (QAP). Similarly, QAP is related to the linear assignment problem (LAP) which also plays a role in GM. The LAP asks how to assign $n$ items (eg. workers or nodes in $G_1$) to $n$ other items (eg. tasks or nodes in $G_2$) with minimum cost. Let $C$ denote an ${n\times n}$ cost matrix, where $C_{ij}$ denotes the cost of matching $i$ to $j$, then the LAP is to find \begin{equation} \label{eq:LAP} \begin{aligned} & \argmin_{P\in \Pi} & \mathrm{trace}(C^TP) \end{aligned} \end{equation} LAP is solvable in $O(n^3)$ time and there are numerous exact and approximate methods for both general (@lap_solver, @Hungarian) and special cases, such as sparse cost matrices (@Volgenant1996-o).
The statement of QAP resembles LAP, except that the cost function is expressed as a quadratic function. Given two $n\text{-by-}n$ matrices $A$ and $B$ which can represent flows between facilities and the distance between locations respectively, or the adjacency matrices of two unaligned graphs, the objective function for QAP is: \begin{equation}\label{eq:qap} \argmin_{P\in \Pi} \mathrm{trace}(APBP^T). \end{equation} This problem is NP-hard (@QAP) leading to a core challenge for any graph matching approach.
As will be illustrated in the rest of the section, some matching algorithms reduce the graph matching problem to solving a LAP. For these algorithms, we include similarity scores $S$ by adding an additional term $\mathrm{trace}(S^TP)$ to the reduced objective function.
Graph matching algorithms {#subsec:gm}
In the iGraphMatch package, we implement three types of prevalent GM algorithms.
The first group uses relaxations of the objective function, including convex, concave, and indefinite relaxations.
The second group consists of algorithms that apply ideas from percolation theory, where matching information is spread from an initial set of matched nodes.
The last group is based on the spectral embedding of vertices.
Relaxation-based algorithms
These approaches relax the constraint that $P$ is a permutation matrix to require only that $P$ is doubly stochastic, optimizing over $\mathcal{D}$, the convex hull of $\Pi$. When $P$ is a permutation matrix \begin{equation} \label{eq:ob_func} \lVert A-PBP^T \rVert_F^2 = \lVert AP-PB \rVert_F^2 = \lVert A \rVert_F^2 + \lVert B \rVert_F^2 - 2\cdot \mathrm{trace}APBP^T. \end{equation} However, these equalities do not hold for all $P\in \mathcal{D}$, leading to different relaxations.
The second term of Eq.\eqref{eq:ob_func} is a convex function and optimizing it over $P\in\mathcal{D}$ gives the convex relaxation, where the gradient at $P$ to the convex relaxed objective function is $-4 APB + 2A^TAP + 2PBB^T$. The last equality in Eq. \eqref{eq:ob_func} shows that minimizing edge disagreements is equivalent to maximizing the number of edge agreements, $\mathrm{trace} APBP^T$, a QAP. Optimizing the indefinite function over $\mathcal{D}$ gives the indefinite relaxation with gradient $-2APB$ (@FAQ).
| Relaxation | Objective Function | Domain | GM Algorithm | Optimization Guarantee | Optimum Form | | ---------- | ---------------- | ------------- | ------------ | ---------------- | -------------- | | None | $\lVert A-PBP^T \rVert_F^2$ | $\Pi$ | NA | | $\Pi$ | | Indefinite | $\mathrm{tr}ADBD^T$ | $\mathcal{D}$ | FW | Local | $\mathcal{D}$ (often $\Pi$) | | Convex | $\lVert AD-DB \rVert_F^2$ | $\mathcal{D}$ | FW, PATH | Global | $\mathcal{D}$ | | Concave | $-\mathrm{tr}(\Delta D)-2\mathrm{tr}(L_1^T D L_2 D^T)$ | $\mathcal{D}$ | PATH | Local | $\Pi$ |
Table: (#tab:relaxation) Summary of relaxation methods for graph matching problem.
Generally, the convex relaxation leads to a solution that is not guaranteed to be near the solution to the original GM. However, @friendly introduced the class of "friendly" graphs based on the spectral properties of the adjacency matrices to characterize the applicability of the convex relaxation. Matching two friendly graphs by using the convex relaxation is guaranteed to find the exact solution to the GM problem. Unfortunately, this class is quite limiting and does not hold for most statistical models or real-world examples.
Another relaxation is the concave relaxation used in the PATH algorithm (@PATH). The concave relaxation uses the Laplacian matrix defined as $L=D-A$, where $D$ is the diagonal degree matrix with diagonal entries $D_{ii}=\sum_{i=1}^N A_{ij}$. Assume $L_i$ and $D_i$, $i=1,2$, are the Laplacian matrices and degree matrices for $G_1$ and $G_2$ respectively, then we can rewrite the objective function as \begin{equation} \label{eq:concave} \begin{split} \lVert A-PBP^T\rVert_F^2 & = \lVert AP-PB \rVert_F^2\ & = \lVert (D_1P-PD_2)-(L_1P-PL_2)\rVert_F^2\ & = -\mathrm{trace}(\Delta P)+\mathrm{trace}(L_1^2)+\mathrm{trace}(L_2^2)-2\mathrm{trace}(L_1^T P L_2 P^T), \end{split} \end{equation} where the matrix $\Delta_{ij}=(D_{2_{jj}}-D_{1_{ii}})^2$. Dropping the terms not dependent on $P$ in equation \ref{eq:concave}, we obtain the concave function $-\mathrm{trace}(\Delta P)-2\mathrm{trace}(L_1^T P L_2 P^T)$ on $\mathcal{D}$.
A summary of the different relaxations is provided in Table \ref{tab:relaxation}. Relaxing the discrete problem to a continuous problem breaks the equivalence to the original formulation of the edge disagreement and enables employing algorithms based on gradient descent.
Frank Wolfe methodology
@FAQ introduced an algorithm for the relaxed graph matching problem, with each iteration computable in polynomial time, that can find local optima for the relaxations above. The \Frank-Wolfe (FW) (@FW) methodology is an iterative gradient ascent approach composed of two steps. The first step finds an ascent direction that maximizes the gradient ascent. In this case the ascent direction is a permutation matrix which is a vertex of the polytope of doubly stochastic matrices. For the convex, indefinite, and concave relaxations, this corresponds to a LAP with the gradient as the cost function. The second step performs a line search along the ascent direction to optimize the relaxed objective function. As the objectives are all quadratic, this line search is simply requires optimizing a single-variable quadratic function along a line segment. After the iterative algorithm converges, the final step of the procedure is to project the doubly stochastic matrix back to the set of permutation matrices, which is also a LAP.
The various relaxed forms can all serve as the objective function $f(\cdot)$ in the FW Methodology, but in all cases a matrix $D^0\in \mathcal{D}$ must be chosen to initialize the procedure. For the convex relaxation, the FW methodology is guaranteed to converge to the global optimum regardless of the $D^0$. On the other hand, the FW algorithm for the indefinite relaxation is not guaranteed to find a global optimum so the initialization is critical.
In many instances, the optimal solution to the convex relaxation lies in the interior of $\mathcal{D}$. This can lead to inaccurate solutions after the last projection step. The local optima for the indefinite relaxation are often at extreme points of $\mathcal{D}$, meaning the final projection often does nothing.
The default initialization for the indefinite problem is at the barycenter matrix, $D^0 = \frac{1}{n}11^T$, but many other initialization procedures can be used. These include randomized initializations, initializations based on similarity matrices, and initializing the indefinite relaxation at the interior point solution of the convex relaxation (@relax_paper). When prior information regarding a partial correspondence is known to be noisy, rather than incorporating this information as seeds, one can incorporate it as "soft" seeds which are used to generate the initialization (@soft_seeding).
\begin{algorithm} \caption{Frank-Wolfe Methodology} \label{alg:FAQ} \SetKwInOut{Input}{Input} \SetKwInOut{Output}{Output} \Input{$A,B$, doubly stochastic matrix $D^0$, tolerance $\epsilon$} \Output{permutation matrix $P$} Set $i=0$;\par \While{$\Vert D^i-D^{i-1} \Vert^2_F\ge \epsilon$}{ $P^i=$ $\argmin_{P\in\Pi}\mathrm{trace}\nabla f(D^i)^TP$;\par $D^{i+1}=$ $\argmin_{D\in\mathcal{D}}f(D)$ over line segment from $D^i$ to $P^i$;\par $i=i+1$;\par } Project $D^{i}$ to the nearest $P$ by maximizing $\mathrm{trace}(P^TD)$;\par \end{algorithm}
When prior information is available in the form of seeds, the seeded graph matching problem (@SGM) works on the objective function \ref{eq:ob_func} with the permutation matrix $P^{n\times n}$ substituted by $I_s\oplus P^{(n-s)\times (n-s)}$, the direct sum of an $s\times s$ identity matrix and an $(n-s)\times (n-s)$ permutation matrix.
Employing the indefinite relaxed objective function incorporating seeds, we formulate the problem as finding \begin{align} \hat{P} % &= \argmin_{P\in \Pi} \mathrm{trace} A(I_s\oplus P)B(I_s\oplus P)^T \ &= \argmax_{P\in\mathcal{D}} 2\cdot\mathrm{trace}P^TA_{21}B_{21}^T+\mathrm{trace}A_{22}PB_{22}P^T \end{align} where the gradient to the objective function is \begin{equation} \label{eq:sgm_gradient} \nabla f(P)=2\cdot A_{21}B_{21}^T+2\cdot A_{22}PB_{22}. \end{equation}
In total, this uses the information between seeded nodes and nonseeded nodes and the nonseed-to-nonseed information. Applying seeded graph matching to the convex relaxation and concave relaxation closely resembles the case of indefinite relaxation.
PATH algorithm
@PATH introduced a convex-concave programming approach to approximately solve the graph matching problem. The concave relaxation has the same solution as the original graph matching problem. The PATH algorithm finds a local optimum to the concave relaxation by considering convex combinations of the convex relaxation $F_0(P)$ and the concave relaxation $F_1(P)$ denoted by $F_{\lambda} =(1 - \lambda) F_0 + \lambda F_1$. Starting from the solution to the convex relaxation ($\lambda=0$) the algorithm iteratively performs gradient ascent using the FW methodology at $F_\lambda$, increasing $\lambda$ after each iteration, until $\lambda = 1$.
Percolation-based algorithms
Under the FW methodology, all the nodes admit a correspondence but the (relaxed) matching correspondence evolves through iterations. On the other hand, percolation approaches start with a set of seeds, adding one new match at each iteration. The new matches are fixed and hence not updated in future iterations.
Each iteration expands the set of matched nodes by propagating the current matching information to neighbors. The guiding intuition is that more matched neighbors are an indicator of a more plausible match, an intuition analogous to the gradient ascent approaches above. We will present two algorithms in this category where the ExpandWhenStuck algorithm is an extension to the Percolation algorithm.
There are some distinctions about the inputs and outputs of percolation methods compared to the above relaxation methods.
Percolation Algorithm
@Percolation provides a simple and fast approach to solve the graph matching problem by starting with a handful of seeds and propagating to the rest of the graphs. At each iteration, the matching information up to the current iteration is encoded in a subpermutation matrix $P$ where $P_{ij}=1$ if $i$ is matched to $j$, and $0$ otherwise.
The Percolation algorithm searches for the most promising new match among the unmatched pairs through the mark matrix, $M=APB$, which is the gradient of the indefinite relaxation when extended to sub-doubly stochastic matrices. When similarity scores are available, they are added to the mark matrix to combine topological structure and similarity scores.
Adopting analogous partitions on the adjacency matrices as in equation \eqref{eq:seed_blocks}, we let $A_{21}, B_{21}$ denote sub-matrix corresponding to potential adjacencies between unmatched and matched nodes. Since all the candidates of matched pairs are permanently removed from consideration, we need only consider $M'=A_{21}B_{21}^T$, the sub-matrix of $M$ corresponding to the unmatched nodes in both graphs. As a result, the Percolation algorithm only uses matched-to-unmatched information to generate new matches.
Moreover, the mark matrix $M$ can also be interpreted as encoding the number of matched neighboring pairs for each pair of nodes $i\in V_1$, $j\in V_2$.
Suppose $u, u'\in V_1$, $v,v'\in V_2$, $[u,u']\in E_1$ and $[v,v']\in E_2$, then $(u',v')$ is a neighboring pair of $(u,v)$.
In each iteration, while there remain unmatched nodes with more than $r$ matched neighboring pairs, the percolation algorithm matches the pair of nodes with the highest score $M_{uv}$,
and adds one mark to all the neighboring pairs of $(u,v)$.
Note that the algorithm may stop before all nodes are matched, leading to the return of a partial match.
There is only one tuning parameter in the Percolation algorithm, the threshold $r$ which controls a tradeoff between quantity of matches and quality of matches. With a small threshold, the algorithm has a larger chance of matching wrong pairs. If $r$ is larger, then the algorithm might stop before matching many pairs (@ExpandWhenStuck).
The Percolation algorithm can be generalized to matching weighted graphs by making an adjustment to how we measure the matching information from the neighbors. Since we prefer to match edges with smaller weight differences and higher absolute weights, we propose to adopt the following update formula for the score associated with each pair of nodes $(i,j)$: $$M_{ij}=M_{ij} + \sum_{u\in N(i)}\sum_{v\in N(j)}1-\frac{|w_{iu}-w_{jv}|}{\max(|w_{iu}|, |w_{jv}|)}.$$ Thus, the score contributed by each neighboring pair of $(i,j)$ is a number in $[0,1]$.
\begin{algorithm}[H] \SetAlgoLined \SetKwInOut{Input}{Input} \SetKwInOut{Output}{Output} \Input{$A$, $B$, $s$ pairs of seeds, threshold $r$} \Output{(sub-)permutation matrix $P$} Initialize the sub-permutation matrix $P$ incorporating seeds;\par Calculate the mark matrix $M=APB$;\par Set the rows and columns of $M$ corresponding to seeds to minus infinity;\par \While{$max(M)\ge r$}{ $P_{ij}\leftarrow 1$, where $[i,j]$ is the index of $max(M)$;\par $M\leftarrow APB$;\par Set $i^{th}$ row and $j^{th}$ column of $M$ to minus infinity;\par } \caption{Percolation Algorithm} \label{alg:perco} \end{algorithm}
ExpandWhenStuck Algorithm
@ExpandWhenStuck extends the Percolation algorithm to a version that can operate with a smaller number of seeds. Without enough seeds, when there are no more unmatched pairs with a score higher or equal to the threshold $r$, the Percolation algorithm would stop even if there are still unmatched pairs. ExpandWhenStuck uses all pairs of nodes with at least one matched neighboring pair, $M_{ij}\geq 1$, as new seeds to restart the matching process by adding one mark to all of the new seeds' neighboring pairs, without updating the matched set. If the updated mark matrix consists of new pairs with marks greater or equal to $r$, then the percolation algorithm continues, leading to larger matched sets.
Spectral-based algorithm
Another class of graph matching algorithms uses the spectral properties of adjacency matrices.
IsoRank Algorithm
@IsoRank propose the IsoRank algorithm that uses neighboring topology and similarity scores and exploits spectral properties of the solution. The IsoRank algorithm is also based on the relaxation-based algorithms by encoding the topological structure of two graphs in $ADB$, which is again proportional to the gradient of the indefinite relaxation. However, the representations of each term of $ADB$ are slightly different. $A$ and $B$ are the column-wise normalized adjacency matrices and $D$ is not necessarily a doubly stochastic matrix yet $D_{ij}$ still indicates how promising it is to match $i\in V_1$ to $j\in V_2$.
Similar to the idea of Percolation algorithm, the intuition is that the impact of a pair of matched nodes is evenly distributed to all of their neighbors to propagate plausible matches. This is achieved by solving the eigenvalue problem \begin{align}\label{eq:Iso1} \mathrm{vec}(D)=(A\otimes B) \mathrm{vec}(D), \end{align} where $\mathrm{vec}(D)$ denotes the vectorization of matrix $D$, and the right hand side is equivalent to $ADB$.
To combine network-topological structure and similarity scores in the objective function, the normalized similarity score $E$ is added to the right hand side of Eq. \ref{eq:Iso1}, where $E=S/\|S\|_1$, and $|\cdot|_1$ denotes the L1 norm.
Note that when similarity score is not available as prior information, we can also construct a doubly stochastic similarity score matrix from seeds by taking $I_{s\times s}\oplus \frac{1}{n-s}11_{(n-s)\times (n-s)}$. To solve the eigenvalue problem \ref{eq:Iso1}, we resort to the power method. Finally, the global alignment is generated by a greedy algorithm or using the algorithms for solving the linear assignment problem (LAP).
\begin{algorithm} \caption{IsoRank Algorithm} \label{alg:IsoRank} \SetKwInOut{Input}{Input} \SetKwInOut{Output}{Output} \Input{$A,B$, similarity scores $S$} \Output{permutation matrix $P$} Column-wise normalize the adjacency matrices $A$ and $B$;\par Normalize similarity scores: $E=S/|S|_1$;\par Initialize $D^0=E$ and tolerance $\epsilon$;\par \While{$| D^i-D^{i-1}|\ge \epsilon$}{ Calculate $D^{i+1}=AD^iB+E$;\par Normalize $D^{i+1}=D^{i+1}/|D^{i+1}|$;\par $i=i+1$;\par } Extract node mapping from $\hat{D}$ using a greedy method or by solving a LAP; \end{algorithm}
Umeyama algorithm
@Umeyama is a spectral approach to find approximate solutions to the graph matching problem. Assuming eigendecompositions of adjacency matrices $A$ and $B$ as $A=U_A\Lambda_AU_A^T$ and $B=U_B\Lambda_BU_B^T$, let $|U_A|$ and $|U_B|$ be matrices which takes absolute values of each element of $U_A$ and $U_B$. Such modification to the eigenvector matrices guarantees the uniqueness of eigenvector selection. The global mapping is obtained by minimizing the differences between matched rows of $U_A$ and $U_B$: \begin{align} \hat{P}=\mathop{\argmin}{P\in\Pi}\lVert |U_A|-P|U_B|\rVert_F=\mathop{\argmax}{P\in\Pi}\mathrm{trace}(|U_B||U_A|^TP) \end{align}
The Umeyama algorithm can be generalized to matching directed graphs by eigendecomposing the Hermitian matrices $E_A$ and $E_B$ derived from the asymmetric adjacency matrices of the directed graphs. The Hermitian matrix for the adjacency matrix $A$ is defined as $E_A=A_S+iA_N$, where $A_S=(A+A^T)/2$ is a symmetric matrix, $A_N=(A-A^T)/2$ is a skew-symmetric matrix and $i$ is the imaginary unit. Similarly, we can define the Hermitian matrix for $B$. Assume the eigendecompositions of $E_A$ and $E_B$ as follows: \begin{align} E_A=W_A\Gamma_AW_A^, \quad E_B=W_B\Gamma_BW_B^ \end{align} and we aim at searching for: \begin{align} \hat{P}=\mathop{\argmax}_{P\in\Pi}\mathrm{trace}(|W_B||W_A|^TP) \end{align} Note that the Umeyama algorithm works on the condition that two graphs are isomorphic or nearly isomorphic.
\begin{algorithm} \caption{Umeyama Algorithm} \label{alg:Umeyama} \SetKwInOut{Input}{Input} \SetKwInOut{Output}{Output} \Input{$A,B$} \Output{permutation matrix $P$} Compute the eigendecompositions of $A$ and $B$: $A=U_A\Lambda_AU_A^T$, $B=U_B\Lambda_BU_B^T$ ;\par Solve the LAP: $P=\mathop{\argmax}_{P\in\Pi}\mathrm{trace}(|U_B||U_A|^TP)$;\par \end{algorithm}
Auxiliary graph matching tools
Centering technique
Instead of encoding the non-adjacencies by zeros in the adjacency matrices, the centering technique (@centering) assigns negative values to such edges. The first approach is encoding non-adjacent node-pairs as $-1$ with centered adjacency matrices $\tilde{A}=2A-\textbf{J}$ and $\tilde{B}=2B-\textbf{J}$, where $\textbf{J}$ is a matrix of all ones. An alternative approach relies on modeling assumptions where the pair of graphs are correlated but do not share a global structure. We match $\tilde{A} = A-\Lambda_A$ and $\tilde{B} = B-\Lambda_B$, where $\Lambda$ is an $n\text{-by-}n$ matrix with $ij$-th entry denoting an estimated marginal probability of an edge. In general, $\Lambda$ is unknown but there are methods in the literature to estimate $\Lambda$.
Matching centered graphs changes the rewards for matching edges, non-edges, and the penalties for mismatches. Adapting the centering technique for the problem at hand can be used to find specific types of correspondences. This can also be combined with constructing multilayer networks out of single layer networks to match according to multiple criteria (@Lin, @GRAMPA). The centering technique can be applied to any of the implemented graph matching algorithm. It is especially useful when padding graphs with differing numbers of vertices to distinguish isolated vertices from padded vertices.
Padding graphs of different orders
Until this section, we have been considering matching two graphs whose vertex sets are of the same cardinality. However, matching graphs with different orders are commonly seen in real-world problems.
Suppose $A\in{0,1}^{n\times n}$ and $B\in{0,1}^{n_c\times n_c}$ with $n_c<n$.
One can then pad the smaller graph with extra vertices to match the order of the larger graph, $\tilde{B}=B\oplus \textbf{0}_{n-n_c}$ and match $A$ and $\tilde{B}$.
Every implemented graph matching algorithm in the iGraphMatch package automatically handles input graphs with a different number of vertices using sparse padding with minimal memory impact.
Since the isolated vertices and the padded vertices share the same topological structure, it can be useful to center the original graphs first then pad the smaller graph in the same manner. This approach serves to differentiate between isolated vertices the padded ones. It's theoretically verified that in the correlated \ErdosRenyi graph model, the centered padding scheme is guaranteed to find the true correspondence between the nodes of $G_1$ and the induced subgraph of $G_2$ under mild conditions even if $|V_1|\ll|V_2|$, but the true alignment is not guaranteed without centering (@centering).
Exploiting sparse and low-rank structure
Many real-world graphs, especially large graphs, are often very sparse with $o(n^2)$ and often $\theta(n)$ edges.
This can increase the difficulty of the graph matching problem due to the fact that there are fewer potential edges to match, but sparse graphs also come with computational advantages.
We rely on igraph and Matrix for efficient storage of these structures as well as the efficient implementation of various matrix operations.
We also use the LAPMOD algorithm for sparse LAP problems (@Volgenant1996-o) (see below).
Similarly, a low-rank structure appears in some of the procedures including starting at the rank-1 barycenter matrix and the different centering schemes.
Since low-rank matrices are generally not sparse and visa-versa we implemented the splr S4 class, standing for sparse plus low-rank matrices.
This class inherits from the Matrix class and includes slots for an $n\times n$ sparse matrix x and $n\times d$ dense matrices a and b, to represent matrices of the form x + tcrossprod(a, b).
This class implements efficient methods for matrix multiplication and other operations that exploit the sparse and low-rank structure of the matrices.
Specifically, these methods often require only $O(\|x\|_0) + O(nd)$ storage as opposed to $O(n^2)$ required for densely stored matrices, and enjoy analogous computational advantages.
While users can also use these matrices explicitly, most use of them is automatic within functions such as init_start and center_graph and the matrices can largely be used interchangeably with other matrices.
LAP methods {#sssec:lap_methods}
Multiple graph matching methods include solving an LAP and so we have included multiple methods for solving LAPs into the package.
Specifically we implement the Jonker-Volgenant algorithm (@lap_solver) for dense cost matrices and the LAPMOD algorithm (@Volgenant1996-o) for sparse cost matrices.
Both algorithms are implemented in C to provide improved performance.
The LAPMOD approach is typically advantageous when the number of non-zero entries is less than 50%.
We also depend on the clue package for the solve_LSAP function which implements the Hungarian algorithm (@Papadimitriou1998-lz) for solving an LAP.
Each of these methods can be used independently of a specific graph matching method using the do_lap function.
Multi-layered graph matching
Frequently, networks edges may have categorical attributes and from these categories, we can construct multilayer graphs (@multilayer), where each layer in the networks contains edges from specific categories. For matching two multilayer graphs, the standard graph matching problem can be extended as $\sum_{l=1}^{m}\|A^{(l)} - PB^{(l)}P^T\|_F^2$ where $m$ denotes the number of categories and $A^{(l)}, B^{(l)}$ are the adjacency matrices for the $l$th layers in each graph. Note, we assume that the layers are aligned, so that layer $l$ corresponds to the same edge-types in both multi-layer networks.
For an igraph object, the function split_igraph can be used to convert a single object with categorical edge attributes into a list with each element only containing the edges with a specific attribute value.
The implemented algorithms can seamlessly match multi-layer graphs, which are encoded as a list of either igraph objects or matrix-like objects.
We also implemented a matrixlist S4 class that implements many standard matrix operations so that algorithms can be easily extended to work with multilayer networks.
Correlated random graph models {#graph-models}
The correlated \ErdosRenyi model (@FAQ) is essential in the theoretical study of graph matching algorithms. In a single \ErdosRenyi graph, each edge is present in the graph independently with probability $p$. The correlated \ErdosRenyi model provides a joint distribution for a pair of graphs, where each graph is marginally distributed as an \ErdosRenyi graph and corresponding edge-pairs are correlated.
\begin{defn}[Correlated \ErdosRenyi]
Suppose $p \in (0,1)$ and $\rho\in [\max{\frac{-p}{1-p}, \frac{p-1}{p}},1]$, a pair of adjacency matrices $(A, B)\sim \text{CorrER}(pJ, \rho J)$ if:
For each $1\le u \le v \le n, A_{uv}$ are independent with $A_{uv} \sim Bernoulli(p)$.\;
For each $1\le u \le v \le n, B_{uv}$ are independent with $B_{uv} \sim Bernoulli(p)$.\;
$A_{uv}$ and $B_{u'v'}$ are independent unless $u=u', v=v'$ in which case the Pearson correlation $\mathrm{corr}(A_{uv},B_{uv})=\rho$.
\end{defn}
To sample a pair of correlated \ErdosRenyi graphs with edge probability $p$, and Pearson correlation $\rho$, we first sample three independent \ErdosRenyi graphs $G_1$, $Z_0$ and $Z_1$ with edge probabilities $p$, $p(1-\rho)$ and $p+\rho(1-p)$ respectively. Let $G_2 = (Z_1 \cap G_1) \bigcup (Z_0\cap G_1^c)$.
@Percolation provide an alternative formulation for the correlated \ErdosRenyi model. First, one samples a single random \ErdosRenyi graph $G$ with edge probability $p'$. Conditioned on $G$, each edge in $G$ is present independently in $G_1,G_2$ with probability $s'$. These two parameterizations are related to each other by the relationship $s'=p+\rho(1-p)$ and $p'=p/(p+\rho(1-p))$. The original parameterization is slightly more general because it allows for the possibility of negative correlation.
In addition to homogeneous correlated \ErdosRenyi random graphs, we also implement heterogeneous generalizations of this model. The stochastic block model (@SBM) and the random dot product graphs (RDPG) model (@rdpg) can both be regarded as extensions of the \ErdosRenyi model. The stochastic block model is useful to represent the community structure of graphs by dividing the graph into $K$ groups. Each node is assigned to a group and the probability of edges is determined by the group memberships of the vertex pair. For the RDPG model, each vertex is assigned a latent position in $\mathbb{R}^d$ and edge probabilities are given by the inner product between the latent positions of the vertex pair.
For both of these models, we can consider correlated graph-pairs where marginally they arise from one of these models and again corresponding edge pairs are correlated.
Measures for goodness of matching {#sec:measure}
The ability to assess the quality of the match when ground truth is unavailable is critical for the usage of the matching approaches. There are various topological criteria that can be applied to measure the quality of matching results. At the graph level, the most frequently used structural measures include matching pairs (MP), edge correctness (EC), and the size of the largest common connected subgraph (LCCS) (@measure). MP counts the number of correctly matched pairs of nodes, thus can only be used when the true alignment is available. Global counts of common edges (CE) and common non-edges (CNE) can be defined as \begin{align} CE=\frac{1}{2}\sum_{i,j}1{A_{ij}=PBP^T_{ij} =1}\quad CNE=\frac{1}{2}\sum_{i,j}1{A_{ij}=PBP^T_{ij} = 0}, \end{align} along with error counts such extra edges (EE) and missing edges (ME), \begin{align} EE=\frac{1}{2}\sum_{i,j}1{A_{ij}=0=1 - PBP^T_{ij}}\quad ME=\frac{1}{2}\sum_{i,j}1{A_{ij}= 1 = 1 - PBP^T_{ij}}, \end{align}
EC measures the percentage of correctly aligned edges, that is the fraction CE$/|E_1|$. The LCCS denotes the largest subset of aligned vertices such that the corresponding induced subgraphs of each graph are connected. Matches with a larger LCCS are often preferable to those with many isolated components.
Another group of criteria measures the goodness of matching at the vertex level. Informally, we aim at testing the hypotheses $$H_0^{(v)}: \text{the vertex } v \text{ is not matched correctly by } P^, $$ $$H_a^{(v)}: \text{the vertex } v \text{ is matched correctly by } P^$$ for each vertex $v$.
The goal is to address if the permutation matrix found by graph matching algorithm is significantly different from the one sampled from uniformly distributed permutation matrices (@row_perm). Unfortunately, vertex-level matching criteria have only received limited attention in the literature, however, we include two test statistics to measure fit. The row difference statistic is the L$1$-norm of the difference between $A$ and $P^B{P^}^T$, namely \begin{equation} T_d(v,P^):=\Vert A{v\cdot}-(P^B{P^}^T){v\cdot}\Vert_1. \end{equation} Intuitively, a correctly matched vertex $v$ should induce a smaller $T_d(v,P^)$, which for unweighted graphs corresponds to the number of edge disagreements induced by matching $v$. Alternatively, the row correlation statistic is defined as \begin{equation} T_c(v,P^):=1 - corr(A{v\cdot},(P^B{P^}^T)_{v\cdot}). \end{equation} We expect the empirical correlation between the neighborhoods of $v$ in $A$ and $P^B{P^*}^T$ to be larger for a correctly matched vertex.
We employ permutation testing ideas to the raw statistics as a normalization across vertices. Let us take the row difference statistic for example. The guiding intuition is that if $v$ is correctly matched, the number of errors induced by $P^$ across the neighborhood of $v$ in $A$ and $B$ (i.e., $T_d(v, P^)$) should be significantly smaller than the number of errors induced by a randomly chosen permutation $P$ (i.e., $T_d(v, P)$).
With this in mind, let $\mathbb{E}_P$ and $\mathrm{Var}_P$ denote the conditional expectation and variance of the raw statistic with $P$ uniformly sampled over all permutation matrices. The normalization is then given by \begin{equation} T_p(v,P^):=\frac{T(v,P^)-\mathbb{E}_PT(v,P)}{\sqrt{Var_PT(v,P)}} \end{equation} where $T(v, P)$ can be either of the two test statistics we introduced earlier.
In addition to measuring match quality, these vertex-wise statistics can also serve as a tool to find which vertices have no valid match in another network, i.e. the vertex entity is present in one network but not the other.
R functions and usage {#sec:usage}
The R package iGraphMatch offers versatile options for graph matching and subsequent analysis.
Here we describe the usage of the package including sampling random correlated graph pairs, graph matching, and evaluating matching results.
Sampling correlated random graph pairs
We first illustrate the usage of functions for sampling correlated random graph pairs. The usage of graph matching will be demonstrated on the graph-pairs sampled using these methods.
Functions of the form sample_correlated_*_pair for sampling random graph pairs have the common syntax:
sample_correlated_*_pair(n, ***model parameters***, permutation = 1:n, directed = FALSE, loops = FALSE)
The argument n specifies the number of nodes in each graph, and the default options are to sample a pair of undirected graphs without self-loop whose true alignment is the identity.
The permutation argument can be used to permute the vertex labels of the second graph.
The model parameters arguments vary according to different random graph models and typically consist of parameters for marginal graph distributions and for correlations between the corresponding edges.
The functions each return a named list of two igraph objects.
For the homogeneous correlated \ErdosRenyi graph model, the model parameters are p, the global edge probability, and corr, the Pearson correlation between aligned vertex-pairs.
For example, to sample a pair of graphs with 5 nodes from $\mathrm{CorrER}(0.5J, 0.7J)$ we run
knitr::opts_chunk$set(warning=FALSE, message=FALSE) options(scipen = 1, digits = 2)
library(iGraphMatch) set.seed(1) cgnp_pair <- sample_correlated_gnp_pair(n = 5, corr = 0.7, p = 0.5) (cgnp_g1 <- cgnp_pair$graph1) cgnp_g1[] cgnp_g2 <- cgnp_pair$graph2
Since we didn't obscure the vertex correspondence by assigning a value to the permutation argument, the underlying true alignment is the identity.
For the more general heterogeneous correlated \ErdosRenyi graph model, one needs to specify an edge probability matrix and a Pearson correlation matrix. To sample a pair of graphs from the heterogeneous correlated \ErdosRenyi model again with 5 nodes in each graph, and with random edge probabilities and Pearson correlations:
set.seed(123) p <- matrix(runif(5^2, .5, .8),5) c <- matrix(runif(5^2, .5, .8),5) ieg_pair <- sample_correlated_ieg_pair(n = 5, p_mat = p, c_mat = c)
Since the default is undirected graphs without self-loops, the entries of p and c along and below the diagonal are effectively ignored.
The stochastic block model requires block-to-block edge probabilities stored in the pref.matrix argument and the block.sizes argument indicates the size of each block, along with the Pearson correlation parameter corr.
Next, we sample a pair of graphs from the stochastic block model with two blocks of size 2 nodes and 3 nodes respectively, within-group edge probabilities of .7 and .5, across-group edge probability of .001, and Pearson correlation equal to .5.
pm <- cbind(c(.7, .001), c(.001, .5)) sbm_pair <- sample_correlated_sbm_pair(n = 5, pref.matrix = pm, block.sizes = c(2,3), corr = 0.5)
The iGraphMatch package also provides functions for sampling a pair of correlated random graphs with junk vertices, i.e vertices that don't have true correspondence in the other graph by specifying the number of overlapping vertices in the argument ncore or overlapping block sizes in the argument core.block.sizes.
The iGraphMatch package offers auxiliary tools for centering graphs to penalize the incorrect matches as well, which is implemented in the center_graph function with syntax:
center_graph(A, scheme = c(-1, 1), use_splr = TRUE)
with the first input being either a matrix-like or igraph object.
The scheme argument specifies the method for centering graphs.
Options include a pair of scalars where the entries of the adjacency matrix are linearly rescaled so that their minimum is min(scheme) and their maximum is max(scheme).
Note, scheme = "center" is the same as scheme = c(-1, 1).
Another option is to pass in a single integer, where the returned value is the adjacency matrix minus its best rank-scheme approximation.
The last argument use_splr is a boolean indicating whether to return a splrMatrix object.
We use use_splr = FALSE here to better display the matrices but use_splr = TRUE will often result in improved performance, especially for large sparse networks.
Here, we center the sampled graph cgnp_g1 using different schemes:
center_graph(cgnp_g1, scheme = "center", use_splr = FALSE) center_graph(cgnp_g1, scheme = 2, use_splr = FALSE)
Users can then use the centered graphs as inputs to the implemented graph matching algorithms, which serve to alter rewards and penalties for common edges, common non-edges, missing edges, and extra edges.
Graph matching methods
The graph matching methods share the same basic syntax:
gm(A, B, seeds = NULL, similarity = NULL, method = "indefinite", ***algorithm parameters***)
| Parameter | Type | Description | Functions |
| ------- | --------- | ---------------- | --------- |
| start | Matrix or character | Initialization of the start matrix for iterations. | FW, convex |
| lap_method | Character | Method for solving linear assignment problem. | FW, convex, PATH, IsoRank |
| max_iter | Number | Maximum number of iterations. | FW, convex, PATH, IsoRank |
| tol | Number | Tolerance of edge disagreements. | FW, convex, PATH |
| r | Number | Threshold of neighboring pair scores. | percolation |
| ExpandWhenStuck | Boolean | TRUE if performs ExpandWhenStuck algorithm. | percolation |
Table: (#tab:gm-alg) Overview of arguments for different graph matching functions.
The first two arguments for graph matching algorithms represent two networks which can be matrices, igraph objects, or two lists of either form in the case of multi-layer matching.
The seeds argument contains prior information on the known partial correspondence of two graphs.
It can be a vector of logicals or indices if the seed pairs have the same indices in both graphs.
In general, the seeds argument takes a matrix or a data frame as input with two columns indicating the indices of seeds in the two graphs respectively.
The similarity parameter is for a matrix of similarity scores between the two vertex sets, with larger scores indicating higher similarity.
Notably, one should be careful with the different scales of the graph topological structure and the vertex similarity information in order to properly address the relative importance of each part of the information.
The method argument specifies a graph matching algorithm to use, and one can choose from "indefinite" (default), "convex", "PATH", "percolation", "IsoRank", "Umeyama", or a self-defined graph matching function which enables users to test out their own algorithms while remaining compatible with the package.
If method is a function, it should take at least two networks, seeds and similarity scores as arguments.
Users can also include additional arguments if applicable.
The self-defined graph matching function should return an object of the "graphMatch" class with matching correspondence, sizes of two input graphs, and other matching details.
As an illustrative example, graph_match_rand defines a new graph matching function which matches by randomly permuting the vertex label of the second graph using a random seed rand_seed.
We then apply this self-defined GM method to matching the correlated \ErdosRenyi graphs sampled earlier with a specified random seed:
graph_match_rand <- function(A, B, seeds = NULL, similarity = NULL, rand_seed){ totv1 <- nrow(A[[1]]) totv2 <- nrow(B[[1]]) nv <- max(totv1, totv2) set.seed(rand_seed) corr <- data.frame(corr_A = 1:nv, corr_B = c(1:nv)[sample(nv)]) graphMatch( corr = corr, nnodes = c(totv1, totv2), detail = list( rand_seed = rand_seed ) ) } match_rand <- gm(cgnp_g1, cgnp_g2, method = graph_match_rand, rand_seed = 123)
Other arguments vary for different graph matching algorithms with an overview given in Table \ref{tab:gm-alg}.
The start argument for the FW methodology with "indefinite" and "convex" relaxations takes any $nns\text{-by-}nns$ matrix or an initialization method including "bari", "rds" or "convex".
These represent initializing the iterations at a specific matrix, the barycenter, a random doubly stochastic matrix, or the doubly stochastic solution from "convex" method on the same graphs, respectively.
Moreover, sometimes we have access to side information on partial correspondence with uncertainty.
If we still treat such prior information as hard seeds and pass them through the seeds argument for "indefinite" and "convex" methods, incorrect information can yield unsatisfactory matching results.
Instead, we provide the option of soft seeding by incorporating the noisy partial correspondence into the initialization of the start matrix.
The core function used for initializing the start matrix with versatile options is the init_start function.
Suppose the first two pairs of nodes are hard seeds and another pair of incorrect seed $(3,4)$ is soft seeds:
hard_seeds <- 1:5 <= 2 soft_seeds <- data.frame(seed_A = 3, seed_B = 4)
We generate a start matrix incorporating soft seeds initialized at the barycenter:
as.matrix(start_bari <- init_start(start = "bari", nns = 3, ns = 2, soft_seeds = soft_seeds))
An alternative is to generate a start matrix that is a random doubly stochastic matrix incorporating soft seeds as follow
set.seed(1) as.matrix(start_rds <- init_start(start = "rds", nns = 3, ns = 2, soft_seeds = soft_seeds))
Then we can initialize the Frank-Wolfe iterations at any of the start matrix by specifying the start parameter.
When there are no soft seeds, we no longer need to initialize the start matrix by using init_start first. Instead we can directly assign an initialization method to the start argument in the gm function:
match_rds <- gm(cgnp_g1, cgnp_g2, seeds = hard_seeds, method = "indefinite", start = "rds")
Below use solution from the convex relaxation as the initialization for the indefinite relaxation.
set.seed(123) match_convex <- gm(cgnp_g1, cgnp_g2, seeds = hard_seeds, method = "indefinite", start = "convex")
Now let's match the sampled pair of graphs \code{sbm_pair} from the stochastic block model by using Percolation algorithm. Apart from the common arguments for all the graph matching algorithms, Percolation has another argument \code{r} representing the minimum number of matched neighbors required for matching a new qualified vertex pair. Here we adopt the default value which is 2. Also, at least one of similarity scores and seeds is required for Percolation algorithm to kick off. Let's utilize the same set of hard seeds and assume there is no available prior information on similarity scores.
sbm_g1 <- sbm_pair$graph1 sbm_g2 <- sbm_pair$graph2 match_perco <- gm(sbm_g1, sbm_g1, seeds = hard_seeds, method = "percolation", r = 2) match_perco
Without enough prior information on partial correspondence, Percolation couldn't find any qualifying matches. Suppose in addition to the current pair of sampled graphs, the above sampled correlated homogeneous and heterogeneous \ErdosRenyi graphs are different layers of connectivity for the same set of vertices. We can then match the nonseed vertices based on the topological information in all of these three graph layers. To be consistent, let's still use the Percolation algorithm with threshold \code{r} equal to 2 and the same set of seeds.
matrix_lA <- list(sbm_g1, ieg_pair$graph1, cgnp_g1) matrix_lB <- list(sbm_g2, ieg_pair$graph2, cgnp_g2) match_perco_list <- gm(A = matrix_lA, B = matrix_lB, seeds = hard_seeds, method = "percolation", r = 2) match_perco_list
With the same amount of available prior information, we are now able to match all the nodes correctly.
Finally, we will give an example of matching multi-layers of graphs using IsoRank algorithm. Unlike the other algorithm, similarity scores are required for IsoRank algorithm. Without further information, we adopt the barycenter as the similarity matrix here.
set.seed(1) sim <- as.matrix(init_start(start = "bari", nns = 5, soft_seeds = hard_seeds)) match_IsoRank <- gm(A = matrix_lA, B = matrix_lB, seeds = hard_seeds, similarity = sim, method = "IsoRank", lap_method = "LAP")
Graph matching functions return an object of class "graphMatch" which contains the details of the matching results, including a list of the matching correspondence, a call to the graph matching function and dimensions of the original two graphs.
match_convex@corr match_convex@call match_convex@nnodes
Additionally, "graphMatch" also returns a list of matching details corresponding to the specified method.
Table \ref{tab:gm-value} provides an overview of returned values for different graph matching methods.
With the seeds information, one can obtain a node mapping for non-seeds accordingly
match_convex[!match_convex$seeds]
| Parameter | Description | Functions |
| :------- | :----------------: | ---------: |
| seeds | A vector of logicals indicating if the corresponding vertex is a seed | All the functions |
| soft | The functional similarity score matrix with which one can extract more than one matching candidates | FW, convex, PATH, IsoRank, Umeyama |
| lap_method | Choice for solving the LAP. | FW, convex, Umeyama, IsoRank |
| iter | Number of iterations until convergence or reaches the max_iter. | FW, convex, PATH |
| max_iter | Maximum number of replacing matches. | FW, convex |
| match_order | The order of vertices getting matched. | percolation, IsoRank |
Table: (#tab:gm-value) Overview of return values for different graph matching functions.
The "graphMatch" class object can also be flexibly used as a matrix. In addition to the returned list of matching correspondence, one can obtain the corresponding permutation matrix in the sparse form.
match_convex[]
Notably, multiplicity is applicable to the "graphMatch" object directly without converting to the permutation matrix. This enables obtaining the permuted second graph, that is $PBP^T$ simply by
match_convex %*% cgnp_g2
Evaluation of goodness of matching
Along with the graph matching methodology, iGraphMatch has many capabilities for evaluating and visualizing the matching performance.
After matching two graphs, the function summary can be used to get a summary of the overall matching result in terms of commonly used measures including the number of matches, the number of correct matches, common edges, missing edges, extra edges and the objective function value.
The edge matching information is stored in a data frame named edge_match_info.
Note that summary outputs the number of correct matches only when the true correspondence is known by specifying the true_label argument with a vector indicating the true correspondence in the second graph.
Applying the summary function on the matching result match_convex with true_label = 1:5, indicating the true correspondence is the identity that provides these summaries.
summary(match_convex, cgnp_g1, cgnp_g2, true_label = 1:5)
Applying the summary function to a multi-layer graph matching result returns edge statistics for each layer.
summary(match_IsoRank, matrix_lA, matrix_lB)
In realistic scenarios, the true correspondence is not available.
As introduced in section \ref{sec:background}, the user can use vertex level statistics to evaluate match performance.
The best_matches function evaluates a vertex-level metric and returns a sorted data.frame of the vertex-matches with the metrics.
The arguments are the two networks, a specific measure to use, the number of top-ranked vertex-matches to output, and the matching correspondence in the second graph if applicable.
As an example here, we apply best_matches to rank the matches from above with the true underlying alignment
best_matches(cgnp_g1, cgnp_g2, match = match_convex, measure = "row_perm_stat", num = 3, true_label = 1:igraph::vcount(cgnp_g1))
Note, best_matches uses seed information from the match parameter and only outputs non-seed matches.
Without the true correspondence, true_label would take the default value and the output data frame only contains the first three columns.
To visualize the matches of smaller graphs, the function plot displays edge discrepancies of the two matched graphs by an adjacency matrix or a ball-and-stick plot, depending on the input format of two graphs.
plot(cgnp_g1, cgnp_g2, match_convex) plot(cgnp_g1[], cgnp_g2[], match_convex)
The plots for visualizing matching performance of match_convex are shown in Figure \@ref(fig:visualization).
Grey edges and pixels indicate common edges, red ones indicate edges only in the second graph.
If they were present, blue pixels and edges represent missing edges that only exist in the first graph.
The corresponding linetypes are solid, short dash, and long dash.
Examples {#sec:example}
In this section, we demonstrate graph matching analysis using iGraphMatch via examples on real datasets, including communication networks, neuronal networks, and transportation networks.
Table \@ref(tab:dataset-overview) presents brief overviews of the first two datasets.
Note that the number of edges doesn't consider weights for weighted graphs, and for directed graphs, an edge from node $i$ to node $j$ and another edge from $j$ to $i$ will be counted as two edges.
Tables \@ref(tab:edge-summary) and \@ref(tab:edge-summary-trans) summarize the edge correspondence between two graphs under the true alignment including the number of common edges, missing edges, and extra edges in two graphs.
In the first Enron email network example, we demonstrate the usage of Frank-Wolfe methodology and how to improve matching performance by using the centering technique and incorporating adaptive seeds. In the second example using C. Elegans synapses networks, we illustrate how to use soft matching for a challenging graph matching task using Frank-Wolfe methodology, PATH algorithm and IsoRank algorithm. Finally, we include an example of matching two multi-layer graphs with similarity scores on the Britain transportation networks.
summarize_graph <- function(gp, name) { nn <- nrow(gp[[1]][]) if (nn != nrow(gp[[2]][])) { nn <- paste(nn, "/", nrow(gp[[2]][])) } ne <- paste(igraph::ecount(gp[[1]]), "/", igraph::ecount(gp[[2]])) data.frame( `Dataset` = name, `Nodes` = as.character(nn), `Edges` = ne, `Correlation` = cor(as.numeric(gp[[1]][]), as.numeric(gp[[2]][])), `Weighted` = ifelse(igraph::is_weighted(gp[[1]]), "Yes", "No") , `Directed` = ifelse(igraph::is_directed(gp[[1]]), "Yes", "No") , `Loop` = paste( ifelse(any(igraph::is.loop(gp[[1]])), "Yes", "No"), "/", ifelse(any(igraph::is.loop(gp[[2]])), "Yes", "No") ) ) } knitr::kable( rbind( summarize_graph(Enron, "Enron"), summarize_graph(C.Elegans, "C. Elegans") ), col.names = c("Dataset", "\\# Nodes", "\\# Edges", "Correlation", "Weighted", "Directed", "Loop"), booktabs = TRUE, caption = "Overview of the Enron and C. Elegans graphs.", escape = FALSE )
summarize_edge <- function(gp1, gp2, name) { nn <- nrow(gp1[]) corr <- graphMatch(data.frame(corr_A = 1:nn, corr_B = 1:nn), nn) edge_info <- summary(corr, gp1, gp2)$edge_match_info data.frame( `Dataset` = name, `Common Edges` = edge_info$common_edges, `Missing Edges` = edge_info$missing_edges, `Extra Edges` = edge_info$extra_edges ) } knitr::kable( rbind( summarize_edge(Enron[[1]], Enron[[2]], "Enron"), summarize_edge(C.Elegans[[1]], C.Elegans[[2]], "C. Elegans") ), col.names = c("Dataset", "Common edges", "Missing edges", "Extra edges"), booktabs = TRUE, caption = "Edge summary under the true alignments of the Enron and C. Elegans graphs. The columns indicate the number of common edges, missing edges in $G_1$, and extra edges in $G_2$. For weighted graphs, we define a pair of corresponding edges as a common edge as long as they both have positive weights. ", escape = FALSE )
Example: Enron Email Network Data {#sec:Enron}
The Enron email network data was originally made public by the Federal Energy Commission during the investigation into the Enron Corporation (@Enron). Each node of Enron network represents an email address and if there is at least one email sent from one address to another address, a directed edge exists between the corresponding nodes.
The iGraphMatch package includes the Enron email network data in the form of a pair of igraph objects derived from the original data where each graph represents one week of emails between 184 email addresses.
The two networks are unweighted and directed with edge densities around
r round((igraph::edge_density(Enron[[1]]) + igraph::edge_density(Enron[[2]]))/2, 3) in each graph and the empirical correlation between two graphs is r round(cor(as.numeric(Enron[[1]][]), as.numeric(Enron[[2]][])), 2).
First, let's load packages required for the following analysis:
library(igraph) library(iGraphMatch) library(purrr) library(dplyr)
Visualization of Enron networks
We visualize the aligned Enron networks using the \code{plot} function with vertices sorted by a community detection algorithm (@community_detection) and degree. For detailed interpretations to figure \@ref(fig:Enron-graph), please refer to figure \@ref(fig:visualization).
g <- igraph::as.undirected(Enron[[1]]) com <- igraph::membership(igraph::cluster_fast_greedy(g)) deg <- rowSums(as.matrix(g[])) ord <- order(max(deg)*com+deg) plot(Enron[[1]][][ord,ord], Enron[[2]][][ord,ord])
Note that r sum(degree(Enron[[1]]) == 0) and r sum(degree(Enron[[2]]) == 0) out of the total r vcount(Enron[[1]]) nodes are isolated from the other nodes in two graphs respectively, indicating the corresponding employees haven't sent or received emails from other employees.
This adds difficulty to matching since it's impossible to distinguish the isolated nodes based on topological structure alone.
We first keep only the largest connected component of each graph.
vid1 <- which(largest_cc(Enron[[1]])$keep) vid2 <- which(largest_cc(Enron[[2]])$keep) vinsct <- intersect(vid1, vid2) v1 <- setdiff(vid1, vid2) v2 <- setdiff(vid2, vid1) A <- Enron[[1]][][c(vinsct, v1), c(vinsct, v1)] B <- Enron[[2]][][c(vinsct, v2), c(vinsct, v2)]
The sizes of largest connect components of two graphs are r dim(A)[1] and r dim(B)[1], which are different.
We reorder two graphs so that the first r length(vinsct) nodes are aligned and common to both graphs.
Matching largest connected components using FW Algorithm
Let's assume the Enron email communication network from the second week is anonymous, and we aim at finding an alignment between the email addresses from the first network and the second one to de-anonymize the latter. Additionally, we want to find the email addresses that are active in both months.
Suppose no prior information on partial alignment is available in this example. We match the two largest connected components using the FW algorithm with indefinite relaxation since seeds and similarity scores are not mandatory for this method.
Without any prior information, \code{seeds} and \code{similarity} arguments take default values which are \code{NULL}. For the \code{start} argument, we assign equal likelihood to all the possible matches by initializing at the barycenter. Since two graphs are of different sizes, the \code{gm} function automatically pads the smaller graph with extra 0's.
set.seed(1) match_FW <- gm(A = A, B = B, start = "bari", max_iter = 200) head(match_FW)
Then, we check the summary of matching performance in terms of matched nodes, matched edges and the graph matching objective function.
summary(match_FW, A, B)
emi <- summary(match_FW, A, B)$edge_match_info
In this example, we can evaluate the matching result based on statistics on matched edges.
Without any seeds or similarity scores, around r paste(round(100*emi$common_edges / (emi$common_edges + emi$missing_edges), 0), "%", sep="") of edges are correctly matched.
Centering the larger graph
We can try to improve performance by centering B by assigning -1 to non-edges, so that we penalize edges that are missing in B but present in A.
A_center <- center_graph(A = A, scheme = "naive", use_splr = TRUE) B_center <- center_graph(A = B, scheme = "center", use_splr = TRUE) set.seed(1) match_FW_center <- gm(A = A_center, B = B_center, start = "bari", max_iter = 200) summary(match_FW_center, A, B)
From the summary tables, we would prefer matching Enron networks with the application of the centering scheme, since we get more matched common edges, as well as fewer missing edges and extra edges.
Matching with adaptive seeds
Supposing we have no access to ground truth, we use the \code{best_matches} function to measure and rank the vertex-wise matching performance. Below shows the 6 matches that minimize the row permutation statistic.
bm <- best_matches(A = A, B = B, match = match_FW_center, measure = "row_perm_stat") head(bm)
Since seeded graph matching enhances the graph matching performance substantially (@SGM), it may be useful to use some of these best matches as seeds to improve matching results. Here, we use adaptive seeds, taking the $ns$ best matches and using them as seeds in a second run of the matching algorithm. The table below displays edge statistics and objective function values for different number of adaptive seeds used. The second column in the table shows the matching precision of the adaptive seeds based on ground truth. Incorporating adaptive seeds and repeating the FW matching procedure on centered graphs further improve the matching results, compared with the case without any adaptive seeds when $ns=0$. The first 40 pairs of matched nodes ranked by \code{best_matches} function are all correctly matched, and this is also when matching is improved the most.
match_w_hard_seeds <- function(ns){ seeds_bm <- head(bm, ns) precision <- mean(seeds_bm$A_best == seeds_bm$B_best) match_FW_center_seeds <- gm(A = A_center, B = B_center, seeds = seeds_bm, similarity = NULL, start = "bari", max_iter = 100) edge_info <- summary(match_FW_center_seeds, A, B)$edge_match_info cbind(ns, precision, edge_info) } set.seed(12345) map_dfr(seq(from = 0, to = 80, by = 20), match_w_hard_seeds)
set.seed(12345) map_dfr(seq(from = 0, to = 80, by = 20), match_w_hard_seeds) %>% knitr::kable(col.names = c("ns", "precision","common", "missing", "extra", "fnorm"), booktabs = TRUE, digits = 2)
As the number of adaptive seeds increases, the precision of adaptive seeds decreases. Note that if they are treated as hard seeds, incorrect matches will remain in the matched set and might cause a cascade of errors. An alternative way is to treat the top-ranked matches as soft seeds embedded in the start matrix to handle the uncertainty. In this way, adaptive seeds not only provide prior information but also evolve over iterations. The table below shows that the soft seeding approach always outperforms or performs as good as the hard seeding approach regardless of the number of adaptive seeds being used.
match_w_soft_seeds <- function(ns){ seeds_bm <- head(bm, ns) precision <- mean(seeds_bm$A_best == seeds_bm$B_best) start_soft <- init_start(start = "bari", nns = max(dim(A)[1], dim(B)[1]), soft_seeds = seeds_bm) match_FW_center_soft_seeds <- gm(A = A_center, B = B_center, start = start_soft, max_iter = 100) edge_info <- summary(match_FW_center_soft_seeds, A, B)$edge_match_info cbind(ns, precision, edge_info) } set.seed(12345) map_dfr(seq(from = 0, to = 80, by = 20), match_w_soft_seeds)
set.seed(12345) map_dfr(seq(from = 0, to = 80, by = 20), match_w_soft_seeds) %>% knitr::kable(col.names = c("ns", "precision","common", "missing", "extra", "fnorm"), booktabs = TRUE, digits = 2)
Core vertices detection
The function \code{best_matches} can also be used to detect core vertices. Suppose the ground truth is known and that the first 145 vertices are core vertices. The mean precision of detecting core vertices and junk vertices using \code{best_matches} function is displayed in figure \ref{fig:core}. A lower rank is a stronger indicator of a core vertex and a higher rank is a stronger indicator of a junk vertex. Let $r^C_i, 1\le i\le n_c$ and $r^J_j, 1\le j\le n_j$ denote the ranks associated with each core vertex and each junk vertex. The figure shows the precision of identifying core vertices at each low rank $r$, i.e. $\frac{1}{r}\sum_{i = 1}^{n_c}1_{r^C_i\le r}$, and the precision of identifying junk vertices at each high rank $r$, i.e. $\frac{1}{r}\sum_{j = 1}^{n_j}1_{r^J_j\ge n_c+n_j-r}$, which are separated by the vertical lines.
nc <- length(vinsct) nj <- max(length(v1), length(v2)) core_precision <- 1:nc %>% map_dbl(~mean(bm$A_best[1:.x]<=nc)) junk_precision <- 1:nj %>% map_dbl(~mean(bm$A_best[(nc+.x):(nc+nj)]>nc))
Core detection performance is substantially better than chance, as represented by the dotted horizontal lines.
The top r max(which(core_precision == 1)) are all core vertices indicating good overall performance for core identification.
For junk identification, the junk vertices are ranked r nc + nj - which(bm$A_best > nc) + 1 according to which have the lowest score, indicating that some junk vertices are difficult to identify.
```rMean precision for identifying core and junk vertices for the Enron networks by using the row permutation test. The vertical lines separate the performance of identifying core vertices with low ranks from junk vertices with high ranks. The horizontal lines indicate the performance of a random classifier."} library(ggplot2)
map_df_core <- data.frame(rank = 1:nc, which = "core", mean_precision = core_precision) map_df_junk <- data.frame(rank = (nc+1):(nc+nj), which = "junk", mean_precision = junk_precision) map_df <- rbind(map_df_core, map_df_junk) map_df <- map_df %>% mutate(random = ifelse(which == "junk", nj/(nc+nj), nc/(nc+nj))) %>% mutate(start=ifelse(which=="core",0,nc-5), end=ifelse(which=="core",nc+5,nc+nj))
map_df %>% ggplot(aes(x=rank,y=mean_precision)) + geom_line(alpha=.5,size=.3) + geom_segment(aes(y=random,yend=random, x=start,xend=end), data=map_df,color="black",alpha=1,linetype=2,size=.1)+ geom_vline(aes(xintercept=nc),size=.1,linetype=2)+ ylab("mean precision")+ theme_minimal()+ theme(text=element_text(family="Palatino"), axis.line=element_line(size=.3,color="black"), panel.grid.major=element_line(color="grey95",size=0), panel.grid.minor=element_line(color="grey98",size=0))
## Example: C. Elegans Network Data {#sec:CE} The _C. Elegans_ networks consist of the chemical synapses network and the electrical synapses network of the roundworm, where each of `r vcount(C.Elegans[[1]])` nodes represents a neuron and each edge represents the intensity of synapse connections between two neurons (@neuro). Matching the chemical synapses network to the electrical synapses network is essential for understanding how the brain functions. These networks are quite sparse with edge densities of `r round(igraph::edge_density(C.Elegans[[1]]), 3)` and `r round(igraph::edge_density(C.Elegans[[2]]), 3)` in each graph and the empirical correlation between two graphs is `r round(cor(as.numeric(C.Elegans[[1]][]), as.numeric(C.Elegans[[2]][])), 2)`. ### A challenging task For simplicity, we made the networks unweighted and undirected for the experiments, and we assume the ground truth is known to be the identity. ```r C1 <- C.Elegans[[1]][] > 0 C2 <- C.Elegans[[2]][] > 0 plot(C1[], C2[]) match <- gm(C1, C2, start = Matrix::Diagonal(nrow(C1))) plot(C1[], C2[], match)
nv <- nrow(C1) id_match <- graphMatch(data.frame(corr_A = 1:nv, corr_B = 1:nv), nv) i_sum <- summary(id_match, C.Elegans[[1]], C.Elegans[[2]]) m_sum <- summary(match, C.Elegans[[1]], C.Elegans[[2]], id_match) i_emi <- i_sum$edge_match_info m_emi <- m_sum$edge_match_info
Matching the C. Elegans networks is a challenging task.
Figures \@ref(fig:C-Elegans-edge) depict the edge discrepancies of two networks under the true alignment and the matching correspondence using FW algorithm initialized at the true alignment.
The alignment found using FW is not the identity with r m_sum$n_true_match out of r nv nodes correctly matched
and improves upon the identity in terms of the number of edge discrepancies.
For the true alignment, there are
r i_emi$common_edges
edge errors and
r round(i_emi$extra_edges + i_emi$missing_edges)
common edges while the alignment yielded by FW initialized at the true correspondence has
r m_emi$common_edges
edge errors and
r round(m_emi$extra_edges + m_emi$missing_edges)
common edges.
Hence, this graph matching object does not have a solution at the true alignment.
One can try to use other objective functions to enhance the matching result, however we do not investigate this here.
Overall, while most performance measures are poor, our results illustrate the spectrum of challenges for graph matching.
Soft matching: MAP\@3
Considering matching C. Elegans graphs is quite challenging, let's assume 20 pairs of vertices are known as seeds, which are chosen at random. Accordingly, we generate a similarity matrix with 1's corresponding to seeds, and the rest being barycenter.
seeds <- sample(nrow(C1), 20) sim <- init_start(start = "bari", nns = nrow(C1), soft_seeds = seeds)
In addition to one-on-one matching, we will also conduct soft matching, which is to find three most promising matches to each non-seed vertex. We achieve the goal of soft matching by finding the top 3 largest values in each row of the doubly stochastic matrix from the last iteration of Frank Wolfe methodology with indefinite relaxation and PATH algorithm, as well as the normalized matrix from the last iteration of the power method for IsoRank algorithm. To evaluate the matching performance, we will look at both matching precision: $precision=\frac{1}{n_m-s}\sum_{i\in V_m\setminus S}P_{ii}$, and Mean Average Precision \@ 3 (MAP\@ 3):$MAP@3 = \frac{1}{n_m-s}\sum_{i\in V_m\setminus S}1_{{i\in T_i}}$, where $T_i$ is the set of 3 most promising matches to node $i$.
set.seed(123) m_FW <- gm(A = C1, B = C2, seeds = seeds, similarity = sim, method = "indefinite", start = "bari", max_iter = 100) m_PATH <- gm(A = C1, B = C2, seeds = seeds, similarity = NULL, method = "PATH", epsilon = 1, tol = 1e-05) m_Iso <- gm(A = C1, B = C2, seeds = seeds, similarity = as.matrix(sim), method = "IsoRank", max_iter = 50, lap_method = "LAP")
match_eval <- function(match){ precision <- mean(match$corr_A == match$corr_B) order <- apply(match$soft, MARGIN = 1, FUN = order, decreasing = TRUE) top3 <- t(order[1:3,]) - 1:ncol(order) MAP3 <- mean(apply(top3, MARGIN = 1, FUN = function(v){0 %in% v})) round(data.frame(precision, MAP3),4) } sapply(list(m_FW, m_PATH, m_Iso), match_eval) %>% knitr::kable(col.names = c("Frank Wolfe", "PATH", "IsoRank"), booktabs = TRUE, digits = 2)
MAP\@ 3 is slightly higher than precision for each method. Soft matching provides an alternative way of matching by generating a set of promising matching candidates.
Example: Britain Transportation Network {#sec:Transp}
library(Matrix) tm <- Transportation[[1]] cm <- Transportation[[2]] candidate <- Transportation[[3]] tn <- nrow(tm[[1]]) wn <- nrow(cm[[1]]) mc <- candidate %>% with(Matrix::sparseMatrix(i = tem, j = wor, x = 1, dims = c(tn,wn))) c_count <- rowSums(mc)
To demonstrate matching multi-layer networks-layers, we consider two graphs derived from the Britain Transportation network (@transportation).
The network reflects the transportation connections in the UK, with five layers representing ferry, rail, metro, coach, and bus.
A smaller template graph was constructed based on a random walk starting from a randomly chosen hub node, a node that has connections in all the layers.
The template graph has r tn nodes and r sum(sapply(tm, sum)) connections in total and is an induced subgraph of the original graph.
Additionally, based on filter methods from \citet{Moorman2018-ha}, the authors of that paper also provided a list of candidate matches for each template node, where the true correspondence is guaranteed to be among the candidates.
The number of candidates ranges from r min(c_count) to r max(c_count) at most, with an average of r round(mean(c_count)) candidates for each template vertex.
Thus, we made an induced subgraph from the transportation network with only candidates, which gave us the world graph with r wn vertices and r sum(sapply(cm, sum)) connections.
tm <- Transportation[[1]] cm <- Transportation[[2]] candidate <- Transportation[[3]]
Figure \ref{Fig:trans_net} visualizes the transportation connections for the induced subgraphs, where means of transportation are represented by different colors.
Note that all edges in the template are common edges shared by two graphs,
where r paste0(round(100*sum(tm[[1]]) / sum(cm[[1]][1:tn, 1:tn]),1), "%"), r paste0(round(100*sum(tm[[2]]) / sum(cm[[2]][1:tn, 1:tn]),1), "%"), r paste0(round(100*sum(tm[[3]]) / sum(cm[[3]][1:tn, 1:tn]),1), "%"), r paste0(round(100*sum(tm[[4]]) / sum(cm[[4]][1:tn, 1:tn]),1), "%") and r paste0(round(100*sum(tm[[5]]) / sum(cm[[5]][1:tn, 1:tn]),1), "%") of edges in the world graph are in template for each layer.
All graphs are unweighted, directed, and do not have self-loops.
Tables \@ref(tab:edge-summary-trans) further displays an overview and edge summary regarding each layer of the Britain Transportation Network.
A true correspondence exists for each template vertex in the world graph, our goal is to locate each template vertex in the Britain Transportation network by matching two multi-layer graphs with different number of vertices.
summarize_graph_trans <- function(gp1, gp2, name) { nn <- nrow(gp1) if (nn != nrow(gp2)) { nn <- paste(nn, "/", nrow(gp2)) } ne <- paste(sum(gp1), "/", sum(gp2)) data.frame( `Layer` = name, `Nodes` = as.character(nn), `Edges` = ne, `Correlation` = cor(as.numeric(gp1), as.numeric(gp2[1:nrow(gp1), 1:nrow(gp1)])) ) } summary_table <- rbind( summarize_graph_trans(tm[[1]], cm[[1]], "Ferry"), summarize_graph_trans(tm[[2]], cm[[2]], "Rail"), summarize_graph_trans(tm[[3]], cm[[3]], "Metro"), summarize_graph_trans(tm[[4]], cm[[4]], "Coach"), summarize_graph_trans(tm[[5]], cm[[5]], "Bus") ) a <- tm b <- cm %>% purrr::map(~.x[][1:tn, 1:tn]) edge_layer <- summary(graphMatch(data.frame(corr_A = 1:tn, corr_B = 1:tn), as.integer(tn)), a, b)$edge_match_info # edge_layer$layer <- c("Ferry", "Rail", "Metro", "Coach", "Bus") knitr::kable( cbind(summary_table, edge_layer[,2:4]), col.names = c("Layer", "\\# Nodes", "\\# Edges", "Correlation","Common", "Missing", "Extra"), booktabs = TRUE, caption = "Overview of the Britain Transportation Network layers. Correlation is calculted using the template graph and the aligned induced subgraph of the world graph. The final three columns indicate the number of common edges, missing edges, and extra edges in the aligned subgraph of the world graph.", escape = FALSE ) # Edge summary of the Britain Transportation Network layers.
```r Visualization of the template graph (left) and the world graph (right) with corresponding vertices, both derived from the Britain Transportation network with five layers: ferry, rail, metro, coach, and bus. Edges represent transportation transactions and each color indicates a different means of transportation from a different layer of network."}
l <- 1:length(tm)
tg_sub <- l %>% map(~ tm[[.x]] %>% graph_from_adjacency_matrix) wg_sub <- l %>% map(~ cm[[.x]][1:tn,1:tn] %>% graph_from_adjacency_matrix)
for(l in 1:5){ tg_sub[[l]] <- set_edge_attr(tg_sub[[l]], "layer", value = l) wg_sub[[l]] <- set_edge_attr(wg_sub[[l]], "layer", value = l) }
tg_colored <- igraph::union(tg_sub[[1]], tg_sub[[2]], tg_sub[[3]], tg_sub[[4]], tg_sub[[5]]) wg_colored <- igraph::union(wg_sub[[1]], wg_sub[[2]], wg_sub[[3]], wg_sub[[4]], wg_sub[[5]])
tg_layer <- rep(0, ecount(tg_colored)) wg_layer <- rep(0, ecount(wg_colored)) for(l in 1:5){ tg_layer <- tg_layer + ifelse(is.na(edge_attr(tg_colored)[[l]]) | tg_layer != 0, 0, edge_attr(tg_colored)[[l]]) wg_layer <- wg_layer + ifelse(is.na(edge_attr(wg_colored)[[l]]) | wg_layer != 0, 0, edge_attr(wg_colored)[[l]]) }
layout <- layout_nicely(wg_colored) plot(tg_colored, edge.color = factor(tg_layer), vertex.size = 5, vertex.label.cex = .666, edge.curved = TRUE, edge.arrow.size = 0.2, layout = layout)
plot(wg_colored, edge.color = factor(wg_layer), vertex.size = 5, vertex.label.cex = .666, edge.curved = TRUE, edge.arrow.size = 0.2, layout = layout)
```r start <- mc %>% rbind(Matrix(0, nrow = wn - tn, ncol = wn)) start[1:tn, ] <- diag(1 / rowSums(start[1:tn, ])) %*% start[1:tn, ] similarity <- start * 1e5
Based on the candidates, we specify a start matrix that is row-stochastic which can be used for the \code{start} argument in the graph matching function for FW methodology. For each row node, its value is either zero or the inverse of the number of candidates for that node. To ensure that template nodes only get matched to candidates, we constructed a similarity score matrix by taking the start matrix $\times 10^5$, so that a high similarity score is assigned to all the template-candidate pairs.
Then we match the template graph with the world graph using Percolation algorithm. The template graph stored in \code{tm} and world graph \code{cm} are lists of 5 matrices of dimensions 53 and 2075 respectively. Since we have no information on seeds, we assign \code{NULL} to the \code{seeds} argument, the Percolation algorithm will initialize the mark matrix using prior information in the similarity score matrix.
match <- gm(A = tm, B = cm, similarity = similarity, method = "percolation", r = 4) summary(match, tm, cm)
The \code{summary} function outputs edge statistics and objective function values for each layer separately. To further improve matching performance, one can replicate all the analysis in the first example on Enron dataset, such as using the centering scheme and adaptive seeds. Finally, one can refer to the match report to compare matching performance and pick the best one.
Conclusions {#sec:conclusion}
In this work, we detail the methods and usage of the R package iGraphMatch for finding and assessing an alignment between the vertex sets of two edge-correlated graphs.
The package implements common steps for the analysis of graph matching: seamless matching of generalized graphs, evaluation of matching performance, and visualization.
For each of the graph matching methodologies, we provide versatile options for the form of input graphs and the specification of available prior information.
Through the discussion in section \ref{sec:example}, we demonstrate the broad functionality and flexibility of the package by analyzing diverse graph matching problems on real data step by step.
The package also provides tools for simulating correlated graphs which can be used in the development and enhancement of graph matching methods.
Methods for graph matching are still under active development. We plan to include other novel methods as the field continues to develop. In the short term we are looking to introduce a suite of additional matching methods that have recently been proposed in the literature.
One of the biggest challenges for graph matching is evaluating the quality of a match, especially at the vertex level. This has received minimal attention in the previous literature. We provide measures of goodness of matching on the vertex level and demonstrate their effectiveness empirically. These baseline methods implement a permutation testing framework for assessing matches that can be readily extended to other metrics.
\section*{Acknowledgments}
The primary authors for the package are Vince Lyzinski, Zihuan Qiao, and Daniel Sussman. Joshua Agterberg, Lujia Wang, and Yixin Kong also provided important contributions. We also want to thank all of our users, especially Youngser Park, for their feedback and patience as we continue to develop the package.
This work was supported in part by grants from DARPA (FA8750-20-2-1001 and FA8750-18-0035) and from MIT Lincoln Labs.
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.