Technical Report: Graph-Structured Sparse
Optimization for Connected Subgraph Detection
Feng Chen
Computer Science Department
University at Albany – SUNY
Albany, USA
fchen5@albany.edu

arXiv:1609.09864v1 [] 30 Sep 2016

Baojian Zhou
Computer Science Department
University at Albany – SUNY
Albany, USA
bzhou6@albany.edu

Abstract—Structured sparse optimization is an important and
challenging problem for analyzing high-dimensional data in a
variety of applications such as bioinformatics, medical imaging,
social networks, and astronomy. Although a number of structured
sparsity models have been explored, such as trees, groups, clusters, and paths, connected subgraphs have been rarely explored
in the current literature. One of the main technical challenges
is that there is no structured sparsity-inducing norm that can
directly model the space of connected subgraphs, and there is
no exact implementation of a projection oracle for connected
subgraphs due to its NP-hardness. In this paper, we explore
efficient approximate projection oracles for connected subgraphs,
and propose two new efficient algorithms, namely, G RAPH -IHT
and G RAPH -GHTP, to optimize a generic nonlinear objective
function subject to connectivity constraint on the support of
the variables. Our proposed algorithms enjoy strong guarantees
analogous to several current methods for sparsity-constrained
optimization, such as Projected Gradient Descent (PGD), Approximate Model Iterative Hard Thresholding (AM-IHT), and
Gradient Hard Thresholding Pursuit (GHTP) with respect to
convergence rate and approximation accuracy. We apply our
proposed algorithms to optimize several well-known graph scan
statistics in several applications of connected subgraph detection
as a case study, and the experimental results demonstrate that
our proposed algorithms outperform state-of-the-art methods.

I. I NTRODUCTION
In recent years, structured sparse methods have attracted
much attention in many domains such as bioinformatics,
medical imaging, social networks, and astronomy [2], [4],
[14], [16], [17]. Structured sparse methods have been shown
effective to identify latent patterns in high-dimensional data
via the integration of prior knowledge about the structure
of the patterns of interest, and at the same time remain
a mathematically tractable concept. A number of structured
sparsity models have been well explored, such as the sparsity
models defined through trees [14], groups [17], clusters [16],
and paths [2]. The generic optimization problem based on a
structured sparsity model has the form
minn f (x) s.t. supp(x) ∈ M

x∈R

(1)

where f : Rn → R is a differentiable cost function, the
sparsity model M is defined as a family of structured supports: M = {S1 , S2 , · · · , SL }, where Si ⊆ [n] satisfies a
certain structure property (e.g., trees, groups, clusters), [n] =

{1, 2, · · · , n}, and the support set supp(x) refers to the set
of indexes of non-zero entries in x. For example, the popular
k-sparsity model is defined as M = {S ⊆ [n] | |S| ≤ k}.
Existing structured sparse methods fall into two main categories: 1) Sparsity-inducing norms based. The methods in
this category explore structured sparsity models (e.g., trees,
groups, clusters, and paths) [4] that can be encoded as structured sparsity-inducing norms, and reformulate Problem (1) as
a convex (or non-convex) optimization problem
minx∈Rn f (x) + λ · Ω(x)

(2)

where Ω(x) is a structured sparsity-inducing norm of M that
is typically non-smooth and non-Euclidean and λ is a trade-off
parameter. 2) Model-projection based. The methods in this
category rely on a projection oracle of M:
P(b) = arg minx∈Rn kb − xk22 s.t. supp(x) ∈ M,

(3)

and decompose the problem into two sub-problems, including
unconstrained minimization of f (x) and the projection problem P(b). Most of the methods in this category assume that
the projection problem P(b) can be solved exactly, including
the forward-backward algorithm [42], the gradient descent
algorithm [38], the gradient hard-thresholding algorithms [6],
[18], [40], the projected iterative hard thresholding [5], [7],
and the Newton greedy pursuit algorithm [41]. However,
when an exact solver of P(b) is unavailable and we have to
apply approximate projections, the theoretical guarantees of
these methods do not hold any more. We note that there is
one recent approach named as G RAPH-C OSAMP that admits
inexact projections by assuming “head” and “tail” oracles for
the projections, but is only applicable to compressive sensing
or linear regression problems [15].
We consider an underlying graph G = (V, E) defined on
the coefficients of the unknown vector x, where V = [n] and
E ⊆ V × V. We focus on the sparsity model of connected
subgraphs that is defined as
M(G, k) = {S ⊆ V | |S| ≤ k, S is connected},

(4)

where k refers to the allowed maximum subgraph size. There
are a wide array of applications that involve the search of
interesting or anomalous connected subgraphs in networks.
The connectivity constraint ensures that subgraphs reflect
changes due to localized in-network processes. We describe

a few applications below.
•

•
•

•

•
•

Detection in sensor networks, e.g., detection of traffic bottlenecks in road networks or airway networks [1]; crime
hot spots in geographic networks [22]; and pollutions in
water distribution networks [27].
Detection in digital signals and images, e.g., detection of
objects in images [15].
Disease outbreak detection, e.g., early detection of disease outbreaks from information networks incorporating
data from hospital emergency visits, ambulance dispatch
calls and pharmacy sales of over-the-counter drugs [36].
Virus detection in a computer network, e.g., detection
of viruses or worms spreading from host to host in a
computer network [24].
Detection in genome-scale interaction network, e.g., detection of significantly mutated subnetworks [21].
Detection in social media networks, e.g., detection and
forecasting of societal events [8], [9].

To the best of our knowledge, there is no existing approach
to Problem (1) for M(G, k) that is computationally tractable
and provides performance bound. First, there is no known
structured sparsity-inducing norm for M(G, k).
P The most relevant norm is fused lasso norm [39]: Ω(x) = (i,j)∈E |xi −xj |,
where xi is the i-th entry in x. This norm is able to enforce the smoothness between neighboring entries in x, but
has limited capability to recover all the possible connected
subsets as described by M(G, k) (See further discussions in
Section V Experiments). Second, there is no exact solver for
the projection oracle of M(G, k):
P(x) = arg minn kb − xk22 s.t. supp(x) ∈ M(G, k),
x∈R

(5)

as this projection problem is NP-hard due to a reduction from
classical Steiner tree problem [19]. As most existing modelprojection based methods require an exact solution to the
projection oracle P(x), these methods are inapplicable to the
problem studied here. To the best of our knowledge, there
is only one recent approach named as G RAPH-C OSAMP that
admits inexact projections for M(G, k) by assuming “head”
and “tail” oracles for the projections, but is only applicable to
compression sensing and linear regression problems [15]. The
main contributions of our study are summarized as follows:
•

•

Design of efficient approximation algorithms. Two new
algorithms, namely, G RAPH -IHT and G RAPH -GHTP,
are developed to approximately solve Problem (1) that
has a differentiable cost function and a sparsity model
of connected subgraphs M(G, k). G RAPH -GHTP is required to minimize f (x) over a projected subspace as
an intermediate step, which could be too costly in some
applications and G RAPH -IHT could be considered as a
fast variant of G RAPH -GHTP.
Theoretical guarantees and connections. The convergence rate and accuracy of our proposed algorithms are
analyzed under a smoothness condition of f (x) that is
more general than popular conditions such as Restricted
Strong Convexity/Smoothness (RSC/RSS) and Stable

