Query-Driven Discovery of Anomalous Subgraphs in Attributed Graphs

1 Introduction

In recent years, graph-structure optimization for anomalous subgraph discovery in high-dimensional data and graphs as an open problem has attracted much attention [Qian et al., 2014; Wu et al., 2016; Chen and Zhou, 2016; Chen and Neill, 2014; Tong et al., 2007; McFowland et al., 2013; Gionis et al., 2015]. In many settings, anomaly in high-dimensional data presents a structure-specific shape, such as the cholera epidemic infection along a river [Patil et al., 2003]. To motivate this scenario, consider the Botnet Infection problem [Choi et al., 2009] as shown in Figure 1. Group activity presented in a "specific" shape is an inherent property of bots attacking on networks [Choi et al., 2009]. Therefore, we employ the prior knowledge about the structure of cyber attacks as the side information regularizing the generic nonlinear function optimization in attributed graphs as follows.

Given the graph G = (V,E) where V = [n] = {1, · · · , n} and E ⊆ V ×V , based on attributes (e.g., "transfer rate") for each vertex v ∈ V , the differentiable function ϕ is formulated as Rn → R. The support set of x is supp(x) = {i|xi ̸= 0}.

The domain for Problem (1) is the structure-specific model M = {Q1, · · · ,QK} where Qi ⊆ [n] is derived from G based on the "query graph". The model M is a family of structured support sets (e.g., its graphs are stars or trees).

Related work. Recent works solving Problem (1) fall into two main categories: 1) Shape induced. By encoding specific-shaped models (e.g., connectivity) [Bach et al., 2011] as structured sparsity-inducing norms, the methods can reformulate Problem (1) as a convex (or non-convex) optimization problem: minx∈Rn ϕ(x) + λ · Ω(x) where λ is a trade-off parameter and Ω(x) is a structured sparsityinducing norm of M that is typically non-smooth and nonEuclidean. These methods just focus on a connecting property without the specific-shape. 2) Model-projection based. The methods depend on a projection oracle of M: P(b) = argminx∈Rn ||x - b||22 subject to supp(x) ∈ M, and decompose Problem (1) into two subproblems: minimizing the unconstrained ϕ(x) and the projection oracle P(b). There are many methods to solve Problem (1) under the assumption that P(b) can be addressed exactly [Bahmani et al., 2013; Jain et al., 2014; Yuan et al., 2014; Bahmani et al., 2016; 2013; Yuan and Liu, 2014; Hegde et al., 2015]. These methods are only applicable to quadratic functions subject to linear constraints. However, few works have investigated combinations of generic nonlinear functions and specific-shaped constraints.

We present an approach for optimizing generic nonlinear cost functions subject to structure-specific constraints, where the constraints are formalized to query graphs. Then the output subgraph is isomorphic to the query graph. By a proper choice of the cost functions ϕ(·), our approach can be implemented easily in different types of attributed graphs, such as water sensor and emergency department networks. Main contributions of this study are summarized as follows:

• Design of an efficient structure-specific subgraph discovery algorithm. A new algorithm is proposed to approximately solve Problem (1) where ϕ is a differentiable nonlinear cost function. Over the structurespecific model M, our proposed algorithm is required to minimize ϕ(x) over a projected subspace supported by the tree-shape projection oracle.
• Theoretical guarantees. The convergence rate and accuracy of our proposed algorithm are analyzed under the Weak Restricted Strong Convexity condition, which is more general than the popular Restricted Strong Convexity (RSC) condition. We prove that our proposed algorithm enjoys rigorous theoretical guarantees.
• Compressive experiments to validate our approach from the effectiveness and efficiency perspectives. Our proposed algorithm is suitable for optimization of a variety of graph scan statistics as the target for anomaly specific-shaped subgraph detection. Extensive experiments on a number of benchmark datasets demonstrate that the algorithm performs better than the representative methods for this task on both accuracy and run time.

2 ProblemFormulation

