arXiv:1803.02423v1 [stat.ML] 6 Mar 2018

Matched Filters for Noisy Induced Subgraph
Detection
Daniel L. Sussman† , Vince Lyzinski‡ , Youngser Park∗ , and
Carey E. Priebe∗
†

Department of Mathematics and Statistics, Boston University
of Mathematics, University of Massachusetts Amherst
∗ Department of Applied Mathematics and Statistics, Johns Hopkins University
‡ Department

March 8, 2018
Abstract
We consider the problem of finding the vertex correspondence between two graphs with different number of vertices where the smaller
graph is still potentially large. We propose a solution to this problem
via a graph matching matched filter: padding the smaller graph in
different ways and then using graph matching methods to align it to
the larger network. Under a statistical model for correlated pairs of
graphs, which yields a noisy copy of the small graph within the larger
graph, the resulting optimization problem can be guaranteed to recover the true vertex correspondence between the networks, though
there are currently no efficient algorithms for solving this problem. We
consider an approach that exploits a partially known correspondence
and show via varied simulations and applications to the Drosophila
connectome that in practice this approach can achieve good performance.
Keywords: multiple graph inference, subgraph detection, graph matching

1

1

Introduction

In many settings, we often want to quantify how multiple networks relate to
each other in order to study how actors jointly use these networks. This may
arise from multiple modalities, such as communications networks, delivery
networks, financial networks, and social networks, or from a time dynamic
setting. Similarly, in neuroscience or biology we may seek to compare brain
networks or protein networks of different individuals or species. Often, these
networks are on different unmatched sets of vertices that are not the same
cardinality. This limits the set of available tools as the adjacency matrices
and other graph statistics may not be directly comparable. In a related fashion, detection and location of particular induced subgraphs that correspond
to a given activity or structure of interest are also paramount problems in
various domains.
When the induced subgraphs are small, with up to ≈10–20 vertices, there
are a number of approaches including color coding [32, 5] and backtracking
algorithms [18]. Parallel algorithms have also been proposed to approximately enumerate the number of isomorphic copies of a subgraph appearing
in a large network [26]. However, it may be the case that these subgraphs
may contain  20 vertices. In this case, many existing approaches will either fail to find the subgraphs of interest or be computationally intractable.
While state space representations [8] have been used to identify these larger
subgraphs, these also have limitations in terms of the applicable graph structure. As a simple example of the challenge of this problem, consider finding
or enumerating cliques. Finding cliques and estimating the number of cliques
of small sizes can be achieved via fast algorithms [25], but the problem of
find the maximum clique is a notoriously challenging problem [3].
Furthermore, while even the subgraph isomorphism is itself NP hard [27],
we often expect the subgraphs may not appear exactly in the graph but
rather only approximately. This may be due to errors and missingness, in
the observation of one or both of the graphs, so finding an exact subgraph
might not be possible.
We view both of these as different points on a continuum of related problems that are all examples of graph matching [7, 15]. In a simpler setting,
when the two graphs adjacency matrices A, B ∈ {0, 1}n×n are of equal size,
the graph matching problem (GMP) seeks the permutation matrix
argmin kA − P BP T k = argmax trace(AP BP T ),
P ∈P

P ∈P

2

(1)

where P is the set of n × n permutation matrices. While this problem is
NP-hard in general, there are numerous approaches in the literature that
attempt to approximately solve the GMP. Conte et al. [7] and EmmertStreib et al. [12] provide reviews of the prescient literature. In particular,
when prior knowledge about the correspondence between vertices the can
be incorporated into the algorithms, the GMP can be approximately solved
efficiently for graphs with more than 105 vertices [21, 29] without the need
for sophisticated modern parallel computing to be brought to bear.
A more challenging problem that we will not consider is to detect anomalous subgraphs within a collection of graphs [2]. In the anomaly detection
setting, the structure of the anomalous subgraph may be only known up
to certain graph characteristics or deviations from the structure of the remaining graph. In our setting, the approximate topology of the graph to be
detected is known explicitly.
To approximately solve the noisy induced subgraph problem, we will use
the machinery of graph matching to construct a graph matching matched
filter, aligning the smaller network to its closest subgraph within the larger
network. Our main goals in this manuscript are to investigate the theoretical
limitations of this noisy graph matching framework when the graphs may be
of very different sizes. To match graphs of radically different sizes in Eq. (1),
we consider a number of padding schemes for the smaller matrix to render
graph matching an appropriate tool for the problem [14]. Under a statistical
model for noisily implanting an induced subgraph into a larger network,
we show that the true induced subgraph will be found by an oracle GM
algorithm using an appropriate padding scheme provided graph correlations
and probabilities satisfy certain parameter assumptions.
We further demonstrate the effectiveness of these strategies when the
vertex correspondence is partially known. By relaxing the graph matching
problem constraints to the set of doubly stochastic matrices, a gradient descent approach can be put to use [28]. The problem is non-convex [22] and
hence the initialization can incorporate the partial correspondence and random restarts to seek multiple local minimums. By choosing the optimal of
these local minima, an exact or approximate recovery of the vertex correspondence can often be found.

3

2

Background for Graph Matching

In this section we provide a brief background on graph matching, some methods to incorporate prior information, and a statistical model for correlated
graphs. Throughout the remainder of this article we will use the following
notation.
Let [n] = {1, 2, . . . , n}. Let Pn and Dn denote the set of n × n permutation matrices and doubly stochastic matrices, respectively. Let Jn and 0n
denote the n × n all-ones and all-zeros matrices, respectively, and Jnm and
0nm for n × m rectangular versions. Let An denote the set of adjacency matrices corresponding to simple undirected graphs, i.e. the set of symmetric
hollow {0, 1}-valued n × n matrices. When clear from context, we may omit
subscripts indicating the sizes of these matrices. Finally, let ⊕ denote the
M 0
direct sum of two matrices, M ⊕ M 0 =
.
0 M0

2.1

Algorithms

