SEVERAL problems related to the isomorphism andmatching of graphs have been an important and enjoyable challenge for the scientific community for a long time, with applications in pattern recognition (see, for example, [1], [2]), computer vision (see, for example, [3], [4], [5]), and machine learning (see, for example, [6], [7]), to name a few. Given two graphs, the graph isomorphism problem consists of determining whether these graphs are isomorphic or not, that is, if there exists a bijection between the vertex sets of the graphs which exactly preserves the vertex adjacency. The graph isomorphism problem is very challenging from a computational complexity point of view. Indeed, its complexity is still unresolved: it is not currently classified as NPcomplete or P [8]. The graph isomorphism problem is contained in the (harder) graph matching problem. The graph matching problem consists of finding the exact isomorphism between two graphs if it exists, or, in general, finding the bijection between the vertex sets that minimizes the number of adjacency disagreements. Graph matching is a very challenging and well-studied problem in the literature with applications in such diverse fields as pattern recognition, computer vision, neuroscience, etc. (see [9]). Although polynomial-time algorithms for solving the graph matching problem are known for certain classes of graphs (e.g., trees [10], [11]; planar graphs [12]; and graphs with some spectral properties [13], [14]), there are no known polynomial-time algorithms for solving the general case. Indeed, in its most general form, the graph matching problem is equivalent to the NP-hard quadratic assignment problem.

Formally, for any two graphs on n vertices with respective adjacency matrices A and B, the graph matching problem is to minimize over all , where denotes the set of permutation matrices, and is the Froebenius matrix norm (other graph matching objectives have been proposed in the literature as well, this being a common one). Note that for any permutation matrix P, counts the number of adjacency disagreements induced by the vertex bijection corresponding to P.

An equivalent formulation of the graph matching problem is to minimize over all , where is the Euclidean inner product, i.e., for all trace. This can be seen by expanding, for any ,

and noting that and are constants for the optimization problem over .

Let D denote the set of doubly stochastic matrices, i.e., nonnegative matrices with row and column sums each equal to 1. We define the convex relaxed graph matching problem to be minimizing over all , and we define the indefinite relaxed graph matching problem to be minimizing over all . Unlike the graph matching problem, which is an integer programming problem, these relaxed graph matching problems are each continuous optimization problems with a quadratic objective function subject to affine constraints. Since the quadratic objective is also convex in the variables D (it is a composition of a convex function and a linear function), there is a polynomial-time algorithm for exactly solving the convex relaxed graph matching problem (see [15]). However, is not convex (in fact, the Hessian has trace zero and is therefore indefinite), and nonconvex quadratic programming is (in general) NP-hard. Nonetheless the indefinite relaxation can be efficiently approximately solved with Frank-Wolfe (F-W) methodology [16], [17].

It is natural to ask how the (possibly different) solutions to these relaxed formulations relate to the solution of the original graph matching problem. Our main theoretical result, Theorem 1, proves, under mild conditions, that convex relaxed graph matching (which is tractable) almost always yields the wrong matching, and indefinite relaxed graph matching (which is intractable) almost always yields the correct matching. We then illustrate via illuminating simulations that this asymptotic result about the trade-off between tractability and correctness is amply felt even in moderately sized instances.

In light of graph matching complexity results (see for example [14], [18], [19]), it is unsurprising that the convex relaxation can fail to recover the true permutation. In our main theorem, we take this a step further and provide an answer from a probabilistic point of view, showing almost sure failure of the convex relaxation for a very rich and general family of graphs when convexly relaxing the graph matching problem. This paints a sharp contrast to the (surprising) almost sure correctness of the solution of the indefinite relaxation. We further illustrate that our theory gives rise to a new state-of-the-art matching strategy.

1.1 Correlated random Bernoulli graphs

Our theoretical results will be set in the context of correlated random (simple) Bernoulli graphs,1 which can be used to model many real-data scenarios. Random Bernoulli graphs are the most general edge independent random graphs, and contain many important random graph families including Erd˝os-R´enyi and the widely used stochastic block model of [21] (in the stochastic block model, is a block constant matrix, with the number of diagonal blocks representing the number of communities in the network). Stochastic block models, in particular, have been extensively used to model networks with inherent community structure (see, for example, [22], [23], [24], [25]). As this model is a submodel of the random Bernoulli graph model here used, our main theorem (Theorem 1) extends to stochastic block models immediately, making it of highly practical relevance.

These graphs are defined as follows. Given , a real number , and a symmetric, hollow matrix , define where [n] := {1, 2, . . . , n}. Two random graphs with respective adjacency matrices A and B are -correlated Bernoullidistributed if, for all , the random variables (matrix entries) are Bernoullidistributed, and all of these random variables are collectively independent except that, for each , the Pearson product-moment correlation coefficient for is . It is straightforward to show that the parameters , and completely specify the random graph pair distribution, and the distribution may be achieved by first, for all , having Bernoulliindependently drawn and then, conditioning on B, have Bernoulli independently drawn. While would imply the graphs are isomorphic, this model allows for a natural vertex alignment (namely the identity function) for 1, i.e. when the graphs are not necessarily isomorphic.

1.2 The main result

We will consider a sequence of correlated random Bernoulli graphs for n = 1, 2, 3, . . . , where is a function of n. When we say that a sequence of events, , holds almost always we mean that almost surely it happens that the events in the sequence occur for all but finitely many m.

Theorem 1: Suppose A and B are adjacency matrices for -correlated Bernoulligraphs, and there is an such that for all . Let , and denote .

a) If , then it almost always holds that

b) If the between graph correlation , then it almost always holds that

This theorem states that: (part a) the unique solution of the indefinite relaxation almost always is the correct permutation matrix, while (part b) the correct permutation is almost always not a solution of the commonly used convex relation. Moreover, as we will show in the experiments section, the convex relaxation can lead to a doubly stochastic matrix that is not even in the Voronoi cell of the true permutation. In this case, the convex optimum is closest to an incorrect permutation, hence the correct permutation will not be recovered by projecting the doubly stochastic solution back onto .