Before the problem formulation for the anomalous subgraph discovery in attributed graphs, we illustrate the three key concepts by a case of anomalous subgraph discovery in Figure 2. Attributed Graph. A network (e.g., computer network) can be modeled as an attributed graph composed of the underlying graph G = (V,E) and attribute matrix W ∈ Rn×p , where n, p denote the number of vertices and attributes, V = [n] and E ⊆ V × V , (e.g., Wv,0 denotes the "transfer rate" of the vertex v at the attribute "0").

Query Graph. Given a query Q = (V, E), we focus on the query having a structure-specific shape (e.g., tree shape).

Match Model. A graph C is a subgraph of G, denoted as C ⊆ G, if VC ⊆ VG, EC ⊆ EG and ∀(u, v) ∈ EC, u, v ∈ VC. A graph C is isomorphic to a query graph Q, denoted as C ∼= Q, if there is a bijection ψ : VC → VQ such that, for every pair of vertices u, v ∈ VC, (u, v) ∈ EC if and only if (ψ(u), ψ(v)) ∈ EQ. Given the underlying graph G and the query graph Q, we present the structure-specific domain

The valid domain M(Q) illustrates all vertex subsets of G whose shapes are same to the query graph Q.

Anomalous Structure-Specific Subgraph Discovery. We aim to obtain S ∈ M(Q) for minimizing the cost function over the domain. Given G, we assume that VG = {1, · · · , n}. The vector form of VG is denoted as x ∈ Rn , and the support set of x in VG is denoted as supp(x) = {i | xi ̸= 0}. Thus this problem is formulated as

where ϕ(·) is a differentiable nonlinear cost function modeled from the attribute matrix W. Figure 2 illustrates our problem through an attributed computed network and a query graph. There exist multiple matchings of the query in the network.

3 Preliminary: Tree-Shape Projection Oracle

For Problem (3), the gradient or solution vector of the function may fall into the invalid domain. We present the work on projecting a vector variable into the valid domain.

Given the query graph Q, the underlying graph G, the match model M(Q), and the vector g, the projection oracle, P : Rn → VG, is defined as

where gS denotes the restriction of g to indices S so that (gS)i = gi for i ∈ S and (gS)i = 0 otherwise. There are many works about solving Problem (4) with high time complexity. For a general query graph, we present a simple oracle to Problem (4) in Figure 3. This oracle consists of two main steps, query decomposition and projection.

Query decomposition. Wu et al., looked at approximating arbitrary query graph metrics by a tree shaped prior (e.g., spanning tree) [Wu et al., 2016]. By a prior, we generate a spanning tree of the query graph. From the root of the tree, we collect a set of query paths terminated with leaf nodes.

Projection. Each vertex v ∈ VG is assigned with the value gv. We search the maximum-weight paths from a root vertex v in G for all of query paths, on condition that the assembled paths are isomorphic to the tree of the query graph. The condition is easily satisfied by caring for "share nodes" (i.e., a share node exists in multiple query paths). Within all paths for v ∈ VG, we select the maximum weight assembled paths as the target structure-specific graph, whose set of nodes is S.

4 Methodology

We propose a new graph-structured optimization method for nonlinear cost functions and discuss its theoretical properties.

4.1 Graph Tree-shape Projection Pursuit Algorithm (Graph-TPP)

The algorithm described here iterates on the solution vector for Problem (3). For each iteration, we first obtain the most significant support set within the valid domain for the gradient of the cost function at the current solution. Then we refine the solution over the significant support set of gradient and the previous solution. Graph-TPP is illustrated in Algorithm 1 in detail, where each iteration consists of six steps.

Projection. For a gradient ∇ϕ(xi ) or a solution vector b, we identify the most significant support set (i.e., the set maximizing Problem (4)) within the valid domain.

Support. First we calculate the projected gradient descent at the current xi with step-size η (1 by default). Then the support set Ω is obtained and can be interpreted as the subspace where the nonconvex set {x | supp(x) ⊆ Ω} is located.

Estimate. Over the support set Ω, let Ωc = VG \ Ω, the function ϕ is minimized to make an intermediate estimate b.