Mode Restricted Hessian (SMRH). We prove that under
mild conditions our proposed G RAPH -IHT and G RAPH GHTP enjoy rigorous theoretical guarantees.
• Compressive experiments to validate the effectiveness and efficiency of the proposed techniques. Both
G RAPH -IHT and G RAPH -GHTP are applied to optimize
a variety of graph scan statistic models for the connected
subgraph detection task. Extensive experiments on a
number of benchmark datasets demonstrate that G RAPH IHT and G RAPH -GHTP perform superior to state-of-theart methods that are designed specifically for this task in
terms of subgraph quality and running time.
Reproducibility: The implementation of our algorithms and
the data sets is open-sourced via the link [11].
The remaining parts of this paper are organized as follows.
Sections II introduces the sparsity model of connected subgraphs and statement of the problem. Sections III presents
two efficient algorithms and their theoretical analysis. Section IV discusses applications of our proposed algorithms to
graph scan statistic models. Experiments on several real world
benchmark datasets are presented in Section V. Section VI
discusses related work and Section VII concludes the paper
and describes future work.
II. P ROBLEM F ORMULATION
Given an underlying graph G = (V, E) defined on the
coefficients of the unknown vector x, where V = [n],
E ⊆ V × V, and n is typically large (e.g., n > 10, 000).
The sparsity model of connected subgraphs in G is defined
in (4), and its projection oracle P(x) is defined in (5). As
this projection oracle is NP-hard to solve, we first introduce
efficient approximation algorithms for P(x) and then present
statement of the problem that will be studied in the paper.
A. Approximation algorithms for the projection oracle P(x)
There are two nearly-linear time approximation algorithms [15] for P(x) that have the following properties:
• Tail approximation (T(x)): Find a S ⊆ V such that

•

kx − xS k2 ≤ cT · 0 min
kx − xS 0 k2 ,
(6)
S ∈M(G,kT )
√
where cT = 7, kT = 5k, and xS is the restriction of
x to indices in S: we have (xS )i = xi for i ∈ S and
(xS )i = 0 otherwise.
Head approximation (H(x)): Find a S ⊆ V such that

kxS k2 ≥ cH · 0 max kxS 0 k2 ,
(7)
S ∈M(G,kH )
p
where cH = 1/14 and kH = 2k.
It can be readily proved that, if cT = cH = 1, then T(x) =
H(x) = P(x), which indicates that these two approximations
(T(x) and H(x)) stem from the fact that cT > 1 and cH < 1.
B. Problem statement
Given a predefined cost function f (x) that is differentiable,
the input graph G, and the sparsity model of connected

Fig. 1. Illustration of G RAPH -IHT and G RAPH -GHTP on the main steps of each iteration. In this example, the gray scale of each node i encodes the weight
of this node wi ∈ R, f (x) = −wT x + 12 kxk2 , and the maximum size of subgraphs is set to k = 6, where w = [w1 , · · · , wn ]T . The resulting problem
tends to find a connected subgraph with the largest overall weight (See discussion about the effect of 12 kxk2 in Section IV). In each iteration, S i is the
connected subset of nodes and its induced subgraph is denoted as GS i . The sequence of intermediate vectors and subgraphs are: (x0 , GS 1 ), · · · , (xi , GS i ).

subgraph M(G, k), the problem to be studied is formulated
as:
min f (x) s.t. supp(x) ∈ M(G, k).

x∈Rn

(8)

Problem (8) is difficult to solve as it involves decision variables
from a nonconvex set that is composed of many disjoint
subsets. In this paper, we will develop nearly-linear time
algorithms to approximately solve Problem (8). The key
idea is decompose this problem to sub-problems that are
easier to solve. These sub-problems include an optimization
sub-problem of f (x) that is independent on M(G, k) and
projection approximations for M(G, k), including T(x) and
H(x). We will design efficient algorithms to couple these subproblems to obtain global solutions to Problem (8) with good
trade-off on running time and accuracy.
III. A LGORITHMS
This section first presents two efficient algorithms, namely,
G RAPH -IHT and G RAPH -GHTP, and then analyzes their time
complexities and performance bounds.
A. Algorithm G RAPH -IHT
The proposed G RAPH -IHT algorithm generalizes the traditional algorithm named as projected gradient descent [5], [7]
that requires a exact solver of the projection oracle P(x). The
high-level summary of G RAPH -IHT is shown in Algorithm 1
and illustrated in Figure 1 (b). The procedure generates a
sequence of intermediate vectors x0 , x1 , · · · from an initial
approximation x0 . At the i-th iteration, the first step (Line
5) first calculates the gradient “∇f (xi )”, and then identifies a
subset of nodes via head approximation that returns a support
set with the head value at least a constant factor of the
optimal head value: “Ω ← H(∇f (xi ))”. The support set Ω
can be interpreted as the subspace where the nonconvex set
“{x | supp(x) ∈ M(G, k)}” is located, and the projected gradient in this subspace is: “∇Ω f (xi )”. The second step (Line
6) calculates the projected gradient descent at the point xi with
step-size η: “b ← xi − η · ∇Ω f (xi )”. The third step (Line 7)

identifies a subset of nodes via tail approximation that returns a
support set with tail value at most a constant times larger than
the optimal tail value: “S i+1 ← T(b)”. The last step (Line 8)
calculates the intermediate solution xi+1 : xi+1 = bS i+1 . The
previous two steps can be interpreted as the projection of b to
the nonconvex set “{x | supp(x) ∈ M(G, k)}” using the tail
approximation.
Algorithm 1 G RAPH -IHT
1: Input: Input graph G, maximum subgraph size k, and step
size η (1 by default).
2: Output: The estimated vector x̂ and the corresponding
connected subgraph S.
3: i ← 0, xi ← 0; S i ← ∅;
4: repeat
5:
Ω ← H(∇f (xi ));
6:
b ← xi − η · ∇Ω f (xi );
7:
S i+1 ← T(b);
8:
xi+1 ← bS i+1 ;
9:
i ← i + 1;
10: until halting condition holds
11: return x̂ = xi and S = GS i ;

B. Algorithm G RAPH -GHTP
The proposed G RAPH -GHTP algorithm generalizes the traditional algorithm named as Gradient Hard Threshold Pursuit
(GHTP) that is designed specifically for the k-sparsity model:
M = {S ⊆ [n] | |S| ≤ k} [40]. The high-level summary
of G RAPH -GHTP is shown in Algorithm 2 and illustrated in
Figure 1 (c). The first two steps (Line 5 and Line 6) in each
iteration is the same as the first two steps (Line 5 and Line
6) of G RAPH -IHT, except that we return the support of the
projected gradient descent: “Ψ ← supp(xi − η · ∇Ω f (xi ))”,
in which pursuing the minimization will be most effective.
Over the support set S, the function f is minimized to
produce an intermediate estimate at the third step (Line 7):

Algorithm 2 G RAPH -GHTP
1: Input: Input graph G, maximum subgraph size k, and step
size η (1 by default).
2: Output: The estimated vector x̂ and the corresponding
connected subgraph S.
3: i ← 0, xi ← 0; S i ← ∅;
4: repeat
5:
Ω ← H(∇f (xi ))
6:
Ψ ← supp(xi − η · ∇Ω f (xi ));
7:
b ← arg minx∈Rn f (x) s.t. supp(x) ⊆ Ψ;
8:
S i+1 ← T(b);
9:
xi+1 ← bS i+1 ;
10:
i ← i + 1;
11: until halting condition holds
12: return x̂ = xi and S = GS i ;

“b ← arg minx∈Rn f (x) s.t. supp(x) ⊆ Ω”. The fourth and
fifth steps (Line 8 and Line 9) are the same as the last two
steps (Line 7 and Line 8) of G RAPH -IHT in each iteration.
C. Relations between G RAPH -IHT and G RAPH -GHTP
These two algorithms are both variants of gradient descent.
In overall, G RAPH -GHTP converges faster than G RAPH -IHT
as it identifies a better intermediate solution in each iteration
by minimizing f (x) over a projected subspace {x | supp(x) ⊆
Ω}. If the cost function f (x) is linear or has some special
structure, this intermediate step can be conducted in nearlylinear time. However, when this step is too costly in some
applications, G RAPH -IHT is preferred.
D. Theoretical Analysis of G RAPH -IHT
In order to demonstrate the accuracy of estimates using
Algorithm 1, we require that the cost function f (x) satisfies
the Weak Restricted Strong Convexity (WRSC) condition as
follows:
Definition III.1 (Weak Restricted Strong Convexity Property
(WRSC)). A function f (x) has the (ξ, δ, M)-model-WRSC
if ∀x, y ∈ Rn and ∀S ∈ M with supp(x) ∪ supp(y) ⊆ S, the
following inequality holds for some ξ > 0 and 0 < δ < 1:
kx − y − ξ∇S f (x) + ξ∇S f (y)k2 ≤ δkx − yk2 .