Solving the GMP problem exactly is NP-hard in general, but there are a
number of approaches which provide good approximate solutions [7, 12].
Some approaches rely on tree based methods for exact branch and bound
algorithms for integer programs [16]. While these approaches are guaranteed to find a global minimum, their computational scaling is very poor. For
larger graphs, the constraints for the GMP are often relaxed so that continuous optimization machinery can be deployed [28, 31, 13]. The relaxed
solutions are then projected back onto P yielding an approximate solution
to the graph matching problem. Some relaxations involve applying spectral
methods which allows the application of fast linear algebra techniques [10].
A convex approach has also been proposed but requires the graph to have
certain properties [1], which do not frequently hold in practice.
For the computational experiments in this manuscript we will rely on
a principled indefinite relaxation of the GMP constraints [22] to the set of
doubly stochastic matrices D, the convex hull of P. Details of this algorithm
are discussed in Section 3.2. This approach has been shown to perform
very well in practice on a number of benchmark problems, especially when
exploiting multiple random restarts to avoid poor performing local minima
[28]. We expect that many of the phenomena observed in our experiments
would appear in similar fashions for other algorithms.
4

Frequently, these approaches can exploit seeded vertices, i.e., a partial list
of known correspondences between the two vertex sets. When sufficiently
many seeds are present, these algorithms, which often have few guarantees
without seeds, can be solved efficiently and the resulting match can be guaranteed to be correct [14, 20] asymptotically almost surely for relatively general random graph models. While the theory we discuss below does not
require seeds, our algorithms will use seeds and the algorithmic performance
relies heavily on a sufficient number of a priori available seeds.
In this work, we employ a soft seeding approach which uses the prior information to initialize the algorithm but does not insist that the correspondence
for seeds be fixed; i.e., if the seeded vertices are [m] in both networks, we initialize the gradient descent graph matching algorithm at Im ⊕ D for suitable
D in Dn−m . Conversely, in the hard-seeded setting we optimize over Πn−m
with the global solution being then of the form Im ⊕ P for P ∈ Πn−m . This
allows for the possibility that seeds may be noisy, allowing for the algorithm
to potentially correct for incorrect initial seeds, a scenario that we do not
explore further here.

2.2

Statistical Models

In order to study the applicability and limitations of a graph matching approach for subgraph detections, we analyze the problem from a statistical
viewpoint in the setting of the correlated heterogeneous Erdős-Rényi model
[22]. The following definition provides a distribution for pairs of random
graphs with different numbers of vertices where there is a positive correlation between corresponding edge-pairs. Without loss of generality, we assume
that the smaller graph (of order nc ) is denoted A and that the corresponding
vertices to A in B are the first nc ≤ n vertices.
Definition 1 (Correlated Erdős-Rényi). Suppose Λ ∈ [0, 1]n×n and R ∈
[0, 1]nc ×nc for 0 < nc ≤ n. Denote by Λc , the order nc principal submatrix
of Λ. A pair of adjacency matrices (A, B) ∼ CorrER(Λ, R) if A ∈ Anc ,
B ∈ An , and for each u < v, Buv are independent with Buv ∼ Bernoulli(Λuv )
and for each u < v ≤ nc , Auv are independent with Auv ∼ Bernoulli(Λuv ).
Additionally, Bu,v and Au0 ,v0 are independent except that for u, v ∈ [nc ], with
u < v, it holds that the Pearson correlation corr(Auv , Buv ) = Ruv .
When nc = n, it can be shown that the solution to the GMP will asymptotically almost surely yield the correct vertex correspondence under mild
5

assumptions on Λ and R, i.e., the unique element in the argmin in Eq. (1) is
the identity matrix I [20, 22].
Note that other authors have considered equivalent reparameterized versions of the model above. In particular, a common viewpoint is that each observed graph is a random subgraph of some common underlying graph, which
also arises from independent edge sampling [29]. By varying the parameters
for the underlying graph and the parameters for the random subgraph, any
correlated Erdős-Rényi can be represented in this way, provided all correlations are non-negative. The reverse is also true, and more distributions
can be constructed by allowing for negative correlations and non-identical
distributions. Another closely related possibility is that one graph is a noisy
observation of an induced subgraph from the other graph, where there may
be extra or missing edges [22]. This can also be related to the correlated
Erdős-Rényi by an appropriate reparameterization.

3

Padding Approaches

In order to match pairs of graphs with differing numbers of vertices, we
propose to pad the smaller graph/matrix with enough rows and columns to
match the order of the larger matrix. We will consider a trio of padding
schemes which will result in differing qualities for the resulting match [14].
e = A ⊕ 0n and match A
e and B.
Naive Padding Let A
j
Centered Padding
e = (2A − J) ⊕ 0n , let B
e = 2B − J, and match A
e and B.
e
Let A
j
Oracle Padding
e = (A − Λnc ) ⊕ 0, let B
e = B − Λ, and match A
e and B.
e
Let A
As we will see in the next section, the naive padding scheme—which finds
the best fitting subgraph of B to match with A—will not be guaranteed
to find the true correspondence between nodes. In contrast, the centered
padding scheme—which finds the best fitting induced subgraph of B to match
with A—and the oracle padding scheme are guaranteed to succeed under mild
model conditions, even in the presence of an exponentially small (in terms
of the size of B) subgraph A. In general, the oracle padding scheme will be

6