Prune. This step (Line 7) calculates the next solution vector xi+1 : xi+1 = bS, which retains the most significant structured-entries in b.

We have two popular options for defining the halting criterion: (1) the difference between the score functions is less than a threshold |ϕ(xi )-ϕ(xi+1 )| < ϵ; and (2) the difference between the vectors is less than a threshold ||xi -xi+1 ||2 < ϵ (e.g., ϵ = 0.01).

Discussion: The projection problem is hard to solve due to the NP-hard subgraph isomorphism problem. For connected subgraph detection, work [Wu et al., 2016] yielded compelling results by tree shape priors (e.g., spanning tree). For a general query graph, we employ the tree-shape projection oracle to perform graph matching, which is different from the previous work that has head and tail projections [Chen and Zhou, 2016]. We focus on the isomorphic subgraph discovery; however, the work [Chen and Zhou, 2016] focuses on the connected subgraph discovery.

4.2 Theoretical Analysis

The theoretical properties of Graph-TPP are examined in two aspects: Studying convergence rate, and time complexity. Before obtaining the theoretical results, we require the following

Algorithm 1: Graph-TPP

Input: Attribute W, underlying graph G, query graph Q and step size η (1 by default)

Result: Vertex set of structure-specific subgraph S

1 Set i = 0, xi = 0;

2 repeat

3 Γ = P(∇ϕ(xi ));

▷ Projection

4 Ω = supp(xi - η∇Γϕ(xi ));

▷ Support

5 b = argminx∈Rn ϕ(x) s.t. xΩc = 0; ▷ Estimate

6 S = P(b);

▷ Projection

7 xi+1 = bS;

8 i ← i + 1;

• 9 until halting condition holds;

10 return S;

▷ Prune

key technical condition, under which the results are guaranteed. Let s = |VQ|.

Definition 1 (Weak Restricted Strong Convexity Property (WRSC)). Let M = {S | S ⊂ [n], |S| ≤ 2s}, and there must be S ∈ M if S ∈ M(Q). The differentiable function ϕ has the condition (ξ, δ,M)-WRSC if ∀x, y ∈ Rn , ∀S ∈ M with supp(x)∪supp(y) ⊆ S, for some ξ > 0 and 0 < δ < 1.

Theorem 1 (Graph-TPP Convergence). Given an attributed graph W, G and a query graph Q, we obtain the domainM. Consider the differentiable cost function ϕ : Rn → R that satisfies the condition (ξ, δ,M)-WRSC. Let xi ∈ Rn . Then for any true x∗ ∈ Rn and supp(x∗ ) ∈M(Q), the iteration of Algorithm 1 holds

The estimation error of Graph-TPP depends on the multipliers of ||∇Sϕ(x∗ )||2. Before reaching the estimation error level, Graph-TPP geometrically converges to the nearoptimal x∗ . Especially, when ||∇Iϕ(x∗ )||2 = 0, Graph-TPP guarantees that the true x∗ is obtained within finite iterations.

Theorem 2 (Graph-TPP Time Cost). Given α < 1, supp(x∗ ) ∈ M(Q) and the function ϕ satisfies the WRSC condition, Graph-TPP in Algorithm 1 returns a ˆx such that ||x∗ - ˆx||2 ≤ (1 + β 1-α )||∇Iϕ(x∗ )||2. The estimation error is less than a constant. Graph-TPP runs in time

where T is the time cost of execution for one iteration in Algorithm 1. If T can scale linearly with n, then Graph-TPP scales nearly linearly with n. (Proof: see Appendix)

4.3 Application to Graph Scan Statistics

We specialize ϕ to be a number of graph scan statistics that are widely employed in pattern detection in graphs.

Graph Scan Statistic. First select a graph C ⊆ G and C ∼= Q. Then a graph scan statistic is defined as F(VC) = log Prob(data|H1(VC)) Prob(data|H0) that is the log-likelihood ratio statistic over all matching subgraphs of G. For the null hypothesis H0, we assume no anomaly structure-specific subgraphs (i.e., all the observed data are generated from the expected distribution). However, for the alternative hypothesis H1(VC), we assume the observed data in the subgraph C show a significant increase in, such as the number of visiting patients, for some multiplicative factors. At last the detection problem is formulated as: minC⊆G -F (VC) s.t. VC ∈ M(Q).