In the above, and are fixed. However, the proofs follow mutatis mutandis if and are allowed to vary in n. If there exist constants such that and then

Theorem 1, part a will hold. Note that also guarantees the corresponding graphs are almost always connected. For the analogous result for part b, let us first define If there exists an such that then the results of Theorem 1, part b hold as proven below.

1.3 Isomorphic versus ρ-correlated graphs

There are numerous algorithms available in the literature for (approximately) solving the graph isomorphism problem (see, for example, [26], [27]), as well as for (approximately) solving the subgraph isomorphism problem (see, for example, [28]). All of the graph matching algorithms we explore herein can be used for the graph isomorphism problem as well.

We emphasize that the -correlated random graph model extends our random graphs beyond isomorphic graph pairs; indeed -correlated graphs and will almost surely have on the order of edgewise disagreements. As such, these graphs are a.s. not isomorphic. In this setting, the goal of graph matching is to align the vertices across graphs whilst simultaneously preserving the adjacency structure as best possible across graphs. However, this model does preserve a very important feature of isomorphic graphs: namely the presence of a latent alignment function (the identity function in the -correlated model).

We note here that in the -correlated Bernoulli() model, both and are marginally Bernoullirandom graphs, which is amenable to theoretical analysis. We note here that real data experiments across a large variety of data sets (see Section 4.3) and simulated experiments across a variety of robust random graph settings (see Section 4.4) also both support the result of Theorem 1. Indeed, we suspect that an analogue of Theorem 1 holds over a much broader class of random graphs, and we are presently investigating this extension.

Without loss of generality, let . We will first sketch the main argument of the proof, and then we will spend the remainder of the section filling in all necessary details of the proof. The proof will proceed as follows. Almost always, for any such that either or . To accomplish this, we count the entrywise disagreements between AQ and PB in two steps (of course, this is the same as the number of entrywise disagreements between A and ). We first count the entrywise disagreements between B and (Lemma 4), and then count the additional disagreements induced by realizing A conditioning on B. Almost always, this two step realization will result in more errors than simply realizing A directly from B without permuting the vertex labels (Lemma 5). This establishes , and Theorem 1, part a is a consequence of the Birkhoff-von Neumann theorem.

We begin with two lemmas used to prove Theorem 1. First, Lemma 2 is adapted from [29], presented here as a variation of the form found in [30, Prop. 3.2]. This lemma lets us tightly estimate the number of disagreements between B and , which we do in Lemma 4.

Lemma 2: For any integer N > 0 and constant , suppose that the random variable X is a function of at most N independent Bernoulli random variables, each with Bernoulli parameter in the interval . Suppose that changing the value of any one of the Bernoulli random variables (and keeping all of the others fixed) changes the value of X by at most . Then for any t such that , it holds that exp.

The next result, Lemma 3, is a special case of the classical Hoeffding inequality (see, for example, [31]), which we use to tightly bound the number of additional entrywise disagreements between AQ and PB when we realize A conditioning on B.

Lemma 3: Let and be positive integers, and and be real numbers in [0, 1]. If Binomialand Binomialare independent, then for any it holds that

exp

Setting notation for the next lemmas, let n be given. Let denote the set of permutation matrices. Just for now, fix any such that they are not both the identify matrix, and let be their respective associated permutations on [n]; i.e. for all it holds that precisely when and, for all , it holds that precisely when . It will be useful to define the following sets:

If we define m to be the maximum of and , then it follows that 2mn. This is clear from noting that . Also, , since for it is necessary that and . Lastly, , since

and , and ∆ : .

We make the following assumption in all that follows: Assumption 1: Suppose that is a symmetric, hollow matrix, there is a real number , and there is a constant such that for all , and . Further, let A, B be the adjacency matrices of two random -correlated Bernoulligraphs.

Define the (random) set

Note that counts the entrywise disagreements induced within the off-diagonal part of B by and

Lemma 4: Under Assumption 1, if n is sufficiently large then

Proof of Lemma 4: For any , note that has a Bernoulli distribution; if then the Bernoulli parameter is either 0 or is in the interval , and if , then the Bernoulli parameter is and this Bernoulli parameter is in the interval since it is a convex combination of values in this interval. Now, , so we obtain that and thus

Next we apply Lemma 2, since is a function of the at-most N := 2mn Bernoulli random variables , which as a set (noting that is counted at most once for each {i, j}) are independent, each with Bernoulli parameter in . Furthermore, changing the value of any one of these random variable would change by at most , thus Lemma 2 can be applied and, for the choice of , we obtain that

Lemma 4 follows from (1) and (2), since

> αmn/