(9)

The WRSC is weaker than the popular Restricted Strong
Convexity/Smoothness (RSC/RSS) conditions that are used in
theoretical analysis of convex optimization algorithms [40].
The RSC condition basically characterizes cost functions
that have quadratic bounds on the derivative of the objective function when restricted to model-sparse vectors. The
RSC/RSS conditions imply condition WRSC, which indicates that WRSC is no stronger than RSC/RSS [40]. In the
special case where f (x) = ky−Axk22 and ξ = 1, the condition
(ξ, δ, M)-model-WRSC reduces to the well known Restricted
Isometry Property (RIP) condition in compressive sensing.
Theorem III.1. Consider the sparsity model of connected
subgraphs M(G, k) for some k ∈ N and a cost function
f : Rn → R that satisfies the (ξ, δ, M(G, 5k))-model-WRSC

condition. If η = cH (1 − δ) − δ then for any x ∈ Rn such
that supp(x) ∈ M(G, k), with η > 0 the iterates of Algorithm
2 obey
kxi+1 − xk2 ≤ αkxi − xk2 + βk∇I f (x)k2

(10)

where
α0 = cH (1 − δ) − δ, β0 = δ(1 + cH ),

q

2(1 + cT )
η
η
2
1 − α0 + (2 − )δ + 1 −
,
α=
1−δ
ξ
ξ
!
√
√
√
√
1 + cT
2β0
2α0 β0
β=
,
(1 + 2 2)ξ + (2 − 2 2)η +
+p
1−δ
α0
1 − α02
√

and I = arg maxS∈M(G,8k) k∇S f (x)k2
Before we prove this result, we give the following two
lemmas III.2 and III.3.
Lemma III.2. [40] Assume that f is a differentiable function.
If f satisfies condition (ξ, δ, M)-WRSC, then ∀x, y ∈ Rn with
supp(x) ∪ supp(y) ⊂ S ∈ M, the following two inequalities
hold
1−δ
1+δ
kx − yk2 ≤ k∇S f (x) − ∇S f (y)k2 ≤
kx − yk2
ξ
ξ
1+δ
kx − yk22
f (x) ≤ f (y) + h∇f (y), x − yi +
2ξ

Lemma III.3. Let α0 = cH (1 − δ) − δ, β0 = ξ(1 + cH ),
ri = xi − x, and Ω = H(∇f (xi )). Then
kriΓc k2

≤

q

#
β0
α0 β0
k∇I f (x)k2
+
+p
α0
1 − α02
"

1−

α02 kri k2

where I = arg maxS∈M(G,8k) k∇S f (x)k2 . We assume that cH
and δ are such that α0 > 0.
Proof: Denote Φ = supp(x) ∈ M(G, k), Ω =
H(∇f (xi )) ∈ M(G, 2k), ri = xi − x, and Λ = supp(ri ) ∈
M(G, 6k). The component k∇Γ f (xi )k2 can be lower bounded
as
k∇Γ f (xi )k2

≥

cH k∇Φ f (xi )k2

≥

cH (k∇Φ f (xi ) − ∇Φ f (x)k2 − k∇Φ f (x)k2
cH (1 − δ) i
kr k2 − cH k∇I f (x)k2 ,
ξ

≥

where the first inequality follows from the definition of
head approximation and the last inequality follows from
Lemma III.2 of our paper. The component k∇Γ f (xi )k2 can
also be upper bounded as
k∇Γ f (xi )k2

≤
≤

≤

≤

1
kξ∇Γ f (xi ) − ξ∇Γ f (x)k2 + k∇Γ f (x)k2
ξ
1
kξ∇Γ f (xi ) − ξ∇Γ f (x) − riΓ + riΓ k2 +
ξ
k∇Γ f (x)k2
1
kξ∇Γ∪Ω f (xi ) − ξ∇Γ∪Ω f (x) − riΓ∪Ω k2 +
ξ
kriΓ k2 + k∇Γ f (x)k2
δ
1
· kri k2 + kriΓ k2 + k∇I f (x)k2 ,
ξ
ξ

where the fourth inequality follows from condition
(ξ, δ, M(G, 8k))-WRSC and the fact that riΓ∪Ω = ri .
Combining the two bounds and grouping terms, we obtain

the inequality:
kriΓ k ≥ α0 kri k2 − ξ(1 + cH )k∇I f (x)k2 .
have kriΓ k ≥ α0 kri k2 − β0 k∇I f (x)k2 . After a

We
number of
algebraic manipulations, we obtain the inequality
kriΓc k2

"
#
q
β0
α0 β0
i
2
k∇I f (x)k2 ,
+p
≤ 1 − α0 kr k2 +
α0
1 − α02

which proves the lemma.
We give the formal proof of III.1.
Proof: From the traingle inequality, we have
kri+1 k2

= kxi+1 − xk2
= kbΨ − xk2
≤

kb − xk2 + kb − bΨ k2

≤ (1 + cT )kb − xk2
=
=

E. Theoretical Analysis of G RAPH -GHTP

i

i

i

i

i

(1 + cT )kx − η∇Ω f (x ) − x k2
(1 + cT )kr − η∇Ω f (x )k2 ,

i

where ∇Ω f (x ) is the projeceted vector of f (xi ) in which
the entries outoside Ω are set to zero and the entries in Ω are
unchanged. kri − η∇Ω f (xi )k2 has the inequalities
kri − η∇Ω f (xi )k2 = kriΩc + riΩ − η∇Ω f (xi )k
≤ kriΩc k2 + kriΩ − η∇Ω f (xi ) + η∇Ω f (x) − η∇Ω f (x)k
≤ kriΩc k2 + kriΩ − η∇Ω f (xi ) + η∇Ω f (x)k + kη∇Ω f (x)k
≤ kriΩc k2 + kriΩ − ξ∇Ω f (xi ) + ξ∇Ω f (x)k+
(ξ − η)k∇Ω f (xi ) − ∇Ω f (x)k2 k + kη∇Ω f (x)k2
≤ kriΩc k2 + (1 − η/ξ + (2 − η/ξ)δ)kri k2 + ηk∇I f (x)k2
where the last inequality follows from condition (ξ, δ, M)WRSC and Lemma III.2. From Lemma III.3, we have
q
hβ
α0 β0 i
0
+p
k∇I f (x)k2
kriΓc k2 ≤ 1 − α02 kri k2 +
α0
1 − α02
Combining the above inequalities, we prove the theorem.
Our proposed G RAPH -IHT generalizes several existing
sparsity-constrained optimization algorithms: 1) Projected
Gradient Descent (PGD) [30]. If we redefine H(b) =
supp(b) and T(b) = supp(P(b)), where P(b) is the projection
oracle defined in Equation (5), then G RAPH -IHT reduces
to the PGD method; 2) Approximated Model-IHT(AMIHT) [13]. If the cost function f (x) is defined as the least
square cost function f (x) = ky − Axk22 , then ∇f (x) has
the specific form −AT (y − A) and G RAPH -IHT reduces to
the AM-IHT algorithm, the state-of-the-art variant of IHT for
compressive sensing and linear regression problems. In particular, let e = y − Ax. The component
k∇f (xi )k2 = kAT ek2
√
is upper bound by bounded by 1 + δkek2 [13], Assume that
ξ = 1 and η = 1. Condition (ξ, η, M)-WRSC then reduces
to the RIP condition in compressive sensing. The convergence
inequality (10) then reduces to
kx

i+1

0

i

0

− xk2 ≤ α kx − xk2 + β kek2 ,
(11)
h
i
p
2
where α = (1 + cT ) δ + 1 − α0 and
√
h (α + β )√1 + δ
α0 β0 ( 1 + δ) i
0
0
0
p
β = (1 + cT )
+
.
α0
1 − α02
0