Let ϕ(x) = -F(VC). Denote the vector form of VC as x ∈ {0, 1}n , such that supp(x) = VC where x can be relaxed to x ∈ [0, 1]n [Chen and Zhou, 2016]. We mainly examine ϕ in Equation (8) to log forms of expectation-based Poisson (EBP) and Kulldorff's (KULL) graph scan statistics. Let c = W:,0 denote the "observed count" attribute values. Similarly, let d =W:,1 denote the "expected count" attribute values.

where x ∈ [0, 1]n . Let ˆ r = min{di/ci} for i ∈ VG. By [Yuan and Liu, 2014], ϕEBP satisfies the WRSC condition that δ = √1 - 2ξ(1 - ˆ r2 ) + ξ 2 for ξ < 2(1 - ˆ r2 ).

5 Experiments

We compare the effectiveness and efficiency of our method Graph-TPP with respect to competitive methods.

5.1 Experimental Design

1) Respiratory Emergency Department (ED) Dataset. Given a grid network which consists of 10,000 nodes and 14,850 edges, we consider each node as a zip code. For each node v ∈ V , we collected the T day period of respiratory ED visit data Dv ∈ RT (e.g., T = 28), where especially the time t = 0 denotes the current day. During non-outbreak period, the number of patients visiting ED in v is simulated with the Poisson distribution [Neill, 2009], i.e., Dtv ← Poisson(µ) for t = 0, · · · , T where µ denotes the expected number of patients visiting ED in v on those days and µ is randomly selected from {1, · · · , 34}. We randomly select the outbreak duration of U days from {1, · · · , 7}. On each day t ∈ {0, · · · , U} of the outbreak, we inject Poisson((T - t)wv∆) cases into each infected node v where wv ∝ ∑t Dtv is normalized so that the total weight equals to 1 in infected zips and ∆ denotes the outbreak severity (e.g., ∆=800), i.e., Dtv ← Dtv + Poisson((T - t)wv∆) for t = 0, · · · ,U. We performed a simulation for this dataset with medium-size injects affecting 10 percent nodes, which consist of the groundtruth connected target subgraph [Neill, 2009]. We simulated 500 graph snapshots. For testing the robustness of methods to noise, we randomly selected K ∈ {2, 4, 6, 8, 10} percent nodes in this network, and flipped their values (i.e., no inject outbreak cases if the nodes are infected, inject outbreak cases otherwise). Let cv = D0 v and dv = 1 T ∑Tt=1 Dtv . 2) Water Pollution Data. These data were collected on the real-world network of 12,527 nodes and 14,831 edges, and there were four nodes with chemical contaminant plumes, distributed in four different areas [Chen and Zhou, 2016]. "The spreads of contaminant plumes were simulated using the water network simulator EPANET for 8 hours" [Chen and Zhou, 2016]. For each hour, the value at each node v was reported with the corresponding sensor, Dv ← 1 if it is polluted and Dv ← 0 otherwise. We randomly selected K percent nodes, and flipped their sensor binary values for testing the robustness of methods to noises, where K ∈ {2, 4, 6, 8, 10}. We considered the assigned attributed graphs as cv ← Dv (i.e., the report of the sensor at the node v) and dv ← K% (i.e., the noise ratio) [Chen and Zhou, 2016]. 3) Real-World Network Dataset. An Internet company 1 provided a traffic network of 131,107 nodes and 358,386 edges derived from ∗edu.cn web sites browsing logs. The total 3,978,073 logs were collected fromMay 31, 2014 to May 13, 2015. For a day t during this period and a node v in this network, we denoted the number of logs within v on that day t as the observed value, and the average number of logs within v before t as the expected value. We have a graph snapshot for each day.