inaccessible as Λ is unknown, but using various methods to estimate Λ [6, 9],
we can approximate Λ in ways that can improve matching performance.
Note that each of these padding approaches corresponds to a particular
graph detection problem. With the naive padding scheme,
P the optimization program is equivalent to the problem argmaxσ:[nc ][n] u,v Aσ(u)σ(v) Buv ,
where  denotes that σ is injective. Hence, the global optimum reveals the
largest common subgraph of B and A, as only edge agreements are rewarded
whereas edge disagreements are not penalized. Note that if B has a large
clique, then the naive approach will tend to match A to this large clique
as all edges in A are present in the large clique. As we will see below, the
naive padding approach can be guaranteed with high probability to yield an
incorrect vertex correspondence under only mild assumptions on R and Λ.
If the graphs have the same number of nodes, both the centered padding
scheme and the naive padding scheme are equivalent [4], they are not equivalent for graphs of different orders. With the centered padding scheme,
the opP
timization program is equivalent to the problem argminσ:[nc ][n] u,v |Aσ(u)σ(v) −
Buv |. In this case, the global optimum reveals the induced subgraph of B
which is nearest to A in terms of edge-edit distance [4]. As will be seen below, the optimum to this program can be guaranteed to correspond the true
vertex correspondence provided that the graphs are not too sparse and the
correlation is sufficiently large.
The oracle padding scheme corresponds to the problem,
X
d σ(u)σ(v) Buv ),
argmin
Cov(A
σ:[nc ][n] u,v

d denotes the empirical covariance when both means are known,
where Cov
d ij , Buv ) = (Aij − Λij )(Buv − Λuv ). As this approach eliminates any
as Cov(A
undue rewards from matching likely edges to even more likely edges or unlikely edges to even more unlikely edges, the global optimum is guaranteed
with high probability to be the true correspondence even under mild assumptions on Λ and R. Indeed, even if A and B were non-identically distributed,
with means ΛA and ΛB respectively, by matching A − ΛA and B − ΛB , the
theoretical results will still hold.
While the theory below shows that under the correlated Erdős-Rényi
model the oracle padding yields stronger guarantees then the centered padding
method and the naive padding method can have very bad performance, this
will not necessarily translate in practice. This can occur especially due to
7

the fact that we do not find the true global optimum of the GMP but rather
a local optimum near our initialization. Our simulations and real data examples illustrate how centering or naive may be preferable in certain scenarios,
and biased estimates of the oracle may also improve performance. In practice, one may try different paddings and even use the solution of one padding
approach to initialize the optimization for another padding approach.

3.1

Theory

For each of the padding scenarios, we will consider nc and n as tending to ∞
in order to understand the ability of these optimization programs to detect
progressively larger subgraphs. Note that we will require that the number
of vertices nc is growing with n. Indeed, if nc is fixed and n grows then
eventually every subgraph of size nc will likely appear multiple times as an
induced subgraph in B just by chance [19].
For each padding scenario let P ∗ ∈ {0, 1}nc ×nc denote the order nc principal submatrix of the solution to the corresponding graph matching problem.
The proofs of each theorem below can be found in Appendix B
The first theorem is a negative result that shows that under weak assumptions on Λc , one can construct a Λ under which the naive padding scheme is
almost surely guaranteed to not detect the errorful version of A in B. Indeed,
it can be guaranteed to recover zero true vertex with high probability.
Theorem 2. Suppose R ∈ [0, 1]nc ×nc satisfies R < ρ < 1 entry-wise. Suppose
nc < n/2, β ∈ (0, 1), and entry-wise Λc ≤ β and Λc = ω(n−1
c log nc ). Then
there exists Λ such that for A, B ∼ CorrER(Λ, R), using the naive padding
scheme,

P[P ∗ 6= I] ≥ 1 − exp −C2 (log nc )2 ,
(2)
for some universal constant C > 0 and  ∈ (0, 1 − β − (1 − β)ρ).
This result hinges on the fact that the naive padding scheme finds the best
matching subgraph, rather than the best matching induced subgraph; indeed,
there is no penalty for matching non-edges in A to edges in B. Hence, if B
has a dense substructure of size nc which does not correspond to A then
the naive padding scheme will match A to that dense structure with high
probability, regardless of the underlying correlation. The proof demonstrates
an extreme case where all vertices in A are matched (incorrectly) to vertices in
B without true matches. Deviations on the proof will yield that any number
8

of vertices may be mismatched with only a mildly denser substructure in B.
Recall that since the centered padding scheme and the naive padding scheme
are equivalent when n = nc , the requirement that nc is sufficiently small is
necessary.
On the other hand, the centered padding scheme can be guaranteed to
result in the correct detection of the subgraph even when the number of
vertices in B is exponentially larger than the number of vertices in A provided
R > 1/2 + .
Theorem 3. Suppose that A, B ∼ CorrER(Λ, R) with Ruv ∈ [1/2 + , 1] and
= o(1) implies
Λuv ∈ [α, 1 − α] for α,  ∈ (0, 1/2). It holds that 2 nclog(n)
α(1−α)2
that using the centered padding scheme

P[P ∗ 6= I] ≤ 2 exp −Θ(nc (α(1 − α))2 ) .
Stated simply, this theorem only requires that the correlations are sufficiently
large to guarantee that large subgraphs of logarithmic size can be found via
an oracle GM algorithm.
Finally, while the oracle padding is inaccessible for general Λ, it represents the optimal padding scheme as it eliminates any empirical correlations
introduced by Λ leaving only the theoretical correlations from R.
Theorem 4. Suppose that A, B ∼ CorrER(Λ, R) with Ruv ∈ [ρ, 1] and Λuv ∈
[α, 1 − α], for some α, ρ ∈ (0, 1). Then, using the centered padding scheme,

P[P ∗ 6= I] ≤ 2 exp −Θ(nc (ρα(1 − α))2 ) .
Remark 5. Note that each while the padding schemes are still not equivalent
in the homogeneous Erdős-Rényi case, the distribution of the matrix P ∗ are
equal under all three padding schemes and hence all are guaranteed to find
the correct solution with high probability under identical scenarios.

3.2

Computation

Our approach to solve the graph matching problem in this setting will be
analogous to the approach for graphs of equal sizes; see the FAQ algorithm of
[28] for full detail. In particular, we will relax the GMP constraints from P to
D and use gradient descent starting at a given D0 . We will also incorporate
seeds, which, without loss of generality, we will assume are the first s nodes.
This gradient ascent approach is then given by Algorithm 1.
9

Data: A, B ∈ A, D0 ∈ D, k = 0
while not converged do
e k BP
e );
1
Pk ← argminP ∈P −trace(AD
e α BD
e α ), where
2
αk ← argminα∈[0,1] −trace(AD
Dα = αDk + (1 − α)Pk ;
3
Dk+1 ← Dαk and k ← k + 1;
end
∗
4 Project Dk onto P yielding P ;
∗
5 Return P
Algorithm 1: Fast approximate quadratic assign program (FAQ) algorithm [28] for graph matching.
Note that using any of the padding schemes, we do not need to store or
compute the entire matrices Dk or Pk as the objective function only depends
on their first nc rows. Hence, lines 1 and 4 can be simplified and accomplished
by searching over the set of nc × n matrices corresponding to injections from
[nc ] to [n], or equivalently the first nc rows of permutation matrices. In this
way, lines 1 and 4 can be solved effectively by variants of the Hungarian
algorithm for non-square matrices [23]. Line 2 is a quadratic equation in α
and is easily solved. Furthermore, by exploiting the fact that Ã and B̃ are
likely formed by the difference between a sparse and low-rank matrix, fast
matrix multiplication and storage can also be exploited for further runtime
improvements.
The convergence criterion is generally easy to check as the optimal doubly stochastic matrix is frequently itself a permutation matrix, which also
means the final projection step can be omitted. While this algorithm is not
guaranteed to converge to a global optimum, if there are enough seeds or if
the matrix D0 is sufficiently close to the identity, the local maximum which
this procedure converges to will likely be the identity.

4

Experiments

In this section we demonstrate the effectiveness of a graph matching matched
filter for errorful subgraph detection in both synthetic and real data settings.
Without loss of generality, we will assume that the seeds are the first s vertices
for some s ≤ nc . The matched filter algorithm is given in Algorithm 2
10

1

2

Data: Template A, network B, number of seeds s, number of MC
replicates M , padding scheme
Result: Matchings {B1 , B2 , · · · , BM }
eB
e ← padded and centered A, B according to padding scheme;
A,
for i ← 1 to M do
Generate a random doubly stochastic matrix


1−α
D0 ← Is ⊕ αP +
Jn−s ,
n−s

where P is a random permutation matrix and α ∼ Unif(0, 1);
e to B;
e
3
Initialize FAQ at D0 ; match A
∗
4
FAQ output is P , assignment of vertices in B to vertices in A
end
Algorithm 2: Matched filters algorithm via random restarts, soft seeding,
and gradient descent.
As described in Line 2, we initialize the FAQ algorithm at a random start
point which initially aligns the seeded vertices across graphs. In particular,
our random starting point has the identity as its principal s × s submatrix,
corresponding to the seeds. The remaining matrix is sampled as a random
1
Jn−s and a random perconvex combination of the the baricenter matrix n−s
mutation matrix sampled uniformly from Pn−s . Repeating this process M
times, we output M potential matches for the subgraph A in B.
To explore the effectiveness of this procedure, we consider applying our
matched filter in the following scenarios: pairs of homogeneous CorrER
graphs where we vary the uniform correlation, a planted partition model
where the goal is to find the dense partition, a heterogeneous model with two
different subgraph sampling schemes, and an application to the Drosophila
connectome.

4.1

Correlated Homogeneous Erdős-Rényi

As a synthetic data example, we will consider subgraph detection in the
CorrER(λJ, ρJ) with λ = 0.5 (i.e., the maximum entropy ER model) for ρ =
0.7, 0.8, 0.9, and 1. When ρ = 1, the induced subgraph of B corresponding to
A is exactly isomorphic to A. We consider n = 500, nc = 40, and M = 200
11

correlation=0.7

correlation=0.8

40

30

correctly matched vertices

20

10

seeds

0

0
correlation=0.9

correlation=1

7

40

15

30

20

10

0
0

50

100

150

200

0

50

100

150

200

objective function rank

Figure 1: We plot the average number of matched vertices, out of 40, which
were correctly matched according to the rank of the objective function. The
different colors correspond to using either 0, 7, or 15 soft seeds and each panel
corresponds to the homogeneous Erdős-Rényi, as described in Section 4.1.
and we replicated the entire experiment for 200 sampled graph pairs. For
each experiment, we rank the objective function values among the M = 200
random restarts, from smallest to largest.
Figure 1 shows the average performance across the replicates at each
rank. As is evident, without seeds very few vertices are matched correctly,
even when the corresponding graph are isomorphic. However, as both the
number of seeds increase and as the correlations increase, the performance
of the procedure becomes very good.
Interpreting the correlation 0.9 panel, if 15 seeds are used then nearly all
restarts will yield a perfect recovery of the original correspondence. Hence,
only one restart is required. Using 7 seeds, most restarts yield suboptimal

12

correlation=0.7

correlation=0.8

34

33

32

30

seeds
30

0

27

objective function

7
15
24

28
correlation=0.9

count

correlation=1

1

35

10
30

50

30

100
200

20
25
10
20
0
0

10

20

30

40

0

10

20

30

40

correctly matched vertices

Figure 2: Each point denotes the conditional mean of the objective function
given the number of matched vertices. Each point is scaled according the
average number of restarts that achieve the given number of matched vertices.
Note that for correlations ≥ 0.8, there is a large gap in the expected objective
function value for when there are ≈ 40 correct matches as compared to < 40.
performance, with only slightly more than the original 7 seeds recovered.
However, if we perform 200 restarts and choose the best among them then
this nearly always will yield a perfect recovery of the original correspondence.
Similar interpretations can be made for the other correlations.
Figure 2 gives a complementary view of this same simulation, showing
the average value of the objective function as a function of the number of
matched vertices. Each point is scaled according to the average number of
restarts (out of 200) which achieve the given number of matched vertices.
For ρ = 0.7, the objective function decreases smoothly with the number of
matched vertices, and the number of matched vertices is relatively evenly

13

seeds ρ = 0.7
0
7
15

4.90 ± 1.86
4.71 ± 1.78
3.62 ± 1.54

ρ = 0.8

ρ = 0.9

ρ=1

4.89 ± 1.84
4.64 ± 1.74
2.33 ± 1.18

4.90 ± 1.86
4.47 ± 1.75
1.58 ± 0.738

4.86 ± 1.83
3.99 ± 1.82
1.13 ± 0.497

Table 1: Runtimes per restart for each setting in terms of the number of
seeds and the correlation ρ in the homogeneous Erdős-Rényi setting.
spread, especially for 15 seeds. On the other hand, for higher correlations
and for 15 seeds, there is a substantial gap in both the objective function
and the matched vertices.
This gap is highly important for practical applications. Indeed, the presence of a gap in the objective function would indicate that the correct matching was likely found. Without the gap, the guarantees that the best match in
terms of the objective function corresponds to a correct match is less certain
(At least to practically ascertain). This gap also allows for a fast and adaptive algorithm which performs random restarts until a sufficient objective
function gap is observed.
Finally, table one shows average run times, on a standard laptop, for each
restart under each parameterization. Even under the worst of the settings,
200 restarts can be performed in around 15 minutes. It is worth noting
that better algorithmic performance seems to correspond to faster algorithmic runtime, a feature which was observed in previous work using the FAQ
algorithm [14].

4.2

Planted Partition Correlated Erdős-Rényi graphs

As a another scenario, we consider the problem where the larger graph has a
denser subgraph which corresponds to the smaller graph. We take n = 500,
nc = 40, A, B ∼ CorrER(Λ, 0.9J) with


qJnc
pJn−nc ,n
Λ=
,
pJn,n−nc pJn−nc
where p = 0.25 and q ∈ {0.25, 0.3, . . . , 0.5}. For 100 random restarts, 7
seeds and using the centered padding, we replicated each setting 200 times.
Figure 3 shows the average number of correctly matched at each objective
14

correctly matched vertices

40

30

q
0.50
0.45
0.40
0.35

20

0.30
0.25

10

0

25

50

75

100

objective function rank

Figure 3: Number of correctly matched vertices as a function of the objective
function rank for the planted partition model described in Section 4.2. Each
line corresponds to a the density of the induced subgraph, ranging from
q = 0.25 to 0.5.
function rank for each value of q. As would be expected, the problem becomes
significantly easier as q increases. A similar gap phenomenon is observed in
this setting, as long as q is sufficiently large, though we do not illustrate it
here.

4.3

Random Dot Product Graph

In order to better evaluate the performance of the procedure under heterogeneous graph models, we will now consider graphs from the random dot
product graph model [30, 24]. Each vertex v is associated with a latent position Xv ∈ R2 , which is sampled uniformly from the set {(x, y) ∈ (0, 1)2 :
x+y < 1}. In this case the matrix Λ is given by Λuv = XuT Xv . We again take
n = 500 and nc = 40 and sample A, B ∼ CorrER(Λ, 0.9J). The matched
subgraphs are sampled in two different ways. Under one sampling design, the
matched nodes simply correspond to a random sample of the latent positions.
Under the “max-angle” design, the matched nodes are selected to be those
with the largest ratio between the second and first elements of the latent po15

1.00

●
●

●
●
●
●

●

0.75

●

●
●
●
● ●
●● ●
●
●
●●
●
●

which

x2

●●
●● ●

●

●

0.50

max angle
non−matched

●
●
●

random

●●

0.25
●
●
●
●
●

0.00
0.00

0.25

0.50

0.75

1.00

x1

Figure 4: A representative set of latent positions for the random dot product
graph model. The red circles correspond to the max-angle sampling and
the points in blue correspond to a possible random sampling (these do not
overlap for the purposes of visualization).
sitions. Figure 4 shows a representative set of latent positions corresponding
to this model. The red circles correspond to the max-angle sampling and
the points in blue correspond to a possible random sampling (these do not
overlap for the purposes of visualization).
We sampled 200 distinct graphs under each subgraph sampling setting
and performed 100 random restarts for each sampled graphs. For each
restart, the matchings were performed with the naive, centered, and two
“oracle”-type padding schemes. As the model above is a low-rank model, we
centered each matrix according to its best rank-1 and rank-2 approximations
to approximate the oracle padding scheme. Letting,
Λ̂(A)
= argmin kA − Λk2F
r

and Λ̂(B)
= argmin kB − Λk2F ,
r

Λ:rank(Λ)=r

Λ:rank(Λ)=r

these schemes match (A − Λ̂(A) ) ⊕ 0n−nc to B − Λ̂(B) for r = 1, 2.
Figure 5 shows the average accuracy at each rank of the objective function. First, the max-angle subgraph is clearly much easier to discover than
the random subgraph. Interestingly, for the max-angle case, the centered
16

max angle

random

correctly matched vertices

40

30

scheme
center
naive

20

rank 1
rank 2
10

0
0

25

50

75

100

0

25

50

75

100

objective function rank

Figure 5: Matching accuracy
padding has by far the best performance, while in the random case, the rankone centering performs the best. Rank-two centering, which corresponds to
the true model dimension, has the third best performance in either case.
Naive padding fails to find any signal.
Under either scenario, the rank-one approximation eliminates most degree
based differences between vertices which will often be sufficient to eliminate
excessive empirical correlations. Rank-two centering however may actually
hurt performance because it eliminates distinctions between vertices in terms
of which nodes they tend to be adjacent to. Note that constant centering
maintains all differential vertex behavior which is a strong indicator of which
vertices to match to. The naive padding fails here as there are always sets
of vertices with the nearly same relative propensities for adjacencies to other
vertices but with strictly larger absolute probabilities of adjacencies.
To understand why the max-angle problem is so much easier, note that,
under either subgraph sampling scheme, there will be many nodes which
have similar behavior to one or more of the sampled nodes due to having
nearby latent positions. Under random sampling, most of these similar nodes
are likely to correspond to nodes that do not have a match and hence will
introduce errors. However, under max-angle sampling the similar nodes have
a better chance of being other matchable nodes which will introduce fewer
errors as the theoretical correlation will strongly distinguish these vertices.
17

max_angle

random
0.8

0.8

objective function

scheme
0.6

0.6

1
2
center

0.4
0.4

count
1
0.2

10
0.2

0.0
0

10

20

30

40

0

10

20

30

40

correctly matched vertices

Figure 6: RDPG objective function gap
Figure 6 shows the average value of the objective function as a function
of the number of matched vertices. In order to more easily compare the objective function values for different centering schemes, we rescaled objectives
functions for each scheme so that the largest for each was equal one and the
smallest was zero. Each point is sized according to the average number of
restarts (out of 100) which achieve the given number of matched vertices.
Note that we again observe a gap in the objective function between the bulk
of solutions with less than 20 correct matches and those with more than 35
correct matches.

4.4

Finding KC cells in a Drosophila connectome

As an application to brain networks, we consider using the matched filter to
locate an induced subgraph of Kenyon cells (KC) in the fully reconstructed
Drosophila mushroom body of [11]. Using the induced subgraph of the KC
cells in the left hemisphere of the mushroom body (i.e., as A), we seek to
find the KC cells on the right hemisphere. Although in this example, the
KC cells are identified across both hemispheres, this was achieved only with
great effort and expenditure. Being able to use one hemisphere to locate
structure in the other hemisphere could potentially allow for faster, cheaper
neuron identification in future connectomes.
18

20
50

40
100
60
150
80

200
50

100

150

200

20

40

60

80

Figure 7: Left: The 209 node graph corresponding to the left hemisphere of
the Drosophila mushroom body. The first 100 nodes (rows/columns) are the
KC cells. Right: The 100 nodes graph of KC cells in the right hemisphere of
the Drosophila mushroom body.
After initial data preprocessing, there are nc = 100 KC cells in each
hemisphere and n = 209 vertices total in the right hemisphere. The pair
of graphs are illustrated in Figure 7. Note that as the true correspondence
between these pair of graphs is unknown, we will evaluate the accuracy in
terms of the number of cells in the right KC graph which are matched to KC
cells in the left graph, rather than the number of correctly matched pairs. We
consider s = 0, 2, 4, 8 seeds and M = 1000 random restarts for our matched
filter. The seeds are random pairs of KC cells across the two graphs.
In the top panel of Figure 8, we plot the performance of our filter with the
centered padding scheme, and in the bottom panel we use an approximation
to the oracle padding scheme, centering by the least-squares optimal rank-1
approximations of Λ and Λc . In the figure, we see that more seeds produces
better performance, and that the approximate oracle centering provides better performance (finding ≈ 85 of the 100 KC cells in the right hemisphere)
than the data-independent centered padding scheme. We expect progressively more accurate estimations of the Λ matrices in the oracle centering to
produce even better results, ideally resulting in an optimization gap as in the
synthetic setting.

19

Figure 8: Detecting KC cells in the Drosophila mushroom body. In the top
panel, we plot the performance of our filter with the centered padding scheme,
and in the bottom panel we use an approximation to the oracle padding
scheme, centering by the least-squares optimal rank-1 approximations of Λ
and Λc . In the first row of each panel, we plot the GMP objective function
for each initialization (note that we ordered the random restarts in the figure
to be monotone in f ); in the second row, we plot the number of vertices
correctly matched. In the third row of the top panel, we plot the number
of edges correctly matched. The number of correctly matched edges are not
reported for the oracle as the objective function does not correspond to edge
errors. Different colors indicate different numbers of seeds.

20

5

Discussion

In this manuscript we have identified a method for noisy subgraph detection using a modification of the graph matching problem. We have proposed a number of padding methods, each of which emphasize a different
aspects of common graph structure. Theoretically, the centered and oracle
padding schemes are guaranteed to recover the original vertex correspondence in the correlated Erdős-Rényi statistical model, provided correlations
and edge probabilities do not decay too quickly and the order of the smaller
graph is at least logarithmic in the order of the larger graph.
By using a partially known vertex correspondence and random restarts,
we are frequently able to recover the full correspondence accurately in a
number of synthetic data settings and in a connectomics setting. Importantly,
we frequently observe a gap in the objective function value between a restart
that achieves good performance and those with poor performance. If this gap
is present, it can used to detect the fact that the correct alignment has been
found rather than a spurious local minimum, and additionally can provide
an adaptive stopping criterion for performing random restarts.
In this manuscript, we have assumed that the induced subgraphs are identically distributed and positively correlated. These assumptions are important for our proof. Indeed, if negative or correlation is allowed, the objective
functions are seeking the wrong solution. Non-identical graphs can be more
easily dealt with by using separate oracle estimates, as was done in the random dot product and Drosophila connectome examples. Further theoretical
issues such as non-independent edge graph may also be possible to deal with
by exploiting the possibility that the graph pair is constructed out of a still
large number of independent Bernoulli trials, none of which can drastically
change graph structure. Such an assumption may allow for the application
of Proposition 6 as in our proofs below.

Acknowledgments
This work was partially supported by a grant from MIT Lincoln Labs and
the Department of Defense, the MAA program of the Defense Advanced
Research Projects Agency (DARPA) via award number FA8750-18-0035.

21

References
[1] Y. Aflalo, A. Bronstein, and R. Kimmel. Graph matching: relax or not?
arXiv:1401.7623, 2014.
[2] Leman Akoglu, Hanghang Tong, and Danai Koutra. Graph based
anomaly detection and description: a survey. Data mining and knowledge discovery, 29(3):626–688, May 2015.
[3] Immanuel M Bomze, Marco Budinich, Panos M Pardalos, and Marcello
Pelillo. The maximum clique problem. In Ding-Zhu Du and Panos M
Pardalos, editors, Handbook of Combinatorial Optimization: Supplement
Volume A, pages 1–74. Springer US, Boston, MA, 1999.
[4] H Bunke. On a relation between graph edit distance and maximum
common subgraph. Pattern recognition letters, 18(8):689–694, August
1997.
[5] V T Chakaravarthy, M Kapralov, P Murali, F Petrini, X Que, Y Sabharwal, and B Schieber. Subgraph counting: Color coding beyond trees.
In 2016 IEEE International Parallel and Distributed Processing Symposium (IPDPS), pages 2–11, May 2016.
[6] S. Chatterjee. Matrix estimation by universal singular value thresholding. The Annals of Statistics, 43(1):177–214, 2014.
[7] 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, 18(03):265–298, 2004.
[8] Luigi P Cordella, Pasquale Foggia, Carlo Sansone, and Mario Vento.
A (sub)graph isomorphism algorithm for matching large graphs. IEEE
transactions on pattern analysis and machine intelligence, 26(10):1367–
1372, October 2004.
[9] Mark A Davenport, Yaniv Plan, Ewout van den Berg, and Mary Wootters. 1-bit matrix completion. Information and Inference: A Journal of
the IMA, 3(3):189–223, September 2014.
[10] Amir Egozi, Yosi Keller, and Hugo Guterman. A probabilistic approach
to spectral graph matching. IEEE transactions on pattern analysis and
machine intelligence, 35(1):18–27, January 2013.
22

[11] K. Eichler, F. Li, A. Litwin-Kumar, Y. Park, I. Andrade, C. M.
Schneider-Mizell, T. Saumweber, A. Huser, C. Eschbach, B. Gerber,
et al. The complete connectome of a learning and memory centre in an
insect brain. Nature, 548(7666), 2017.
[12] F. Emmert-Streib, M. Dehmer, and Y. Shi. Fifty years of graph matching, network alignment and network comparison. Information sciences,
346–347:180–197, 2016.
[13] M. Fiori, P. Sprechmann, J. Vogelstein, P. Mus, and G. Sapiro. Robust
multimodal graph matching: Sparse coding meets graph matching. Advances in Neural Information Processing Systems, pages 127–135, 2013.
[14] D.E. Fishkind, S. Adali, H. G. Patsolic, L. Meng, V. Lyzinski, and C.E.
Priebe. Seeded graph matching. arXiv preprint arXiv:1209.0367, 2017.
[15] P. Foggia, G. Percannella, and M. Vento. Graph matching and learning
in pattern recognition in the last 10 years. International Journal of
Pattern Recognition and Artificial Intelligence, 28(01):1450001, 2014.
[16] D E Ghahraman, A K C Wong, and T Au. Graph optimal monomorphism algorithms. IEEE transactions on systems, man, and cybernetics,
10(4):181–188, April 1980.
[17] J. H. Kim, B. Sudakov, and V. H. Vu. On the asymmetry of random
regular graphs and random graphs. Random Structures and Algorithms,
21:216–224, 2002.
[18] M Kuramochi and G Karypis. Frequent subgraph discovery. In Proceedings 2001 IEEE International Conference on Data Mining, pages
313–320, 2001.
[19] László Lovász. Large networks and graph limits, volume 60. American
Mathematical Society Providence, 2012.
[20] V. Lyzinski, D. E. Fishkind, and C. E. Priebe. Seeded graph matching for correlated Erdos-Renyi graphs. Journal of Machine Learning
Research, 15:3513–3540, 2014. URL http://jmlr.org/papers/v15/
lyzinski14a.html.

23

[21] V. Lyzinski, D. L. Sussman, D. E. Fishkind, H. Pao, L. Chen, J. T.
Vogelstein, Y. Park, and C. E. Priebe. Spectral clustering for divideand-conquer graph matching. Parallel Computing, 47:70–87, 2015.
[22] V Lyzinski, D E Fishkind, M Fiori, J T Vogelstein, C E Priebe, and
G Sapiro. Graph matching: Relax at your own risk. IEEE transactions
on pattern analysis and machine intelligence, 38(1):60–73, 2016.
[23] James Munkres. Algorithms for the assignment and transportation problems. Journal of the Society for Industrial and Applied Mathematics, 5
(1):32–38, 1957.
[24] C. L. M. Nickel. Random dot product graphs: A model for social networks. PhD thesis, Johns Hopkins University, 2006.
[25] Matthew C Schmidt, Nagiza F Samatova, Kevin Thomas, and ByungHoon Park. A scalable, parallel algorithm for maximal clique enumeration. Journal of parallel and distributed computing, 69(4):417–428, April
2009.
[26] G M Slota and K Madduri. Fast approximate subgraph counting and
enumeration. In 2013 42nd International Conference on Parallel Processing, pages 210–219, October 2013.
[27] J R Ullmann. An algorithm for subgraph isomorphism. Journal of the
ACM, 23(1):31–42, January 1976.
[28] J. T. Vogelstein, J. M. Conroy, V. Lyzinski, L. J. Podrazik, S. G. Kratzer,
E. T. Harley, D. E. Fishkind, R. J. Vogelstein, and C. E. Priebe. Fast
Approximate Quadratic Programming for Graph Matching. PLoS ONE,
10(04), 2014.
[29] L. Yartseva and M. Grossglauser. On the performance of percolation
graph matching. In Proceedings of the first ACM conference on Online
social networks, pages 119–130. ACM, 2013.
[30] S. Young and E. Scheinerman. Random dot product graph models for
social networks. In Proceedings of the 5th international conference on
algorithms and models for the web-graph, pages 138–149, 2007.

24

[31] M. Zaslavskiy, F. Bach, and J.-P. Vert. A path following algorithm for
the graph matching problem. IEEE Transactions on Pattern Analysis
and Machine Intelligence, 31(12):2227–2242, 2009.
[32] Z Zhao, M Khan, V S A Kumar, and M V Marathe. Subgraph enumeration in large social contact networks using parallel color coding and
streaming. In 2010 39th International Conference on Parallel Processing, pages 594–603, 2010.

A

Auxiliary Results

We use the following McDiarmid-like concentration result in the proof of each
of our theorems.
Proposition 6 (Proposition 3.2 from [17]). Let X1 , . . . , Xm be a sequence of
independent Bernoulli random variables where E[Xi ] = pi . Let f : {0, 1}m 7→
R be such that changing any Xi to 1 − Xi changes f by at most
M = sup sup |f (X1 , . . . , Xi , . . . , Xn ) − f (X1 , . . . , 1 − Xi , . . . , Xn )|.
i

Let σ 2 = M 2
Then

X1 ,...,Xn

P

i

pi (1 − pi ) and let Y = f (X1 , . . . , Xn ).
2 /4

P[|Y − E[Y ]| ≥ tσ] ≤ 2e−t
for all 0 < t < 2σ/M .

The following proposition characterizes how a pair of correlated Bernoulli
random variables can be generated from a set of three independent random
variables. Along with a bound on the variances of these three Bernoullis this
allows for the use of Proposition 6.
Proposition 7. If X, Y are marginally Bernoulli random variables with parameter Λ and correlation ρ, then the bivariate Bernoulli distribution of
(X, Y ) is given by
(X, Y )
Y =1
Y =0
Total

X=1
X=0
Total
Λ[Λ + ρ(1 − Λ)]
(1 − Λ)Λ(1 − ρ)
Λ
Λ(1 − ρ)(1 − Λ) (1 − Λ + Λρ)(1 − Λ) 1 − Λ
Λ
1−Λ
1
25

Let Z0 , Z1 and Z2 be independent Bernoulli random variables with Z0 ∼
Bern(Λ), Z1 ∼ Bern (Λ(1 − ρ)) and Z2 ∼ Bern (Λ + ρ(1 − Λ)). Then
d

(X, Y ) ∼ (Z0 , (1 − Z0 )Z1 + Z0 Z2 ).
It also holds that
Var(Z0 ) + Var(Z1 ) + Var(Z2 ) ≤ 3Λ(1 − Λ) + ρ(1 − ρ)
Var(Z0 ) + Var(Z1 ) + Var(Z2 ) ≥ 3(1 − ρ)Λ(1 − Λ).

B

Proofs of Main Results

In this section we prove our three main theorems.
Proof of Theorem 2. Let Q ∈ Pn (with associated permutation τ ) map all
vertices in the small graph to vertices without matches in the large graph,
so that kj = |{i ≤ nc : i 6= τ (i) > nc }| = nc . Let EQ be the set of edges in A
permuted by Q, and note that EQ = [n2c ] .
eB
e − AQ
e BQ
e T ). We will construct an appropriate Λ such
Let XQ = trace(A
that P[XQ ≥ 0] ≤ exp{−C2 (log nc )2 }, which implies the desired result.
Note that for u, v ∈ [nc ], u 6= v, w, r ∈ [n], and w 6= r, we have
(
eu,v B
ew,r ) = Λu,v (Λu,v + (1 − Λu,v )Ru,v ) if {u, v} = {w, r},
E(A
Λu,v Λw,r
if {u, v} =
6 {w, r},
so that



X

EXQ = 4 

Λu,v (Λu,v + (1 − Λu,v )Ru,v ) − Λu,v Λτ (u),τ (v) 

{u,v}∈EQ


≤ 4


X

Λu,v (Λu,v + (1 − Λu,v )%) − Λu,v Λτ (u),τ (v) 

{u,v}∈EQ

Letting each Λτ (u),τ (v) > β + (1 − β)% +  be chosen to keep Λτ (u),τ (v) ∈ (0, 1)
(which is possible by the assumption that β and % are both strictly less than

26

1), we have



X

−EXQ ≥ 4 

Λu,v (β + (1 − β)% + ) − Λu,v (Λu,v + (1 − Λu,v )%)

{u,v}∈EQ

= ω( nc log nc ).
Applying Proposition 6 with M = 2, σ 2 = Θ(n2c ), we have
P [XQ ≥ 0] = P [XQ − EXQ ≥ −EXQ ]
≤ P [|XQ − EXQ | ≥ C nc log nc ]

≤ exp −C( log nc )2 ,
where C is a positive constant that can vary line-to-line.
Proof of Theore 3. Let M = {Q ∈ P : Qc = Inc }, where Qc denotes the nc ×
nc principal submatrix of Q, denote the set of permutation which correctly
e − Bk
e F = kA
e−
match the nc core vertices. Note that if P ∈ M, then kA
e T kF .
P BP
Let Q ∈
/ M (with associated permutation τ ) satisfy kj = |{i ≤ nc : i 6=
τ (i) > nc }| and kc = |{i ≤ nc : i 6= τ (i) ≤ nc |, so that Q correctly matches
the labels of nc − kc − kj of the core nc vertices across A and B. Define
eB
e − AQ
e BQ
e T ) and
XQ = 12 trace(A

 

[nc ]
EQ = {u, v} ∈
: {τ (u), τ (v)} =
6 {u, v} .
2
Note that |EQ | ≥

(nc −1)(kc +kj )
.
2

We then have


1 
T e
T
T
e
e
e
EXQ = E trace(A B − A QBQ )
2X
=
E(Au,v Bu,v − Au,v Bτ (u),τ (v) ).
{u,v}∈EQ

For u, v ∈ [nc ], u 6= v, w, r ∈ [n], and w 6= r, we have
(
2
eu,v B
ew,r ) = (1 − 2Λu,v ) + 4Λu,v (1 − Λu,v )Ru,v
E(A
(1 − 2Λu,v )(1 − 2Λw,r )
27

if {u, v} = {w, r},
if {u, v} =
6 {w, r}.

Hence,
EXQ =


X 
(1 − 2Λu,v )2(Λτ (u),τ (v) − Λu,v ) + 4Λu,v (1 − Λu,v )Ru,v .
{u,v}∈EQ

Letting γ = α(1 − α), under the assumptions of the theorem we have that
2(1 − 2Λu,v )2(Λτ (u),τ (v) − Λu,v ) + 4Λu,v (1 − Λu,v )Ru,v
≥ 2Λτ (u),τ (v) − 2Λu,v − 4Λu,v Λτ (u),τ (v) + 4Λ2u,v + 2Λu,v − 2Λ2u,v + 4γ
≥ 2Λ2τ (u),τ (v) + 2Λ2u,v − 4Λu,v Λτ (u),τ (v) + 4γ ≥ 4γ,

(3)

so that
EXQ ≥ 2(nc − 1)(kc + kj )γ
Note that XQ is a function of

  
 
kc
kj
NQ := 3 (nc − kc − kj )kc +
+2 (nc − kj )kj +
≤ 3nc (kc +kj )
2
2
independent Bernoulli random variables, and in the language of Proposition
6, we have that M = 2 and σ 2 satisfies σ 2 ≤ 3nc (kc + kj ) and
  

nc
nc − kc − kj
2
σ ≥
−
γ
2
2
Setting
2(nc − 1)(kc + kj ) · γ
σ
yields t < σ as required, and letting C be a positive constant that can vary
line by line we have that,
t=

P (XQ ≤ 0) ≤ P (|XQ − E(XQ )| ≥ E(XQ ))


C(nc − 1)(kc + kj )γ
≤ P |XP − E(XP )| ≥
·σ
σ

≤ 2 exp −C2 nc (kc + kj )γ 2 .

(4)

Define the equivalence relation “∼” on P n via P ∼ Q if P·,v = Q·,v for all
core vertices v. Note that P ∼ Q implies that XP = XQ . To prove that
there is no P ∈
/ M satisfying XP ≤ 0, it suffices to consider a single P
from each equivalence class under “∼”. Summing Eq. (4) over the (at most)
2(k +k ) 2k
nc j c nj j equivalence classes for each kc , kj and then over kc and kj yields
the desired result.
28

Proof of Theorem 4. The proof for the oracle padding scheme proceeds almost identically to the proof for the centered padding scheme. In particular,
for each Q we let XQ , EQ , kj and(kc be defined as in the proof of Theorem 3.
euv B
ewr ] = Ruv Λuv (1 − Λuv ), if {u, v} = {w, r} we
Hence, using that E[A
0,
if {u, v} =
6 {w, r},
have,
X
euv B
ewr − A
euv B
ewr ]
E[XQ ] =
E[A
{u,v}∈EQ

=

X

Ruv Λuv (1 − Λuv ) ≥ (nc − 1)(kc + kj )γ.

{u,v}∈EQ

The remainder of the proof follows mutatis mutandis.

29