Surprisingly, the above convergence inequality is identical
to the convergence inequality of AM-IHT derived in [13]
based on the RIP condition, which indicates that G RAPH -IHT
has the same convergence rate and approximation error as
AM-IHT, although we did not make any attempt to explore
the special properties of the RIP condition. We note that the
convergence properties of G RAPH -IHT hold in fairly general
setups beyond compressive sensing and linear regression. As
we consider G RAPH -IHT as a fast variant of G RAPH -GHTP,
due to space limit we ignore the discussions about the convergence condition of G RAPH -IHT. The theoretical analysis
of G RAPH -GHTP to be discussed in the next subsection can
be readily adapted to the theoretical analysis of G RAPH -IHT.

Theorem III.4. Consider the sparsity model of connected
subgraphs M(G, k) for some k ∈ N and a cost function
f : Rn → R that satisfies the (ξ, δ, M(G, 5k))-model-WRSC
condition. If η = cH (1 − δ) − δ then for any x ∈ Rn such
that supp(x) ∈ M(G, k), with η > 0 the iterates of Algorithm
2 obey
kxi+1 − xk2 ≤ αkxi − xk2 + βk∇I f (x)k2

(12)

where
α0 = cH (1 − δ) − δ, β0 = δ(1 + cH ),

q

2(1 + cT )
η
η
1 − α02 + (2 − )δ + 1 −
,
α=
1−δ
ξ
ξ
!
√
√
√
√
1 + cT
2β0
2α0 β0
(1 + 2 2)ξ + (2 − 2 2)η +
+p
β=
,
1−δ
α0
1 − α02
√

and I = arg maxS∈M(G,8k) k∇S f (x)k2
Proof: Denote Ω = H(∇f (xi )) and Ψ = supp(xi − η ·
∇Ω f (xi )). Let ri+1 = xi+1 − x. kri+1 k2 is bounded as
kri+1 k2 = kxi+1 − xk2

≤
≤
≤

kxi+1 − bk2 + kx − bk2
cT kx − bk2 + kx − bk2
(1 + cT )kx − bk2 ,
(13)

where the second inequality follows from the definition of tail
approximation. The component k(x − b)Ψ k22 is bounded as
k(x − b)Ψ k22 = hb − x, (b − x)Ψ i
= hb − x − ξ∇Ψ f (b) + ξ∇Ψ f (x), (b − x)Ψ i −
hξ∇Ψ f (x), (b − x)Ψ i
≤ δkb − xk2 k(b − x)Ψ k2 + ξk∇Ψ f (x)k2 k(b − x)Ψ k2 ,

where the second equality follows from the fact that
∇S f (b) = 0 since b is the solution to the problem in the
third Step (Line 7) of G RAPH -GHTP, and the last inequality
can be derived from condition (ξ, δ, M(G, 8k))-WRSC. After
simplification, we have
k(x − b)Ψ k2 ≤ δkb − xk2 + ξk∇Ψ f (x)k2 .

It follows that

kx − bk2 ≤ k(x − b)Ψ k2 + k(x − b)Ψc k2
≤ δkb − xk2 + ξk∇Ψ f (x)k2 + k(x − b)Ψc k2 .

After rearrangement we obtain
kb − xk2

≤

k(b − x)Ψc k2
ξk∇Ψ f (x)k2
+
,
1−δ
1−δ

(14)

where this equality follows from the fact that supp(b) ⊆ S.
Let Φ = supp(x) ∈ M(G, k).
k(xi − η∇Ω f (xi ))Φ k2 ≤ k(xi − η∇Ω f (xi ))Ψ k2 ,

as Ψ = supp(xi −η·∇Ω f (xi )). By eliminating the contribution
on Φ ∩ Ψ, we derive
k(xi − η∇Ω f (xi ))Φ\Ψ k2 ≤ k(xi − η∇Ω f (xi ))Ψ\Φ k2

For the right-hand side, we have

k(xi − η∇Ω f (xi ))Ψ\Φ k2 ≤

k(xi − x − η∇Ω f (xi ) + η∇Ω f (x))Ψ\Φ k2 + ηk∇Ω∪Ψ f (x)k2 ,

where the inequality falls from the fact that Φ = supp(x).
From the left-hand
side,i we have
i
k(x − η∇Ω f (x ))Φ\Ψ k2 ≤ −ηk∇Ω∪Φ f (x)k2 +

k(xi − x − η∇Ω f (xi ) + η∇Ω f (x))Φ\Ψ + (x − b)Ψc k2

where the inequality follows from the fact that bΨc = 0,
xΦ\Ψ = xΨc , and −xΦ\Ψ + (x − b)Ψc = 0. Let Φ∆Ψ be
the symmetric difference of the set Φ and Ψ. It follows that
≤

k(b − x)Ψc k2
√
2k(xi − x − η∇Ω f (xi ) + η∇Ω f (x))Φ∆Ψ k2 + 2ηk∇I f (x)k2
√
≤ 2k(xi − x − ξ∇Ω f (xi ) + ξ∇Ω f (x))Φ∆Ψ k2 +
√
2(ξ − η)k(∇Ω f (xi ) + ∇Ω f (x))Φ∆Ψ k + 2ηk∇I f (x)k2
√
≤ 2k(riΩc + riΩ − ξ∇Ω f (xi ) + ξ∇Ω f (x))Φ∆Ψ k2 +
√
2(ξ − η)k(∇Ω f (xi ) − ∇Ω f (x))Φ∆Ψ k + 2ηk∇I f (x)k2
√
√
≤ 2kriΩc k + 2k(riΩ − ξ∇Ω f (xi ) + ξ∇Ω f (x))Ψ∆Ψ k2 +
√
2(ξ − η)k(∇Ω f (xi ) − ∇Ω f (x))Ψ∆Ψ k + 2ηk∇I f (x)k2
√
√
i
≤ 2krΩc k + 2kri − ξ∇Ω∪Ψ∪Φ f (xi ) + ξ∇Ω∪Ψ∪Φ f (x)k2 +

√
2(ξ − η)k(∇Ω∪Ψ∪Φ f (xi ) − ∇Ω∪Ψ∪Φ f (x))Ψ∆Ψ k + 2ηk∇I f (x)k2



√
√
η
η
i
c k2 +
≤ 2krΩ
δ+1−
kri k +
2
2−
ξ
ξ
√
√ 
2
2ξ + (1 − 2)η k∇I f (x)k2 ,

where the first inequality follows from the fact that
ηk∇Ω∪Φ f (x)k2 + ηk∇Ψ∪Φ∪Ω f (x)k2 ≤ 2ηk∇I f (x)k2 ,
the third inequality follows as xi − x = ri = riΩc + riΩ , the
fourth inequality follows from the fact that k(riΩc )Φ∆Ψ k2 ≤
kriΩc k2 , the fifth inequality follows as ri ⊆ Ω ∪ Ψ ∪ Φ, and the
last inequality follows from condition (ξ, δ, M(G, 8k))-WRSC
and Lemma III.2. From Lemma
III.3, we have
"
#
kriΩc k2 ≤

p
ξ(1 + cH )
ξη(1 + cH )
k∇I f (x)k2
1 − η 2 kri k2 +
+ p
η
1 − η2

Combining (14) and above inequalities, we prove the theorem.
Theorem III.4 shows the estimator error of G RAPH -GHTP
is determined by the multiple of k∇S f (x)k2 , and the convergence rate is geometric. Specifically, if x is an uncontrained
minimizer of f (x), then ∇f (x) = 0. It means G RAPH -GHTP
is guaranteed to obtain the true x to arbitrary precision. The
estimation error is negligible when x is sufficiently close to
an unconstrained minimizer of f (x) as k∇S f (x)k2 is a small
value. The parameter
√

q

2(1 + cT )
η
η
1 − α02 + (2 − )δ + 1 −
< 1,
α=
1−δ
ξ
ξ
controls the convergence rate of G RAPH -GHTP. Our algorithm allows an exact recovery if α < 1. As δ is an arbitrary

constant parameter, it can be an arbitrary small positive value.
Let η be ξ and δ be an arbitrary small positive value, the
parameters cH and cT satisfy the following inequality
c2H > 1 − 1/(1 + cT )2 .