Comparison Methods. We considered two state-of-theart methods: Top-K [Gupta et al., 2014] and Fast-K [Yang et al., 2016]. The main idea behind the methods is as follows: Given a query graph Q and an attributed graph G, the goal is to find top k subgraphs that are matched to Q and have the highest ranking scores by a ranking function (e.g., summing vertex weights for Fast-K and edge weights for Top-K). From the k subgraphs, we selected the top one subgraph as the anomalous subgraph. By author recommendations, set k = 10 for Top-K, and k = 20 for Fast-K. Let the vertex weight w(v) ← Wv,0. As the Top-k method is applied to edge weights, we replicated each node v to v′ and added edge (v, v′ ) to G. We set the weight of edge (v, v′ ) to w(v) for a pair of original node and replicate node, and set the weight of edge (v, v′ ) to 0 otherwise. For Top-k method, return the subgraph by removing replicate nodes and corresponding edges.

Performance Metrics. 1) Effectiveness: Precision. We compute the precision of the target subgraph (i.e., the ratio of the number of correct anomalous nodes and the number of nodes). As the size of returned target subgraph is fixed, we ignore the recall metric. 2) Efficiency: Graph Scan Statistic Score and Running Time. The optimization power of our method is analyzed through the scores and running times. 3) Case Study. We present some cases to illustrate our method.

5.2 Effectiveness

As the size of detected subgraph is fixed for the graph isomorphism constraint, we focus on the precision of methods in Figure 4. By evaluating methods to three typical query graphs in the left sub-figure of Figure 4 [Du and Yang, 2011; Yang et al., 2014; Kim et al., 2015; Pan et al., 2013], we can observe precision in the middle sub-figure of Figure 4. The query graph Q1 shows the star-shape, which means that a node infects its neighbors, and Q2 denotes a many-infectmany case. The query graph Q3 denotes the infection from one star-shape area to another neighbor star-shape area.

1 An Internet security company in China with more than 0.6 billion users.

In Figure 4, we present a comparison of precision for methods in these two data sets in detail. For the emergency data, at the 2 percent noise level, our proposed Graph-TPP (i.e., ϕEBP and ϕKULL) achieved higher precision (close to 1) than competitive baselines (close to 0.9). However, even at 10 percent noise level, our methods achieved at least 0.58 precision, and baselines achieved the best precision to 0.5. Especially for Q1 and Q3, the precision of baselines decreases sharply to about 0.43 with the noise increasing, but the precision of ϕEBP is greater than 0.68. For the water pollution data, Graph-TPP and Fast-K perform similarly in Q1 and Q2, but perform better than Top-K. However, in Q3, the worst precision 0.88 for Graph-TPP (ϕEBP ) is greater than the best precision 0.87 for baselines. Our methods outperformed all the baselines on precision in both of datasets.

5.3 Efficiency

Right sub-figures of Figure 4 depicts the scores of graph scan statistics (EBP and KULL) changing with the iteration. The score corresponds to the graph scan statistic F(VC). We can observe that our proposed Graph-TPP has converged at most five iterations. According to Theorem 1, our methods geometrically converge to the near-optimal stationary point.

The table of Figure 4 shows the time cost of all competitive methods on the two datasets. We can observe that the running times of Graph-TPP (EBP and KULL) are approximately equal in the two datasets. Although Fast-K performed well on the run time, it is a heuristic algorithm. Our methods enjoy a rigorous theoretical guarantee. In particular, for the water pollution data, the run time of Graph-TPP is close to the Fast-K method.

5.4 Case Study

We tested our methods Graph-TPP on the real-world network dataset for the star and bipartite shaped attack patterns.

Star-shaped attack case. From the left star-shaped subfigures of Figure 5(b-e), the left subfigure of Figure 5(b) is a Dedecms Attack 2 network, and the left subfigure of Figure 5(e) is a FckEditor Attack 3 network. Our methods successfully discovered the cyber-attack networks without innocent nodes. The difference between the two networks is that the client X.X.74.39 attacked the other server sites ∗.hlu.edu.cn and 0517edu.cn; however, the server site www.hlu.edu.cn was attacked by four clients. These cyberattack patterns are the most common forms in networks. For the left subfigures of Figure 5(c, d), our method detected the attack client X.X.176.66 and most of the attacked server sites.