and when n is sufficiently large (e.g.

With the above bound on the number of (nondiagonal) entrywise disagreements between B and , we next count the number of additional disagreements introduced by realizing A conditioning on B. In Lemma 5, we prove that this two step realization will almost always result in more entrywise errors than simply realizing A from B without permuting the vertex labels.

Lemma 5: Under Assumption 1, it almost always holds that, for all such that either or , .

Proof of Lemma 5: Just for now, let us fix any such that either or , and say and are their respective associated permutations on [n]. Let and

be defined as before. For every , a combinatorial argument, combined with A and B being binary valued, yields (where for an event is the indicator random variable for the event C)

(3)

Note that

Summing Eq. (3) over the relevant indices then yields that

where the sets and are defined as

Now, partition into sets , and partition into sets where

Note that all (i, j) such that i = j are not in . Also note that can be partitioned into the disjoint union .

Equation (4) implies

Now, conditioning on B (hence, conditioning on ), we have, for all , that (see Section 1.1), Bernoulli Thus has a Bernoulli distribution with parameter bounded above by . Thus, is stochastically dominated by a Binomialrandom variable, and the independent random variable is stochastically dominated by a Binomialrandom variable. An application of Lemma 3 with , , and t := , yields (recall that we are conditioning on B here)

exp

2exp

No longer conditioning (broadly) on B, Lemma 4, equations (5) and (6), and , imply that

∥≤ ∥

exp

Until this point, P and Q—and their associated permutations and —have been fixed. Now, for each , define to be the event that for any with the property that their associated permutations are such that the maximum of and is exactly m. There are at mostsuch permutation pairs.

By (7), for every , setting

we have exp expfor some positive constant (the last inequality holding when n is large enough). Thus, for sufficiently large n, expdecays exponentially in n, and is thus finitely summable over n = 1, 2, 3, . . .. Lemma 5 follows from the Borel-Cantelli Lemma.

Proof of Theorem 1, part a: By Lemma 5, it almost always follows that for every not both the identity, . By the Birkhoff-von Neuman Theorem, D is the convex hull of , i.e., for every , there exists constants such that and . Thus, if D is not the identity matrix, then almost always

and almost always argmin

The proof will proceed as follows: we will use Lemma 6 to prove that the identity is almost always not a KKT (Karush-Kuhn-Tucker) point of the relaxed graph matching problem. Since the relaxed graph matching problem is a constrained optimization problem with convex feasible region and affine constraints, this is sufficient for the proof of Theorem 1, part b.

First, we state Lemma 6, a variant of Hoeffding’s inequality, which we use to prove Theorem 1, part b.

Lemma 6: Let N be a positive integer. Suppose that the random variable X is the sum of N independent random variables, each with mean 0 and each taking values in the real interval . Then for any , it holds that

Again, without loss of generality, we may assume I. We first note that the convex relaxed graph matching problem can be written as

where (8) is a convex function (of D) subject to affine constraints (9)-(11) (i.e., ). It follows that if I is the global (or local) optimizer of the convex relaxed graph matching problem, then I must be a KKT (Karush-Kuhn-Tucker) point (see, for example, [32, Chapter 4]).

Hence, a satisfying (9)-(11) (i.e., is primal feasible) is a KKT point if it satisfies

where and are as follows:

noting that the dual variables are not restricted. They correspond to the equality primal constraints (9) that the row-sums of a primal feasible D are all one;

noting that the dual variables are not restricted. They correspond to the equality primal constraints (10) that the column-sums of a primal feasible D are all one;

noting that the dual variables are restricted to be nonnegative. They correspond to the inequality primal constraints (11) that the entries of a primal feasible D be nonnegative. Complementary slackness further constrains the requiring that for all i, j.

At the identity matrix I, the gradient , denoted , simplifies to and I being a KKT point is equivalent to:

where and are as specified above. At the identity matrix, complimentary slackness translates to having .

Now, for Equation (13) to hold, it is necessary that there exist such that

Adding equations (16), (17) and subtracting equations (14), (15), we obtain

Note that , hence Equation (18) is equivalent to (where )

Next, referring back to the joint distribution of A and B (see Section 1.1), we have, for all ,

Now, since

is the sum of Bernoulli random variables which are collectively independent—besides the two of them which are equal, namely and —we have that is stochastically greater than or equal to a Binomialrandom variable. Also note that

is the sum of independent random variables (namely, the ’s) each with mean 0 and each taking on values in . Applying Lemma 3 and Lemma 6, respectively, to and to , with , yields

for some positive constant c (the last inequality holds when n is large enough). Hence the probability that Equation (19) holds is seen to decay exponentially in n, and is finitely summable over n = 1, 2, 3, . . .. Therefore, by the Borel-Cantelli Lemma we have that almost always Equation (19) does not hold. Theorem 1, part b is now shown, since Equation (19) is a necessary condition for .

In the preceding section, we presented a theoretical result exploring the trade-off between tractability and correctness when relaxing the graph matching problem. On one hand, we have an optimistic result (Theorem 1, part a) about an indefinite relaxation of the graph matching problem. However, since the objective function is nonconvex, there is no efficient algorithm known to exactly solve this relaxation. On the other hand, Theorem 1, part b, is a pessimistic result about a commonly used efficiently solvable convex relaxation, which almost always provides an incorrect/non-permutation solution.

TABLE 1. Notation

After solving (approximately or exactly) the relaxed problem, the solution is commonly projected to the nearest permutation matrix. We have not theoretically addressed this projection step yet. It might be that, even though the solution in D is not the correct permutation, it is very close to it, and the projection step fixes this. We will numerically illustrate this not being the case.

We next present simulations that corroborate and illuminate the presented theoretical results, address the projection step, and provide intuition and practical considerations for solving the graph matching problem. Our simulated graphs have n = 150 vertices and follow the Bernoulli model described above, where the entries of the matrix are i.i.d. uniformly distributed in with . In each simulation, we run 100 Monte Carlo replicates for each value of . Note that given this value, the threshold in order to fulfill the hypothesis of the first part of Theorem 1 (namely that is . As in Theorem 1, for a fixed , we let , so that the correct vertex alignment between and B is provided by the permutation matrix .

We then highlight the applicability of our theory and simulations in a series of real data examples. In the first set of experiments, we match three pairs of graphs with known latent alignment functions. We then explore the applicability of our theory in matching graphs without a pre-specified latent alignment. Specifically, we match 16 benchmark problems (those used in [17], [33]) from the QAPLIB library of [34]. See Section 4.3 for more detail. As expected by the theory, in all of our examples a smartly initialized local minimum of the indefinite relaxation achieves best performance.

We summarize the notation we employ in Table 1. To find , we employ the F-W algorithm ([16], [33]), run to convergence, to exactly solve the convex relaxation. We also use the Hungarian algorithm ([35]) to compute , the projection of to . To find a local minimum of , we use the FAQ algorithm of [17]. We use FAQ:, FAQ:, and FAQ:J to denote the FAQ algorithm initialized at , and (the barycenter of D). We compare our results to the GLAG and PATH algorithms, implemented with off-the-shelf code provided by the algorithms’ authors. We restrict our focus to these algorithms (indeed, there are a multitude of graph matching algorithms present in the literature) as these are the prominent relaxation algorithms; i.e., they all first relax the graph matching problem, solve the relaxation, and then project the solution onto .

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0

Fig. 1: For , we plot (red /green) and (blue/gray). Red/blue dots correspond to simulations where , and grey/green dots to . Black dots correspond to . For each we ran 100 MC replicates.

4.1 On the convex relaxed graph matching problem

Theorem 1, part b, states that we cannot, in general, expect . However, is often projected onto , which could potentially recover . Unfortunately, this projection step suffers from the same problems as rounding steps in many integer programming solvers, namely that the distance from the best interior solution to the best feasible solution is not well understood.

In Figure 1, we plot versus the cor- relation between the random graphs, with 100 replicates per value of . Each experiment produces a pair of dots, either a red/blue pair or a green/grey pair. The energy levels corresponding to the red/green dots correspond to , while the energies corresponding to the blue/grey dots correspond . The colors indicate whether was (green/grey pair) or was not (red/blue pair) . The black dots correspond to the values of .

Note that, for correlations , as expected from Theorem 1, part b. Also note that, even for correlations greater than , we note after projecting to the closest permutation matrix, even though with high probability is the solution to the unrelaxed problem.

We note the large gap between the pre/post projection energy levels when the algorithm fails/succeeds in recovering , the fast decay in this energy (around in Figure 1), and the fact that the value for can be easily predicted from the correla- tion value. These together suggest that can be used a posteriori to assess whether or not graph matching recovered . This is especially true if is known or can be estimated.

How far is from ? When the graphs are isomorphic (i.e., in our setting), then for a large class of graphs, with certain spectral constraints, then is the unique solution of the convex relaxed graph matching

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0

Fig. 2: Distance from to (in blue), to (in red), and to a random permutation (in black). For each value of , we ran 100 MC replicates.

problem [14]. Indeed, in Figure 1, when we see that as expected. On the other hand, we know from Theorem 1, part b that if it is often the case that . We may think that, via a continuity argument, if the correlation is very close to one, then will be very close to , and will probably recover .

We empirically explore this phenomena in Figure 2. For , with 100 MC replicates for each , we plot the (Frobenius) distances from to (in blue), from to (in red), and from to a uniformly random permutation in (in black). Note that all three distances are very similar for , implying that is very close to the barycenter and far from the boundary of D. With this in mind, it is not surprising that the projection fails to recover for in Figure 1, as at the barycenter, the projection onto is uniformly random.

For very high correlation values (), the distances to and to sharply decrease, and the distance to a random permutation sharply increases. This suggests that at these high correlation levels moves away from the barycenter and towards . Indeed, in Figure 1 we see for that is the closest permutation to , and is typically recovered by the projection step.

4.2 On indefinite relaxed graph matching problem

The continuous problem one would like to solve, (since its optimum is with high probability), is indefinite. One option is to look for a local minimum of the objective function, as done in the FAQ algorithm of [17]. The FAQ algorithm uses F-W methodology ([16]) to find a local minimum of . Not surprisingly (as there are many local minima), the performance of the algorithm is heavily dependent on the initialization. Below we study the effect of initializing the algorithm at the non-informative barycenter, at (a principled starting point), and at . We then compare performance of the different FAQ initializations to the PATH algorithm [33] and to the GLAG algorithm [36].

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0

Fig. 3: Success rate in recovering . In gray, FAQ starting at, from left to right, , and J; in black, ; in red, PATH; in blue, GLAG. For each we ran 100 MC replicates.

The GLAG algorithm presents an alternate formulation of the graph matching problem. The algorithm convexly relaxes the alternate formulation, solves the relaxation and projects it onto . As demonstrated in [36], the algorithm’s main advantage is in matching weighted graphs and multimodal graphs. The PATH algorithm begins by finding , and then solves a sequence of concave and convex problems in order to improve the solution. The PATH algorithm can be viewed as an alternative way of projecting onto . Together with FAQ, these algorithms achieve the current best performance in matching a large variety of graphs (see [36], [17], [33]). However, we note that GLAG and PATH often have significantly longer running times than FAQ (even if computing for FAQ:); see [17], [37].

Figure 3 shows the success rate of the graph matching methodologies in recovering . The vertical dashed red line at corresponds to the threshold in Theorem 1 part a (above which is optimal whp) for the parameters used in these experiments, and the solid lines correspond to the performance of the different methods: from left to right in gray, FAQ:, FAQ:, FAQ:J; in black, the success rate of ; the performance of GLAG and PATH are plotted in blue and red respectively.

Observe that, when initializing with , the fact that FAQ succeeds in recovering means that is a local minimum, and the algorithm did not move from the initial point. From the theoretical results, this was expected for , and the experimental results show that this is also often true for smaller values of . However, this only means that is a local minimum, and the function could have a different global minimum. On the other hand, for very lowly correlated graphs (is not even a local minimum.

The difference in the performance illustrated by the gray lines indicates that the resultant graph matching solution can be improved by using as an initialization to find a local minimum of the indefinite relaxed problem. We see in the figure that FAQ:achieves

Fig. 4: Average run time for FAQ:(note that this does not include the time to find ) and FAQ:J in gray; finding (first finding ) in black; PATH in red; and GLAG in blue. For each , we average over 100 MC replicates. Note that the runtime of PATH drop precipitously at which corresponds to the performance increase in Figure 3.

best performance, while being computationally less intensive than PATH and GLAG, see Figure 4 for the runtime result. This amalgam of the convex and indefi-nite methodologies (initialize indefinite with the convex solution) is an important tool for obtaining solutions to graph matching problems, providing a computationally tractable algorithm with state-of-the-art performance.

However, for all the algorithms there is still room for improvement. In these experiments, for theory guarantees that with high probability the global minimum of the indefinite problem is , and we cannot find it with the available methods.

When FAQ:fails to recover , how close is the objective function at the obtained local minima to the objective function at ? Figure 5 shows for the true permutation, , and for the pre-projection doubly stochastic local minimum found by FAQ:. For , the state-of-the-art algorithm not only fails to recover the correct bijection, but also the value of the objective function is relatively far from the optimal one. There is a transition (around ) where the algorithm moves from getting a wrong local minimum to obtaining (without projection!). For low values of , the objective function values are very close, suggesting that both and the pre-projection FAQ solution are far from the true global minima. At we see a separation between the two objective function values (agreeing with the findings in Figure 3). As , we expect that is the global minima and the pre-projection FAQ solution is far from until the phase transition at .

4.3 Real data experiments

We further demonstrate the applicability of our theory in a series of real data examples. First we match three pairs of graphs where a latent alignment is known. We

Local and global minima of the indefinite relaxation

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -1.7

Fig. 5: Value of for (black) and for the output of FAQ:(red/blue indicating failure/success in recovering the true permutation). For each we ran 100 MC replicates.

further compare different graph matching approaches on a set of 16 benchmark problems (those used in [17], [33]) from the QAPLIB QAP library of [34], where no latent alignment is known a priori. Across all of our examples, an intelligently initialized local solution of the indefinite relaxation achieves best performance.

Our first example is from human connectomics. For 45 healthy patients, we have DT-MRI scans from one of two different medical centers: 21 patients scanned (twice) at the Kennedy Krieger Institute (KKI), and 24 patients scanned (once) at the Nathan Kline Institute (NKI) (all data available at http://openconnecto.me/ data/public/MR/MIGRAINE v1 0/). Each scan is identically processed via the MIGRAINE pipeline of [38] yielding a 70 vertex weighted symmetric graph. In the graphs, vertices correspond to regions in the Desikan brain atlas, which provides the latent alignment of the vertices. Edge weights count the number of neural fiber bundles connecting the regions. We first average the graphs within each medical center and then match the averaged graphs across centers.

For our second example, the graphs consist of the twohop neighborhoods of the “Algebraic Geometry” page in the French and English Wikipedia graphs. The 1382 vertices correspond to Wikipedia pages with (undirected) edges representing hyperlinks between the pages. Page subject provides the latent alignment function, and to make the graphs of commensurate size we match the intersection graphs.

Lastly, we match the chemical and electrical connectomes of the C. elegans worm. The connectomes consist of 253 vertices, each representing a specific neuron (the same neuron in each graph). Weighted edges representing the strength of the (electrical or chemical) connection between neurons. Additionally, the electrical graph is directed while the chemical graph is not.

The results of these experiments are summarized in Table 2. In each example, the computationally inexpensive FAQ:procedure achieves the best performance compared to the more computationally expensive GLAG

TABLE 2. for the P given by each al- gorithm together with the number of vertices correctly matched () in real data experiments

and PATH procedures. This reinforces the theoretical and simulation results presented earlier, and again points to the practical utility of our amalgamated approach. While there is a canonical alignment in each example, the results point to the potential use of our proposed procedure (FAQ:) for measuring the strength of this alignment, i.e., measuring the strength of the correlation between the graphs. If the graphs are strongly aligned, as in the KKI-NKI example, the performance of FAQ:will be close to the truth and a large portion of the latent alignment with be recovered. As the alignment is weaker, FAQ:will perform even better than the true alignment, and the true alignment will be poorly recovered, as we see in the C. elegans example.

What implications do our results have in graph matching problems without a natural latent alignment? To test this, we matched 16 particularly difficult examples from the QAPLIB library of [34]. We choose these particular examples, because they were previously used in [17], [33] to assess and demonstrate the effectiveness of their respective matching procedures. Results are summarized in Table 3. We see that in every example, the indefinite relaxation (suitably initialized) obtains the best possible result. Although there is no latent alignment here, if we view the best possible alignment as the “true” alignment here, then this is indeed suggested by our theory and simulations. As the FAQ procedure is computationally fast (even initializing FAQ at both J and is often comparatively faster than GLAG and PATH; see [17] and [37]), these results further point to the applicability of our theory. Once again, theory suggests, and experiments confirm, that approximately solving the indefinite relaxation yields the best matching results.

4.4 Other random graph models

While the random Bernoulli graph model is the most general edge-independent random graph model, in this section we present analogous experiments for a wider variety of edge-dependent random graph models. For these models, we are unaware of a simple way to exploit pairwise edge correlation in the generation of these graphs, as was present in Section 1.1. Here, to simulate aligned non-isomorphic random graphs, we

TABLE 3. for the different tested algo- rithms on 16 benchmark examples of the QAPLIB library.

proceed as follows. We generate a graph from the appropriate underlying distribution, and then model as an errorful version of ; i.e., for each edge in , we randomly flip the edge (i.e., bit-flip from or ) independently with probability . We then graph match and , and we plot the performance of the algorithms in recovering the latent alignment function across a range of values of p.

We first evaluate the performance of our algorithms on power law random graphs [39]; these graphs have a degree distribution that follows a power law, i.e., the proportion of vertices of degree d is proportional to for some constant . These graphs have been used to model many real data networks, from the Internet [40], [41], to social and biological networks [42], to name a few. In general, these graphs have only a few vertices with high degree, and the great majority of the vertices have relatively low degree.

Figure 6 shows the performance comparison for the methods analyzed above: FAQ:, FAQ:, FAQ:, PATH, and GLAG. For a range of , we generated a 150 vertex power law graph with , and subsequently graph matched this graph and its errorful version. For each p, we have 100 MC replicates. As with the random Bernoulli graphs, we see from Figure 6 that the true permutation is a local minimum of the non-convex formulation for a wide range of flipping probabilities (), implying that in this range of and share significant common structure. Across all values of p < 0.5, FAQ:outperforms all other algorithms considered (with FAQ:being second best across this range). This echoes the results of Sections (4.1)–(4.3), and suggests an analogue of Theorem 1 may hold in the power law setting. We are presently investigating this.

We next evaluate the performance of our algorithms on graphs with bounded maximum degree (also called bounded valence graphs). These graphs have been extensively studied in the literature, and for bounded valence graphs, the graph isomorphism problem is in P [43]. For the experiments in this paper we generate a random graph from the model in [44] with maximum degree equal to 4, and vary the graph order from 50 to 350

Fig. 6: Success rate in recovering for 150 vertex power law graphs with for: In gray, from right to left, FAQ:, FAQ:, and FAQ:J; in black, ; in red, PATH; in blue, GLAG. For each value of the bit-flip parameter p, we ran 100 MC replicates.

vertices. Figure 7 shows the comparison of the different techniques and initializations for these graphs, across a range of bit-flipping parameters .

It can be observed that even for isomorphic graphs (p = 0), all but FAQ:fail to perfectly recover the true alignment. We did not see this phenomena in the other random graph models, and this can be explained as follows. It is a well known fact that convex relaxations fail for regular graphs [13], and also that the bounded degree model tends to generate almost regular graphs [45]. Therefore, even without flipped edges, the graph matching problem with the original graphs is very illconditioned for relaxation techniques. Nevertheless, the true alignment is a local minimum of the non-convex formulation for a wide range of values of p (shown by FAQ:performing perfectly over a range of p in Figure 7). We again note that FAQ:outperforms , PATH and GLAG across all graph sizes and bit-flip parameters p. This suggests that a variant of Theorem 1 may also hold for bounded valence graphs as well, and we are presently exploring this.

We did not include experiments with any random graph models that are highly regular and symmetric (for example, mesh graphs). Symmetry and regularity have two effects on the graph matching problem. Firstly, it is well known that for non-isomorphic regular graphs (indeed, J is a solution of the convex relaxed graph matching problem). Secondly, the symmetry of these graphs means that there are potentially several isomorphisms between a graph and its vertex permuted analogue. Hence, any flipped edge could make permutations other than into the minima of the graph matching problem.

4.5 Directed graphs

All the theory developed above is proven in the undirected graph setting (i.e., A and B are assumed symmetric). However, directed graphs are common in numerous

Fig. 7: Success rate in recovering for bounded degree graphs (max degree 4). In gray, from right to left, FAQ:, FAQ:, and FAQ:J; in black, ; in red, PATH; in blue, GLAG. For each probability we ran 100 MC replicates.

applications. Figure 8 repeats the analysis of Figure 3 with directed graphs, all other simulation parameters being unchanged. The PATH algorithm is not shown in this new figure because it is designed for undirected graphs, and its performance for directed graphs is very poor. Recall that in Figure 3, i.e., in the undirected setting, FAQ:J performed significantly worse than . In Figure 8, i.e., the directed setting, we note that the performance of FAQ:J outperforms over a range of ]. As in the undirected case, we again see significant performance improvement (over FAQ:, and GLAG) when starting FAQ from (the convex solution). Indeed, we suspect that a directed analogue of Theorem 1 holds, which would explain the performance increase achieved by the nonconvex relaxation over . Here, we note that the setting for the remainder of the examples considered is the undirected graphs setting.

4.6 Seeded graphs

In some applications it is common to have some a priori information about partial vertex correspondences, and seeded graph matching includes these known partial matchings as constraints in the optimization (see [46], [47], [14]). However, seeds do more than just reducing the number of unknowns in the alignment of the vertices. Even a few seeds can dramatically increase performance graph matching performance, and (in the -correlated Erd˝os-R´enyi setting) a logarithmic (in n) number of seeds contain enough signal in their seed–to– nonseed adjacency structure to a.s. perfectly align two graphs [47]. Also, as shown in the deterministic graph setting in [14], very often is closer to .

In Figure 9, the graphs are generated from the -correlated random Bernoulli model with random (entrywise uniform over [0.1, 0.9]). We run the Frank-Wolfe method (modified to incorporate the seeds) to solve the convex relaxed graph matching problem, and the method in [46], [47] to approximately solve the non-

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0

Fig. 8: Success rate for directed graphs. We plot (black), the GLAG method (blue), and the nonconvex relaxation starting from different points in green, from right to left: FAQ:J, FAQ:, FAQ:.

convex relaxation, starting from , and . Note that with seeds, perfect matching is achieved even below the theoretical bound on provided in Theorem 1 (for ensuring is the global minimizer). This provides a potential way to improve the theoretical bound on in Theorem 1, and the extension of Theorem 1 for graphs with seeds is the subject of future research. With the exception of the nonconvex relaxation starting from , each of the different FAQ initializations and the convex formulation all see significantly improved performance as the number of seeds increases. We also observe that the nonconvex relaxation seems to benefit much more from seeds than the convex relaxation. Indeed, when comparing the performance with no seeds, the performs better than FAQ:J. However, with just five seeds, this behavior is inverted. Also of note, in cases when seeding returns the correct permutation, we’ve empirically observed that merely initializing the FAQ algorithm with the seeded start, and not enforcing the seeding constraint, also yields the correct permutation as its solution (not shown).

Figure 10 shows the running time (to obtain a solution) when starting from for the nonconvex relaxation, using different numbers of seeds. For a fixed seed level, the running time is remarkably stable across when FAQ does not recover the true permutation. On the other hand, when FAQ does recover the correct permutation, the algorithm runs significantly faster than when it fails to recover the truth. This suggests that, across all seed levels, the running time might, by itself, be a good indicator of whether the algorithm succeeded in recovering the underlying correspondence or not. Also note that as seeds increase, the overall speed of convergence of the algorithm decreases and, unsurprisingly, the correct permutation is obtained for lower correlation levels.

4.7 Features

Features are additional information that can be utilized to improve performance in graph matching methods,

Fig. 9: Success rate of different methods using seeds. We plot (top left), FAQ:J (top right), FAQ:(bottom left), and FAQ:(bottom right). For each method, the number of seeds increases from right to left: 0 (black), 5 (green), 10 (blue) and 15 (red) seeds. Note that more seeds increases the success rate across the board.

and often these features are manifested as additional vertex characteristics besides the connections with other vertices. For instance, in social networks we may have have a complete profile of a person in addition to his/her social connections.

We demonstrate the utility of using features with the nonconvex relaxation, the standard convex relaxation and the GLAG method, duely modified to include the features into the optimization. Namely, the new objective function to minimize is tracewhere F(P) is the original cost function (in the nonconvex setting, for the convex relaxation and for the GLAG method), the matrix C codes the features fitness cost, and the parameter balances the trade-off between pure graph matching and fit in the features domain. For each of the matching methodologies, the optimization is very similar to the original featureless version.

For the experiments, we generate -correlated Bernoulli graphs as before, and in addition we generate a Gaussian random vector (zero mean, unit variance) of 5 features for each node of one graph, forming a matrix of features; we permute that matrix according to to align new features vectors with the nodes of the second graph. Lastly, additive zero-mean Gaussian noise with a range of variance values is added to each feature matrix independently. If for each vertex the resulting noisy feature for , is , then the entries of C are defined to be for . Lastly, we set

Figure 11 shows the behavior of the methods when using features for different levels of noise in the feature matrix. Even for highly noisy features (recalling that both feature matrices are contaminated with noise), this external information still helps in the graph matching problem. For all noise levels, all three methods improve their performance with the addition of features, and of

Fig. 10: Running time for the nonconvex relaxation when starting from , for different number of seeds. A red “x” indicates the algorithm failed to recover , and a black “o” indicates it succeeded. In each, the algorithm was run to termination at discovery of a local min.

course, the improvement is greater when the noise level decreases. Note that, as before, FAQ outperforms both and GLAG across all noise levels. It is also worth noting that for low noise, FAQ:performs comparably to FAQ:, which we did not observe in the seeded (or unseeded) setting.

Even for modestly errorful features, including these features improves downstream matching performance versus the setting without features. This points to the utility of high fidelity features in the matching task. Indeed, given that the state-of-the-art graph matching algorithms may not achieve the optimal matching for even modestly correlated graphs, the use of external information like seeds and features can be critical.

In this work we presented theoretical results showing the surprising fact that the indefinite relaxation (if solved exactly) obtains the optimal solution to the graph matching problem with high probability, under mild conditions. Conversely, we also present the novel result that the popular convex relaxation of graph matching almost always fails to find the correct (and optimal) permutation. In spite of the apparently negative statements presented here, these results have an immediate practical implication: the utility of intelligently initializing the indefinite matching algorithm to obtain a good approximate solution of the indefinite problem.

The experimental results further emphasize the trade-off between tractability and correctness in relaxing the graph matching problem, with real data experiments and simulations in non edge-independent random graph models suggesting that our theory could be extended to more general random graph settings. Indeed, all of our experiments corroborate that best results are obtained via approximately solving the intractable indefinite problem. Additionally, both theory and examples point to

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0

Fig. 11: Success rate of different methods using features: (in black), GLAG (in blue), FAQ:(in red), and FAQ:(in green). For each method, the noise level (variance of the Gaussian random noise) increases from left to right: 0.3, 0.5, and 0.7. In dashed lines, we show the success of the same methods without features.

the utility of combining the convex and indefinite approaches, using the convex to initialize the indefinite.

[1] T. Caelli and S. Kosinov, “An eigenspace projection clustering method for inexact graph matching,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 26, no. 4, pp. 515–519, 2004.

[2] A. C. Berg, T. L. Berg, and J. Malik, “Shape matching and object recognition using low distortion correspondences,” in 2005 IEEE Conference on Computer Vision and Pattern Recognition, 2005, pp. 26–33.

[3] B. Xiao, E. R. Hancock, and R. C. Wilson, “A generative model for graph matching and embedding,” Computer Vision and Image Understanding, vol. 113, no. 7, pp. 777–789, 2009.

[4] M. Cho and K. M. Lee, “Progressive graph matching: Making a move of graphs via probabilistic voting,” in 2012 IEEE Conference on Computer Vision and Pattern Recognition, 2012, pp. 398–405.

[5] F. Zhou and F. De la Torre, “Factorized graph matching,” in 2012 IEEE Conference on Computer Vision and Pattern Recognition, 2012, pp. 127–134.

[6] B. Huet, A. D. J. Cross, and E. R. Hancock, “Graph matching for shape retrieval,” in Advances in Neural Information Processing Systems, 1999, pp. 896–902.

[7] T. Cour, P. Srinivasan, and J. Shi, “Balanced graph matching,” in Advances in Neural Information Processing Systems, vol. 19. MIT; 1998, 2007, pp. 313–320.

[8] M. Garey and D. Johnson, Computers and Intractability: A Guide to the Theory of NP-completeness. W.H. Freeman, 1979.

[9] D. Conte, P. Foggia, C. Sansone, and M. Vento, “Thirty years of graph matching in pattern recognition,” International Journal of Pattern Recognition and Artificial Intelligence, vol. 18, no. 03, pp. 265–298, 2004.

[10] A. Torsello, D. Hidovic-Rowe, and M. Pelillo, “Polynomial-time metrics for attributed trees,” Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol. 27, no. 7, pp. 1087–1099, 2005.

[11] J. D. Ullman, A. V. Aho, and J. E. Hopcroft, “The design and analysis of computer algorithms,” Addison-Wesley, Reading, vol. 4, pp. 1–2, 1974.

[12] J. E. Hopcroft and J.-K. Wong, “Linear time algorithm for isomor- phism of planar graphs (preliminary report),” in Proceedings of the sixth annual ACM symposium on Theory of computing. ACM, 1974, pp. 172–184.

[13] M. Fiori and G. Sapiro, “On spectral properties for graph matching and graph isomorphism problems,” arXiv preprint arXiv:1409.6806, 2014.

[14] Y. Aflalo, A. Bronstein, and R. Kimmel, “Graph matching: Relax or not?” arXiv:1401.7623, 2014.

[15] D. Goldfarb and S. Liu, “An oprimal interior point algorithm for convex quadratic programming,” Mathematical Programming, vol. 49, no. 1-3, pp. 325–340, 1990.

[16] M. Frank and P. Wolfe, “An algorithm for quadratic program- ming,” Naval Research Logistics Quarterly, vol. 3, no. 1-2, pp. 95– 110, 1956.

[17] J. Vogelstein, J. Conroy, V. Lyzinski, L. Podrazik, S. Kratzer, E. Harley, D. Fishkind, R. Vogelstein, and C. Priebe, “Fast approximate quadratic programming for graph matching,” arXiv:1112.5507, 2012.

[18] R. O’Donnell, J. Wright, C. Wu, and Y. Zhou, “Hardness of robust graph isomorphism, Lasserre gaps, and asymmetry of random graphs,” arXiv:1401.2436, 2014.

[19] A. Atserias and E. Maneva, “Sherali–Adams relaxations and in- distinguishability in counting logics,” SIAM Journal on Computing, vol. 42, no. 1, pp. 112–137, 2013.

[20] B. Bollob´as, Random Graphs. Springer, 1998.

[21] P. W. Holland, K. B. Laskey, and S. Leinhardt, “Stochastic block- models: First steps,” Social Networks, vol. 5, no. 2, pp. 109–137, 1983.

[22] T. A. Snijders and K. Nowicki, “Estimation and prediction for stochastic blockmodels for graphs with latent block structure,” Journal of classification, vol. 14, no. 1, pp. 75–100, 1997.

[23] K. Nowicki and T. A. B. Snijders, “Estimation and prediction for stochastic blockstructures,” Journal of the American Statistical Association, vol. 96, no. 455, pp. 1077–1087, 2001.

[24] M. E. Newman and M. Girvan, “Finding and evaluating commu- nity structure in networks,” Physical review E, vol. 69, no. 2, p. 026113, 2004.

[25] E. M. Airoldi, D. M. Blei, S. E. Fienberg, and E. P. Xing, “Mixed membership stochastic blockmodels,” in Advances in Neural Information Processing Systems, 2009, pp. 33–40.

[26] L. P. Cordella, P. Foggia, C. Sansone, and M. Vento, “A (sub) graph isomorphism algorithm for matching large graphs,” Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol. 26, no. 10, pp. 1367–1372, 2004.

[27] S. Fankhauser, K. Riesen, H. Bunke, and P. Dickinson, “Subopti- mal graph isomorphism using bipartite matching,” International Journal of Pattern Recognition and Artificial Intelligence, vol. 26, no. 06, 2012.

[28] J. R. Ullmann, “An algorithm for subgraph isomorphism,” Journal of the ACM (JACM), vol. 23, no. 1, pp. 31–42, 1976.

[29] N. Alon, J. Kim, and J. Spencer, “Nearly perfect matchings in regular simple hypergraphs,” Israel Journal of Mathematics, vol. 100, pp. 171–187, 1997.

[30] J. H. Kim, B. Sudakov, and V. H. Vu, “On the asymmetry of random regular graphs and random graphs,” Random Structures and Algorithms, vol. 21, pp. 216–224, 2002.

[31] F. Chung and L. Lu, “Concentration inequalities and Martingale inequalities: A survey,” Internet Mathematics, vol. 3, pp. 79–127, 2006.

[32] M. Bazaraa, S. Mokhtar, H. Sherali, and C. M. Shetty, Nonlinear Programming: Theory and Algorithms. John Wiley & Sons, 2013.

[33] M. Zaslavskiy, F. Bach, and J. Vert, “A path following algorithm for the graph matching problem,” Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol. 31, no. 12, pp. 2227–2242, 2009.

[34] R. E. Burkard, S. E. Karisch, and F. Rendl, “Qaplib–a quadratic assignment problem library,” Journal of Global Optimization, vol. 10, no. 4, pp. 391–403, 1997.

[35] H. W. Kuhn, “The Hungarian method for the assignment prob- lem,” Naval Research Logistic Quarterly, vol. 2, pp. 83–97, 1955.

[36] M. Fiori, P. Sprechmann, J. Vogelstein, P. Mus´e, and G. Sapiro, “Robust multimodal graph matching: Sparse coding meets graph matching,” Advances in Neural Information Processing Systems 26, pp. 127–135, 2013.

[37] V. Lyzinski, D. L. Sussman, D. E. Fishkind, H. Pao, L. Chen, J. T. Vogelstein, Y. Park, and C. E. Priebe, “Spectral clustering for divide-and-conquer graph matching,” stat, vol. 1050, p. 22, 2014.

[38] W. R. Gray, J. A. Bogovic, J. T. Vogelstein, B. A. Landman, J. L. Prince, and R. J. Vogelstein, “Magnetic resonance connectome automated pipeline: an overview,” Pulse, IEEE, vol. 3, no. 2, pp. 42–48, 2012.

[39] A.-L. Barab´asi and R. Albert, “Emergence of scaling in random networks,” Science, vol. 286, no. 5439, pp. 509–512, 1999.

[40] R. Albert, H. Jeong, and A.-L. Barab´asi, “Internet: Diameter of the world-wide web,” Nature, vol. 401, no. 6749, pp. 130–131, 1999.

[41] M. Faloutsos, P. Faloutsos, and C. Faloutsos, “On power-law rela- tionships of the internet topology,” in ACM SIGCOMM Computer Communication Review, vol. 29, no. 4. ACM, 1999, pp. 251–262.

[42] M. Girvan and M. E. Newman, “Community structure in social and biological networks,” Proceedings of the National Academy of Sciences, vol. 99, no. 12, pp. 7821–7826, 2002.

[43] E. M. Luks, “Isomorphism of graphs of bounded valence can be tested in polynomial time,” Journal of Computer and System Sciences, vol. 25, no. 1, pp. 42–65, 1982.

[44] K. Balinska and L. Quintas, “Algorithms for the random f- graph process,” Communications in Mathematical and in Computer Chemistry/MATCH, no. 44, pp. 319–333, 2001.

[45] V. Koponen, “Random graphs with bounded maximum degree: asymptotic structure and a logical limit law,” arXiv preprint arXiv:1204.2446, 2012.

[46] D. E. Fishkind, S. Adali, and C. E. Priebe, “Seeded graph match- ing,” arXiv:1209.0367, 2012.

[47] V. Lyzinski, D. E. Fishkind, and C. E. Priebe, “Seeded graph matching for correlated Erdos-Renyi graphs,” Journal of Machine Learning Research, vol. 15, pp. 3513–3540, 2014. [Online]. Available: http://jmlr.org/papers/v15/lyzinski14a.html