(15)

It is noted that the head and tail approximation algorithms
described in [15] do not meet the inequality (15). Nonetheless,
the approximation factor cH of any given head approximation
algorithm can be boosted to any arbitrary constant c0H < 1,
which leads to the satisfaction of the above condition as shown
in [15]. Boosting the head-approximation algorithm, though
strongly suggested by [13], is not empirically necessary.
Our proposed G RAPH -GHTP has strong connections to the
recently proposed algorithm named as Gradient Hard Thresholding Pursuit (GHTP) [40] that is designed specifically for
the k-sparsity model: M = {S ⊆ [n] | |S| ≤ k}. In particular,
if we redefine H(b) = supp(b) and T(b) = supp(P(b)),
where P(b) is the projection oracle defined in Equation (5),
and assume that there is an algorithm that solves the projection
oracle exactly, in which the sparsity model does not require
to be the k-sparsity model. It then follows that the upper
bound of kriΩc k2 stated in Lemma III.2 in Appendix is updated
as kriΩc k2 ≤ 0, since supp(ri ) = Ω and riΩc = 0. In
addition, the multiplier (1 + cT ) is replaced as 1 as the first
inequality (13) in the proof of Theorem III.4 in Appendix is
updated as kri+1 k2 ≤ kx−bk2 , instead of the original version
kri+1 k2 ≤ (1 + cT )kx − bk2 . After these two changes, the
shrinkage rate α is updated as
√ 

2
η
η
(2 − )δ + 1 −
,
(16)
α=
1−δ
ξ
ξ
which is the same as the shrinkage rate of G RAPH -GHTP
as derived in [40] specifically for the k-sparsity model. The
above shrinkage rate α (16) should satisfy the condition α < 1
to ensure the geometric convergence of G RAPH -GHTP, which
implies that
√
√
√
√
(17)
η > ((2 2 + 1)δ + 2 − 1)ξ/( 2 + 2δ).
√
It follows that if δ < 1/( 2 + 1), a step-size η < ξ can
always be found to satisfy the above inequality. This constant
condition of δ is analogous to the constant condition of stateof-the-art compressive sensing methods that consider noisy
measurements [23] under the assumption of the RIP condition.
We derive the analogous constant using the WRSC condition
that weaker than the RIP condition.
As discussed above, our proposed G RAPH -GHTP has
connections to GHTP on the shrinkage rate of geometric
convergence. We note that the shrinkage rate of our proposed
G RAPH -GHTP stated in Theorem III.4 is derived based on
head and tail approximations of the sparsity model of connected subgraphs M(G, k), instead of the k-sparsity model
that has an exact projection oracle solver. Our convergence
properties hold in fairly general setups beyond k-sparsity
model, as a number of popular structured sparsity models such
as the “standard” k-sparsity, block sparsity, cluster sparsity,
and tree sparsity can be encoded as special cases of M(G, k).

Theorem III.5. Let x ∈ Rn such that supp(x) ∈ M(G, k),
and f : Rn → R be cost function that satisfies condition
(ξ, δ, M(8k, g))-WRSC. Assuming that α < 1, G RAPH -GHTP
(or G RAPH -IHT) returns a x̂ such that, supp(x̂) ∈ M(5k, g)
β
) is a fixed
and kx − x̂k2 ≤ ck∇I f (x)k2 , where c = (1 + 1−α
constant. Moreover, G RAPH -GHTP runs in time

O (T + |E| log3 n) log(kxk2 /k∇I f (x)k2 ) ,
(18)
where T is the time complexity of one execution of the
subproblem in Step 6 in G RAPH -GHTP (or Step 5 in G RAPH IHT). In particular, if T scales linearly with n, then G RAPH GHTP (or G RAPH -IHT) scales nearly linearly with n.
Proof: The i-th iterate of G RAPH -GHTP (or G RAPH IHT) satisfies
β
k∇I f (x)k2 .
1 − α m

kx − xi k2 ≤ αi kxk2 +

(19)

l 
1
2
iterations, G RAPH After t = log k∇Ikxk
f (x)k2 / log α
GHTP (or G RAPH -IHT) returns an estimate x̂ satisfying kx−
β
x̂k2 ≤ (1 + 1−α
)k∇I f (x)k2 . The time complexities of both
head approximation and tail approximation are O(|E| log3 n).
The time complexity of one iteration in G RAPH -GHTP (or
3
G RAPH -IHT)l is 
(T + |E| log
 n), and
m the total number of
kxk2
iterations is log k∇I f (x)k2 / log α1 , and the overall time
follows. G RAPH -GHTP and G RAPH -IHT are only different
in the definition of α and β in this Theorem.
As shown in Theorem 18, the time complexity of G RAPH GHTP is dependent on the total number of iterations and the
time cost (T ) to solve the subproblem in Step 6. In comparison,
the time complexity of G RAPH -IHT is dependent on the total
number of iterations and the time cost T to calculate the
gradient ∇f (xi ) in Step 5. It implies that, although G RAPH GHTP converges faster than G RAPH -IHT, the time cost to
solve the subproblem in Step 6 is often much higher than the
time cost to calculate a gradient ∇f (xi ), and hence G RAPH IHT runs faster than G RAPH -GHTP in practice.
IV. A PPLICATIONS ON G RAPH S CAN S TATISTICS
In this section, we specialize G RAPH-IHT and G RAPHGHTP to optimize a number of well-known graph scan
statistics for the task of connected subgraph detection, including elevated mean scan (EMS) statistic [29], Kulldorff’s
scan statistic [25], and expectation-based Poisson (EBP) scan
statistic [26]. Each graph scan statistic is defined as the
generalized likelihood ratio test (GLRT) statistic of a specific
hypothesis testing about the distributions of features of normal
and abnormal nodes. The EMS statistic corresponds to the
following GLRT test: Given a graph G = (V, E), where
V = [n] and E ⊆ V × V, each node i is associated with a
random variable xi :
xi = µ · 1(i ∈ S) + i , i ∈ V,

(20)

where |µ| represents the signal strength and i ∈ N (0, 1). S
is some unknown anomalous cluster that forms as a connected
subgraph. The task is to decide between the null hypothesis
(H0 ): ci ∈ N (0, 1), ∀i ∈ V and the alternative (H1 (S)): ci ∈
N (µ, 1), ∀i ∈ S and ci ∈ N (0, 1), ∀i ∈
/ S. The EMS statistic

is defined as the GLRT function under this hypothesis testing:
Prob(Data|H1 (S))
1 X
F (S) =
ci .
(21)
=p
Prob(Data|H0 )
|S|
i∈S

The problem of connected subgraph detection based on the
EMS statistic is then formulated as
1 X 2
ci ) s.t. S ∈ M(G, k),
(22)
min − (
S⊆V
|S|
i∈S

where the square of the EMS scan statistic is considered to
make the function smooth, and this transformation does not
infect the optimum solution. Let the {0, 1}-vectors form of S
be x ∈ {0, 1}n , such that supp(x) = S. Problem (22) can be
reformulated as
min

x∈{0,1}n

−(cT x)2 /(1T x) s.t. supp(x) ∈ M(G, k),

(23)

where c = [c1 , · · · , cn ]T . To apply our proposed algorithms,
we relax the input domain of x and maximize the strongly
convex function [3]:
1
min −(cT x)2 /(1T x) + xT x s.t. supp(x) ∈ M(G, k). (24)
x∈Rn
2
The connected subset of nodes can be found as the subset of
indexes of positive entries in x̂, where x̂ refers to the solution
of the Problem (24). Assume that c is normalized and ci ≤ 1,
∀i. Let ĉ = max{c1 , · · · , cn }. The Hessian matrix of the above
objective function satisfies the following conditions
cT x
cT x
(1 − ĉ2 ) · I  I − (c − T 1)(c − T 1)T  1 · I. (25)
1 x
1 x
According to Lemma 1 (b) in [40]), the objective function
f (x) satisfies condition (ξ, δ, M(G, 8k))-WRSC that
p
δ = 1 − 2ξ(1 − ĉ2 ) + ξ 2 ,
for any ξ such that ξ < 2(1 − ĉ2 ). The geometric convergence
of G RAPH -GTHP as shown in Theorem III.4 is guaranteed.
Different from the EMS statistic that is defined for numerical features based on Gaussian distribution, the Kulldorff’s
scan statistic and Expectation Based Poisson statistic (EBP)
are defined for count features based on Poisson distribution. In
particular, each node i is associated with a feature ci , the count
of events (e.g., crimes, flu infections) observed at the current
time, and a feature bi , the expected count (or ‘baseline’) of
events by using historical data. Let c = [c1 , · · · , cn ]T and
b = [b1 , · · · , bn ]T . The Kulldorff’s scan statistic and EBP
scan statistics are described Table I. We note that these two
scan statistics do not satisfy the WRSC condition, but as
demonstrated in our experiments, our proposed algorithms
perform empirically well for all the three scan statistics, and
in particular, our proposed G RAPH -GHTP converged in less
than 10 iterations in all the settings.
V. EXPERIMENTS
This section evaluates the performance of our proposed
methods using four public benchmark data sets for connected
subgraph detection. The experimental code and data sets are
available from the Link [11] for reproducibility.