Bipartite-shaped attack case. In Figure 5(b-e), the right subfigures are bipartite-shaped cyber-attack networks. In the right subfigures of Figure 5(b-c), Graph-TPP (EBP and KULL) detected the same cyber-attack network. In the right subfigures of Figure 5(d-e), Graph-TPP (KULL) detected the cyber-attack network without innocent nodes, and GraphTPP (EBP) detected the cyber-attack network with one innocent node X.X.55.47. As the client X.X.55.47 is replaced by X.X.176.22, a new server site www.hlu.edu.cn attacked by the four clients is discovered. For this case, Graph-TPP (KULL) performs better than Graph-TPP (EBP).

In Figure 5(b-e), the variations of X.X.176.X contributed to many cyber-attacks. The server site www.ncyzedu.cn is the innocent node in the right subfigures of Figure 5(b-c), and this server site is not successfully attacked for removing the common vulnerability. From Figure 5, we can observe that Graph-TPP (EBP and KULL) works well on different scenarios. In this paper, we present a feasible framework GraphTPP to optimize a generic nonlinear function on the structurespecific query graph constraints.

6 Conclusions

We present an efficient algorithm to optimize a general nonlinear function on the structure-specific constraints. For future work, we will extend our work on large graphs, other similarity measures (e.g., exploiting domain semantics), and other types of attributes (e.g., categorical attributes).

Acknowledgments

The work is supported by China 973 Fundamental R&D Program (No.2014CB340300), NSFC program (No.61472022, 61421003), SKLSDE-2016ZX-11, and partly by the Beijing Advanced Innovation Center for Big Data and Brain Computing. The corresponding author is Jianxin Li.

A Proof of Theorem 1

2 Generate sql injection attacks by Common Vulnerability.

3 Generate FckEditor attacks by a vulnerable FckEditor version.

where as the support set of optimal x∗ can be built from the space M(Q), bS is restricted to the elements in b with the maximal score. Thus we have ||bS - b||2 ≤ ||b - x∗ ||2. Next we compute the upper bound of the component ||x∗Ω - bΩ||2, ||x∗Ω - bΩ||22 =< x∗ - b, x∗Ω - bΩ >

where ∇Ωϕ(b) = 0 since b is minimized at Line 5 of Graph-TPP, and Property (5) is applied. Then we have ||x∗Ω - bΩ||2 ≤ δ||x∗ - b||2 + ξ||∇Ωϕ(x∗ )||2.

Let Φ = supp(x∗ ) ∈ M(Q). As Ω = supp(xi - η∇Γϕ(xi )), we obtain the fact as follows. Eliminate the intersection Φ ∩ Ω, we obtain

For the right-hand formula, we have

where ||x∗Ω\Φ ||2 = 0. For the left-hand formula, we have

where x∗Φ\Ω = x∗Ωc and bΩc = 0. Let Π = Φ ∪ Ω - Φ ∩ Ω be the symmetric difference of the sets Φ and Ω. We have that

where the first inequality follows from the fact that, for a ≥ 0 and b ≥ 0, (√a + √b)2 ≤ a + b + 2√ab ≤ 2(a + b)

Let Ψ = supp(ri ). Inequalities on the other side,

where the third inequality satisfies the condition-WRSC and riΓ∪Ψ = ri . With the two bounds, we obtain the inequality.

By applying the similar algebraic techniques in [Hegde et al., 2014] Page 11, we obtain

B Proof of Theorem 2

Proof. The i-th solution xi in Algorithm 1 satisfies

After the t = ⌈ log (||x∗ ||2/||∇Iϕ(x∗ )||2 )/ log(1/α)⌉ iterations, Graph-TPP returns a desired estimate ˆx. As T is the time cost of execution for one iterate in Algorithm 1, and the total iterations are ⌈ log (||x∗ ||2/||∇Iϕ(x∗ )||2 )/ log(1/α)⌉, we prove this theorem over.