TABLE I
The three typical graph scan statistics that are tested in our experiments (The vectors x, c, and b are defined in Section IV)
Score Functions

Definition

Kulldorff’s Scan Statistic [25]

cT x log

Expectation-based Poisson Statistic (EBP) [26]
Elevated Mean Scan
Statsitic (EMS) [29]

(1T c

−

cT x log

Applications

T
cT x
c
− 1T c log 11T b
bT x
T
T
1 c−c x
T
c x) log 1T b−bT x

cT x
bT x

+

The statistic is used for anomalous pattern detection in graphs with count
features, such as detection of traffic bottlenecks in sensor networks [1], [22],
detection of anomalous regions in digitals and images [10], detection of attacks
in computer networks [24], disease outbreak detection [36], and various others.

+ bT x − cT x

This statistic is used for the same applications as above [1], [12], [25], [36],
[37], but has different assumptions on data distribution [25].
cT x/1T x
This statistic is used for anomalous pattern detection in graphs with numerical
features, such as event detection in social networks, network surveillance,
disease outbreak detection, biomedical imaging [29], [35]
TABLE II
Summary of dataset settings in the experiments. (For the network of each dataset, we only use its maximal connected component if it is not fully connected).
Dataset

Application

Training & Testing Time Periods

BWSN

Detection of
contaminated nodes

CitHepPh

Detection of emerging
research areas
Detection of most
congested subgraphs
Detection of crime
hot spots

Training: Hours 3 to 5 with 0% noise
Testing: Hours 3 to 5 with
2%, 4%, · · · , 10% noises
Testing: 1999 to 2002

Traffic
ChicagoCrime

Observed value at
node v: cv
Sensor value (0 or 1)

Count of citations

Testing: Mar.2014, 5AM to 10PM

− log(pt (v)/µ) 1

Testing: Year of 2015

Count of burglaries

Baseline value at node v: bv

# of Nodes

# of Edges

# of snapshots

Average sensor value for EBP
Constant ‘1’ for Kulldorff

12,527

14,323

hourly: 3 × 6

Average count of citations for EBP
Maximum count of citations for Kulldorff
None
(EBP and Kulldorff are not applicable)
Average count of burglaries for EBP
Maximum count of burglaries for Kulldorff

11,895

75,873

yearly: 1 × 11

1,723

5,301

per 15 min: 68 × 304

46,357

168,020

yearly: 1 × 15

pt (v) refers to the statistical p-value of node v at time t that is calculated via empirical calibration based on historical speed values of v from 2013 June.
1st to 2014 Feb. 29; and µ is a significance level threshold and is set to 0.15. The larger this value − log(pt (v)/µ), the more congested in the region near v.

1

A. Experiment Design
Datasets: 1) BWSN Dataset. A real-world water network is offered in the Battle of the Water Sensor Networks
(BWSN) [28]. That has 12,527 nodes and 14,323 edges.
In order to simulate a contaminant sub-area, 4 nodes with
chemical contaminant plumes, which were distributed in this
sub-area, were generated. We use the water network simulator
EPANET [31] that was employed in BWSN for a period of 3
hours to simulate the spreads for contaminant plumes on this
graph. If a node is polluted by the chemical, then its sensor
reports 1, otherwise, 0, in each hour. To test the tolerance of
noise of our methods, K ∈ {2, 4, 6, 8, 10} percent vertices
were selected randomly, and their sensor reports were set to
0 if their original reports were 1 and vice versa. Each hour
has a graph snapshot. The snapshots corresponding to the 3
hours that have 0% noise are considered for training, and the
snapshots that have 2%, · · · , 10% noise reports for testing. The
goal is to detect a connected subgraph that is corresponding to
the contaminant sub-area. 2) CitHepPh Dataset. We downloaded the high energy physics phenomenology citation data
(CitHepPh) from Stanford Network Analysis Project (SNAP)
[20]. This citation graph contains 11,897 papers corresponding
to graph vertices and 75,873 edges. An undirected edge
between two vertices (papers) exists , if one paper is cited by
another. The period of these papers published is from January
1992 to April 2002. Each vertex has two attributes for each
specific year (t = 1992, · · · , t = 2002). We denote the number
of citations of each specific year as the first attribute and
the average citations of all papers in that year as the second
attribute. The goal is to detect a connected subgraph where the
number of citations of vertices (papers) in this subgraph are abnormally high compared with the citations of vertices that are
not in this subgraph. This connected subgraph is considered
as a potential emerging research area. Since the training data

is required for some baseline methods, the data before 1999 is
considered as the training data, and the rest from 1999 to 2002
as the testing data. 3) Traffic Dataset. Road traffic speed data
from June 1, 2013 to Mar. 31, 2014 in the arterial road network
of the Washington D.C. region is collected from the INRIX
database (http://inrix.com/publicsector.asp), with 1,723 nodes
and 5,301 edges. The database provides traffic speed for each
link at a 15-minute rate. For each 15-minute interval, each
method identities a connected subgraph as the most congested
region. 4) ChicagoCrime Dataset. We collected crime data
from City of Chicago [https://data.cityofchicago.org/PublicSafety/Crimes-2001-to-present/ijzp-q8t2] from Jan. 2001 to
Dec. 2015. There are 46,357 nodes (census blocks) and
168,020 edges (Two census blocks are connected with each
other if they are neighbours). Specifically, we collected all
records of burglaries from 2001 to 2015. The data covers
burglaries in the period from Jan. 2001 to Dec. 2015. Each
vertex has an attribute denoting the number of burglaries in
sepcific year and average number of burglaries over 10 years.
We aim to detect connected census areas which has anomaly
high burglaries accidents. The data before 2010 is considered
as training data, and the data from 2011 to 2015 is considered
as testing data.
Graph Scan Statistics: As shown in Table I, three graph
scan statistics were considered as the scoring functions of
connected subgraphs, including Kulldorff’s scan statistic [25],
expectation-based Poisson (EBP) scan statistic [26], and elevated mean scan (EMS) statistic [29]. The first two statistic
functions require that each vertex v has a count cv representing
the count of events observed at that vertex, and an expected
count (‘baseline‘) bv . For EMS statistic, only cv is used. We
need to normalize cv for EMS as it is defined based on the
assumptions of standard normal distribution for normal values
and shifted-mean normal distribution for abnormal values.

Table II provides details about the calculations of cv and bv
for each data set.
Comparison Methods: We compared our proposed
methods with four state-of-the-art baseline methods that
are designed specifically for connected subgraph detection, namely, GraphLaplacian
[33], EventTree
[32], DepthFirstGraphScan [36] and NPHGS [8].
DepthFirstGraphScan is an exact search algorithm
based on depth-first search and takes weeks to run on graphs
that have more than 1000 nodes. We imposed a maximum limit
on the depth of the search to 10 to reduce its time complexity.
The basic ideas of these baseline methods are summarized
as follows: NPHGS starts from random seeds (nodes) as
initial candidate clusters and gradually expends each candidate
cluster by including its neighboring nodes that could help
improve its BJ statistic score until no new nodes can be added.
The candidate cluster with the largest BJ statistic score is
returned. DepthFirstGraphScan adopts a similar strategy
to NPHGS but expands the initial clusters based on depth-first
search. GraphLaplacian uses a graph Laplacian penalty
function to replace the connectivity constraint and converts
the problem to a convex optimization problem. EventTree
reformulates the connected subgraph detection problem as a
prize-collecting steiner tree (PCST) problem [19] and apply
the Goemans-Williamson (G-W) algorithm for PCST [19]
to detect anomalous subgraphs. We also implemented the
generalized fused lasso model (GenFusedLasso) for graph
scan statistics using the framework of alternating direction
method of multipliers (ADMM). GenFusedLasso method
solves the following minimization problem
X
minn −f (x) + λ
|xi − xj |,
(26)
x∈R

(i,j)∈E

where f (x) is a predefined graph scan statistic and the
trade-off parameter λ controls the degree of smoothness of
neighboring entries in x. We applied the heuristic rounding
step proposed in [29] to x to generate connected subgraphs.
Parameter Tunning: We strictly followed strategies recommended by authors in their original papers to tune the
related model parameters. Specifically, for EventTree, we
tested the set of λ values: {0.02, 0.04, · · · , 2.0, 3.0, · · · , 20}.
For Graph-Laplacian, we tested the set of λ values:
{0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1} and returned the best result. For GenFusedLasso, we tested the set of λ values:
{0.02, 0.04, · · · , 2.0, 3.0, · · · , 20}. For NPHGS, we set the
suggested parameters by the authors: αmax = 0.15 and K =
5. Our proposed methods G RAPH -IHT and G RAPH -GHTP
have a single parameter k, an upper bound of the subgraph
size. We tested the set of k values: {50, 100, · · · , 1000}. As
the BWSN dataset has the ground truth about the contaminated
nodes, we identified the best parameter for each method that
has the largest F-measure in the training data. For the other
data sets, as we do not have ground truth about the true
subgraphs, for each specific scan statistic, we identified the
best parameter for each method that was able to identify the
connected subgraph with the largest statistic score.
Performance Metrics: 1) Optimization Power. The over-

all scores of the three graph scan statistic functions of the
connected subgraphs returned by the comparison methods are
compared and analyzed. The objective is to identify methods
that could find the connected subgraphs with the largest graph
scan statistic scores. 2) Precision, Recall, and F-Measure.
For the BWSN dataset, as the true anomalous subgraphs are
known,we use F-measure that combines precision and recall to
evaluate the quality of detected subgraphs by different methods. 3) Run Time. The running times of different methods
are compared.
B. Evolving Curves of Graph Scan Statistics
Figure 3 compares our methods Graph-IHT and G RAPH GHTP with GenFusedLasso on the scores of two graph
scan statistics (Kulldorff’s scan statistic and elevated mean
scan statistic (EMS)) based on the best connected subgraphs
identified by both methods in different iterations. Note that,
a heuristic rounding process as proposed in [29] is applied to
continuous vector xi estimated by GenFusedLasso at each
iteration i in order to identify the best connected subgraph at
the current iteration. As the setting of the parameter λ will
influence the quality of the detected connected subgraph, the
results under different λ values are also shown in Figure 2.
The results indicate that our method G RAPH -GHTP converges
in less than 10 iterations, and G RAPH -IHT converges in more
steps. The qualities (scan statistic scores) of the connected subgraphs identified at different iterations by our two methods are
consistently higher than those returned by GenFusedLasso.
C. Optimization Power
The comparisons between our method and the other baseline
methods are shown in Table III and Table IV. The scores
of the three graph scan statistics based on the connected
subgraphs returned by these methods are reported in these two
tables. The results in indicate that our method outperformed
all the baseline methods on the three graph scan statistics,
except that EventTree achieved the highest Kulldorff score
(16738.43) on the CitHepPh dataset, but is only 2.71% larger
than the returned score of our method G RAPH -GHTP. We note
EventTree is a heuristic algorithm and does not provide
theoretical guarantee on the quality of the connected subgraph
returned, as measured by the scan statistic scores.
D. Water Pollution Detection
Figure 4 shows the precision, recall, and F-measure of all
the comparison methods on the detection of polluted nodes in
the water distribution network in BWSN with respect to different noise ratios. The results indicate that our proposed method
G RAPH -GHTP and DepthFirstGraphScan were the best
methods on all the three measures for most of the settings.
However, DepthFirstGraphScan spent 5929 seconds to
finish, and G RAPH -GHTP spent only 166 seconds, 35.8
times faster than DepthFirstGraphScan. EventTree
achieved high recalls but low precisions consistently in different settings. In contrast, GraphLaplacian and NPHGS
achieved high precisions but low recalls in most settings.

Fig. 2. Scability of G RAPH -GHTP with respect to k Fig. 3. Evolving curves of graph scan statistic scores between our methods( G RAPH -IHT and
(the upper bound of subgraph size).
G RAPH -GHTP) and GenFusedLasso in different iterations.

(a) BWSN(precision)
Fig. 4.

(b) BWSN(recall)

(c) BWSN(fmeasure)

Precision, Recall, and F-measure curves for water pollution detection in BWSN with respect to different noise ratios.

TABLE III
Comparison on scores of the three graph scan statistics based on connected subgraphs returned by comparison methods. EMS and EBP are Elevated Mean
Scan Statistic and Expectation-Based Poisson Statistic, respectively.

G RAPH -GHTP
GraphLaplacian
EventTree
DepthFirstGraphScan
NPHGS

Kulldorff
1097.15
474.96
834.59
735.85
541.13

BWSN
EMS
EBP
21.56
79.71
14.89
49.91
20.25
32.13
20.41
79.30
16.90
58.59

E. Scalability Analysis
Table III and Table IV also show the comparison between
our proposed method G RAPH -GHTP and other baseline methods on the running time. The results indicate that our proposed method G RAPH -GHTP ran faster than all the baseline
methods in most of the settings, except for EventTree.
EventTree was the fastest method but was unable to detect
subgraphs with high qualities. As our method has a parameter
on the upper bound (k) of the subgraph returned, we also conducted the scalability of our method with respect to different
values of k as shown in Figure 2. The results indicate that the
running time of our algorithm is insensitive to the setting of
k, which is consistent with the time complexity analysis of
G RAPH -GHTP as discussed in Theorem III.5.

Run Time
165.86
55315.94
441.74
5929.00
256.91

Kulldorff
16296.40
2585.44
16738.43
9531.19
11965.14

CitHepPh
EMS
EBP
337.90
9342.94
202.38
2305.05
335.34
9061.56
260.06
5561.66
326.23
9098.22

Run Time
155.74
22424.24
124.28
12183.88
175.08

VI. R ELATED W ORK
A. Structured sparse optimization. The methods in this
category have been briefly reviewed in the introduction section. The most relevant work is by Hegde et al. [15]. The
authors present G RAPH-C OSAMP, a variant of C OSAMP [23],
for compressive sensing and linear regression problems based
on head and tail approximations of M(G, k).
B. Connected subgraph detection. Existing methods in this
category fall into three major categories:1) Exact algorithms.
The most recent method is a brunch-and-bounding
algorithm DepthFirstGraphScan [36] that runs in
exponential time in the worst case; 2) Heuristic algorithms.
The most recent methods in this category include
EventTree [32], NPHGS [8], AdditiveScan [37],

TABLE IV
Comparison on scores of the three graph scan statistics based on connected subgraphs returned by comparison methods. EMS and EBP are Elevated Mean
Scan Statistic and Expectation-Based Poisson Statistic, respectively. GraphLaplacian failed to run on ChicagoCrime due to out-of-memory error.

G RAPH -GHTP
GraphLaplacian
EventTree
DepthFirstGraphScan
NPHGS

EMS
20.45
5.40
12.40
8.13
6.28

Traffic
Run Time
22.25
291.75
5.02
47.73
0.22

GraphLaplacian
[33], and EdgeLasso
[34]; 3)
Approximation algorithms that provide performance bounds.
The most recent method is presented by Qian et al.. The
authors reformulate the connectivity constraint as linear matrix
inequalities (LMI) and present a semi-definite programming
algorithm based on convex relaxation of the LMI [18, 19] with
a performance bound. However, this method is not scalable
to large graphs (≥ 1000 nodes). Most of the above methods
are considered as baseline methods in our experiments and
are briefly summarized in Section V-A.
VII. CONCLUSION
This paper presents, G RAPH -IHT and G RAPH -GHTP, two
efficient algorithms to optimize a general nonlinear optimization problem subject to connectivity constraint on the
support of variables. Extensive experiments demonstrate the
effectiveness and efficency of our algorithms. For the future
work, we plan to explore graph-structured constraints other
than connectivity constraint and extend our proposed methods
such that good theoretical properties of the cost functions that
do not satisfy the WRSC condition can also be analyzed.
R EFERENCES
[1] B. Anbaroğlu, T. Cheng, and B. Heydecker. Non-recurrent traffic
congestion detection on heterogeneous urban road networks. TRANSPORTMETRICA, 11(9):754–771, 2015.
[2] M. Asteris and e. a. Kyrillidis. Stay on path: Pca along graph paths. In
ICML, pages 1728–1736, 2015.
[3] F. Bach. Learning with submodular functions: A convex optimization
perspective. arXiv:1111.6453, 2011.
[4] F. Bach, R. Jenatton, J. Mairal, G. Obozinski, et al. Structured sparsity
through convex optimization. Stat Sci, 27(4):450–468, 2012.
[5] S. Bahmani, P. T. Boufounos, and B. Raj. Learning model-based sparsity
via projected gradient descent. IT, 2016.
[6] S. Bahmani, B. Raj, and P. T. Boufounos. Greedy sparsity-constrained
optimization. JMLR, 14(1):807–841, 2013.
[7] T. Blumensath. Compressed sensing with nonlinear observations and
related nonlinear optimization problems. IT, 59(6):3466–3474, 2013.
[8] F. Chen and D. B. Neill. Non-parametric scan statistics for event
detection and forecasting in heterogeneous social media graphs. In ACM
SIGKDD, pages 1166–1175, 2014.
[9] F. Chen and D. B. Neill. Human rights event detection from heterogeneous social media graphs. Big Data, 3(1):34–40, 2015.
[10] J. W. Coulston and K. H. Riitters. Geographic analysis of forest health
indicators using spatial scan statistics. EM, 31(6):764–773, 2003.
[11] Dataset and code. https://github.com/baojianzhou/Graph-GHTP.
[12] W. L. Gorr and Y. Lee. Early warning system for temporary crime hot
spots. Journal of Quantitative Criminology, 31(1):25–47, 2015.
[13] C. Hegde, P. Indyk, and L. Schmidt. Approximation-tolerant modelbased compressive sensing. In SODA, pages 1544–1561. SIAM, 2014.
[14] C. Hegde, P. Indyk, and L. Schmidt. A fast approximation algorithm
for tree-sparse recovery. In ISIT, pages 1842–1846. IEEE, 2014.
[15] C. Hegde, P. Indyk, and L. Schmidt. A nearly-linear time framework
for graph-structured sparsity. In ICML, pages 928–937, 2015.

Kulldorff
6386.08
4388.42
1123.49
966.70

ChicagoCrime
EMS
EBP
5.45
5172.54
4.91
3965.96
2.56
1094.21
2.43
948.23

Run Time
3177.60
226.50
12133.50
701.40

[16] J. Huang, T. Zhang, and D. Metaxas. Learning with structured sparsity.
JMLR, 12:3371–3412, 2011.
[17] L. Jacob, G. Obozinski, and J.-P. Vert. Group lasso with overlap and
graph lasso. In ICML, pages 433–440. ACM, 2009.
[18] P. Jain, A. Tewari, and P. Kar. On iterative hard thresholding methods
for high-dimensional m-estimation. In NIPS, pages 685–693, 2014.
[19] D. S. Johnson, M. Minkoff, and S. Phillips. The prize collecting steiner
tree problem: theory and practice. In SODA, pages 760–769, 2000.
[20] J. e. a. Leskovec. Graphs over time: densification laws, shrinking
diameters and possible explanations. In KDD, pages 177–187, 2005.
[21] J. Mairal and B. Yu. Supervised feature selection in graphs with path
coding penalties and network flows. JMLR, 14(1):2449–2485, 2013.
[22] R. Modarres and G. Patil. Hotspot detection with bivariate data. Journal
of Statistical planning and inference, 137(11):3643–3654, 2007.
[23] D. Needell and J. A. Tropp. Cosamp: Iterative signal recovery from
incomplete and inaccurate samples. ACHA, 26(3):301–321, 2009.
[24] J. Neil and H. et al. Scan statistics for the online detection of locally
anomalous subgraphs. Technometrics, 55(4):403–414, 2013.
[25] D. B. Neill. An empirical comparison of spatial scan statistics for
outbreak detection. IJHG, 8(1):1, 2009.
[26] D. B. Neill. Fast subset scan for spatial pattern detection. JRSS: Series
B (Statistical Methodology), 74(2):337–360, 2012.
[27] D. Oliveira and D. P. et al. Detection of patterns in water distribution
pipe breakage using spatial scan statistics for point events in a physical
network. JCCE, 25(1):21–30, 2010.
[28] A. Ostfeld and J. G. e. a. Uber. The battle of the water sensor networks
(bwsn): A design challenge for engineers and algorithms. JWRPM,
134(6):556–568, 2008.
[29] J. Qian, V. Saligrama, and Y. Chen. Connected sub-graph detection. In
AISTATS, 2014.
[30] R. T. Rockafellar. Monotone operators and the proximal point algorithm.
SIAM journal on control and optimization, 14(5):877–898, 1976.
[31] L. A. Rossman and W. Supply. Epanet 2 users manual, 2000.
[32] P. Rozenshtein, A. Anagnostopoulos, A. Gionis, and N. Tatti. Event
detection in activity networks. In SIGKDD, pages 1176–1185, 2014.
[33] J. Sharpnack, A. Rinaldo, and A. Singh. Changepoint detection over
graphs with the spectral scan statistic. arXiv:1206.0773, 2012.
[34] J. Sharpnack, A. Rinaldo, and A. Singh. Sparsistency of the edge lasso
over graphs. In AISTATS, pages 1028–1036, 2012.
[35] J. L. Sharpnack, A. Krishnamurthy, and A. Singh. Near-optimal anomaly
detection in graphs using lovász extended scan statistic. In NIPS, pages
1959–1967, 2013.
[36] S. Speakman, E. McFowland Iii, and D. B. Neill. Scalable detection of
anomalous patterns with connectivity constraints. JCGS, 0(ja):00–00, 0.
[37] S. Speakman, Y. Zhang, and D. B. Neill. Dynamic pattern detection
with temporal consistency and connectivity constraints. In ICDM, pages
697–706. IEEE, 2013.
[38] A. Tewari and R. et al. Greedy algorithms for structurally constrained
high dimensional problems. In NIPS, pages 882–890, 2011.
[39] B. Xin and K. et al. Efficient generalized fused lasso and its application
to the diagnosis of alzheimers disease. In AAAI, pages 2163–2169, 2014.
[40] X. Yuan, P. Li, and T. Zhang. Gradient hard thresholding pursuit for
sparsity-constrained optimization. In ICML, pages 127–135, 2014.
[41] X.-T. Yuan and Q. Liu. Newton greedy pursuit: A quadratic approximation method for sparsity-constrained optimization. In CVPR, pages
4122–4129. IEEE, 2014.
[42] T. Zhang. Adaptive forward-backward greedy algorithm for sparse
learning with linear models. In NIPS, pages 1921–1928, 2009.

