Concept Drift and Anomaly Detection in Graph Streams

TL;DR

This paper introduces a general methodology for detecting concept drift and anomalies in graph streams by embedding graphs into a vector space and applying change detection, validated through theoretical proofs and experiments.

Contribution

The paper proposes a novel, general approach for change detection in attributed graph streams without requiring vertex correspondence, supported by theoretical analysis and empirical validation.

Findings

Effective detection of concept drift in graph streams.

Outperforms baseline methods in biological and drawing graph datasets.

Theoretical guarantees for the proposed change detection method.

Abstract

Graph representations offer powerful and intuitive ways to describe data in a multitude of application domains. Here, we consider stochastic processes generating graphs and propose a methodology for detecting changes in stationarity of such processes. The methodology is general and considers a process generating attributed graphs with a variable number of vertices/edges, without the need to assume one-to-one correspondence between vertices at different time steps. The methodology acts by embedding every graph of the stream into a vector domain, where a conventional multivariate change detection procedure can be easily applied. We ground the soundness of our proposal by proving several theoretical results. In addition, we provide a specific implementation of the methodology and evaluate its effectiveness on several detection problems involving attributed graphs representing biological molecules and drawings. Experimental results are contrasted with respect to suitable baseline methods, demonstrating the effectiveness of our approach.

Tables

Table 1. Table I: Experimental settings. The first column contains an identifier for each experiment; in particular, in the Letter dataset, ‘D’, ‘O’, and ‘S’ stand for disjoint, overlapping, and subset, respectively. The second column reports the dataset involved, and the third and fourth columns show the set of classes from which nominal/non-nominal graphs are extracted. The collections of letters are selected in alphabetical order.

		Collection
Experiment ID	Dataset	Nominal	Non-nominal
L-D2	Letter	A,E	F,H
L-D5	Letter	A,E,F,H,I	K,L,M,N,T
L-O	Letter	A,E,F,H	F,H,I,K
L-S	Letter	A,E,F,H,I	F,H,I
MUT	Mutagenicity	non-mutagenic	mutagenic
AIDS	AIDS	inactive	active

Table 2. Table II: Results attained by our method varying number of prototypes ( M 𝑀 M ) and window size ( n 𝑛 n ). Symbol ∗ *\; indicate DCR statistically larger than 0.99, with a confidence of 95%. Symbols ∘ \circ\; and ∙ ∙ \bullet\; highlight significant trends in the DCR increasing M 𝑀 M and n 𝑛 n , respectively.

Experiment			DCR		ARL0		DoD		FA1000
dataset	$M$	$n$	mean	95CI	mean	95CI	mean	95CI	mean	std
L-D2	4	5	$\ast$1.000	[1.000, 1.000]	189	[103, 333]	26	[8, 66]	1.135	0.368
L-D5	4	5	0.380	[0.290, 0.480]	175	[91, 271]	411	[9, 2149]	1.207	0.390
L-D5	4	25	0.970	[0.930, 1.000]	183	[88, 306]	41	[1, 184]	0.233	0.086
L-D5	4	125	$\ast$1.000	[1.000, 1.000]	305	[61, 1148]	2	[1, 5]	0.045	0.038
L-D5	8	25	0.990	[0.970, 1.000]	175	[93, 372]	10	[1, 54]	0.245	0.078
L-D5	8	125	$\ast$1.000	[1.000, 1.000]	255	[50, 928]	1	[1, 1]	0.054	0.047
L-O	4	5	0.230	[0.150, 0.310]	191	[120, 344]	410	[21, 1136]	1.096	0.318
L-O	4	25	$\circ$0.790	[0.710, 0.870]	205	[101, 355]	123	[3, 723]	0.205	0.070
L-O	4	125	0.940	[0.890, 0.980]	291	[55, 725]	37	[1, 405]	0.045	0.035
L-O	8	25	$\circ$0.980	[0.950, 1.000]	163	[92, 306]	24	[2, 132]	0.260	0.080
L-O	8	125	$\ast$1.000	[1.000, 1.000]	238	[52, 620]	2	[1, 5]	0.052	0.053
L-S	4	5	0.720	[0.630, 0.810]	164	[90, 359]	132	[33, 349]	1.321	0.418
L-S	4	25	$\ast$1.000	[1.000, 1.000]	193	[103, 372]	28	[4, 123]	0.223	0.075
AIDS	4	5	$\ast$1.000	[1.000, 1.000]	175	[100, 337]	1	[1, 1]	1.221	0.364
MUT	4	5	$\bullet$0.050	[0.010, 0.100]	37	[18, 89]	70	[39, 124]	6.294	2.159
MUT	4	25	$\bullet$0.800	[0.720, 0.880]	52	[19, 145]	53	[4, 348]	1.075	0.521
MUT	4	125	$\bullet$0.980	[0.950, 1.000]	153	[8, 820]	5	[1, 22]	0.276	0.284
MUT	8	25	0.920	[0.860, 0.970]	42	[15, 91]	15	[2, 75]	1.242	0.680
MUT	8	125	$\ast$1.000	[1.000, 1.000]	112	[9, 540]	2	[1, 4]	0.292	0.267

Table 3. Table III: Results attained by baseline methods (index column) based on the graph density (den), the spectral gap of the Laplacian (SG), and the degenerate implementation of the methodology with M = 1 𝑀 1 M=1 prototype (M1). Symbol † † \dagger\; indicates significantly better (95% confidence) results.

Experiment			DCR		ARL0		DoD		FA1000
dataset	index	$n$	mean	95CI	mean	95CI	mean	95CI	mean	std
L-D2	M1	25	1.000	[1.000, 1.000]	238	[106, 594]	15	[2, 87]	0.194	0.082
L-D2	den	1	1.000	[1.000, 1.000]	194	[68, 556]	41	[22, 104]	6.475	3.337
L-D2	SG	1	0.850	[0.780, 0.920]	210	[75, 394]	133	[55, 274]	6.460	3.342
L-D5	M1	25	0.780	[0.700, 0.860]	222	[96, 407]	85	[2, 542]	0.200	0.098
L-D5	den	1	0.180	[0.110, 0.260]	219	[70, 485]	369	[91, 1147]	5.775	3.346
L-D5	SG	1	$†$1.000	[1.000, 1.000]	209	[74, 674]	33	[18, 53]	6.058	2.971
L-O	M1	25	0.500	[0.400, 0.600]	218	[100, 459]	304	[7, 1199]	0.203	0.082
L-O	den	1	$†$1.000	[1.000, 1.000]	194	[70, 549]	5	[4, 5]	6.853	3.428
L-O	SG	1	0.130	[0.070, 0.200]	221	[89, 613]	352	[115, 973]	5.558	2.555
L-S	M1	25	0.880	[0.810, 0.940]	230	[116, 469]	93	[4, 636]	0.189	0.066
L-S	den	1	$†$1.000	[1.000, 1.000]	188	[79, 385]	32	[20, 48]	6.237	2.777
L-S	SG	1	0.080	[0.030, 0.140]	189	[91, 455]	308	[124, 771]	6.100	2.247
AIDS	M1	25	1.000	[1.000, 1.000]	229	[95, 489]	1	[1, 1]	0.198	0.081
AIDS	den	1	1.000	[1.000, 1.000]	207	[67, 505]	4	[2, 6]	6.312	3.245
AIDS	SG	1	1.000	[1.000, 1.000]	197	[79, 443]	93	[48, 152]	6.077	2.744
MUT	M1	25	0.260	[0.180, 0.350]	154	[24, 575]	484	[13, 2208]	0.506	0.441
MUT	den	1	0.230	[0.150, 0.310]	194	[69, 381]	267	[96, 741]	6.225	3.175
MUT	SG	1	0.030	[0.000, 0.070]	208	[72, 437]	688	[134, 2342]	5.940	2.997

Equations

$$g_t\sim \{\begin{array}{ll}Q& t<\tau ,\\ \tilde{Q}& t\ge \tau .\end{array}$$

$$y=\zeta (g):={\left[d(g,r_1),\dots ,d(g,r_M)\right]}^{\top }.$$

$$F(B)=Q({\zeta }^{-1}(B)),\forall B\in \mathcal{B}.$$

$$\begin{array}[]{ll}H_{0}:\quad\mathbb{E}[S_{t}]=0\\ H_{1}:\quad\mathbb{E}[S_{t}]>0\end{array}$$

$$S_t>h_t\Rightarrow \text{a change is detected}$$

$$\hat{\tau }=\inf \{t:S_t>h_t\}.$$

$$\begin{array}[]{ll}H_{0}:&d^{\prime}(\mathbb{E}[F],\overline{\mathbf{\boldsymbol{y}}})\leq\delta\\ H_{1}:&d^{\prime}(\mathbb{E}[F],\overline{\mathbf{\boldsymbol{y}}})>\delta\end{array}$$

$$\begin{array}[]{ll}\alpha&=\mathbb{P}(d(\mu[Q],\mu[\mathbf{\boldsymbol{g}}])>\gamma|H_{0}),\\ \alpha^{\prime}&=\mathbb{P}(d^{\prime}(\mathbb{E}[F],\overline{\mathbf{\boldsymbol{y}}})>\gamma^{\prime}|H_{0}).\end{array}$$

$$d(g,f)$$

$$\begin{array}[]{cl}\left\|\mathbb{E}[F]-\zeta\left(\mu[Q]\right)\right\|_{{2}}^{2}&\leq\ v_{2}\\ \mathbb{P}\left(\left\|\overline{\mathbf{\boldsymbol{y}}}-\zeta(\mu[\mathbf{\boldsymbol{g}}])\right\|_{{2}}^{2}\geq\delta\right)&\leq\ \frac{v_{2}}{\delta}\qquad\forall\delta>0.\end{array}$$

$$\mathbb{P}\left(d_1(x)\le u(d_2(x))\right)\ge p\Rightarrow {\Phi }_1(\cdot )\ge p\cdot {\Phi }_2(u^{-1}(\cdot )).$$

$$p_{\delta }\cdot {\rm Y}\left(\frac{\gamma }{C}-b_{\delta }\right)\le \Psi (\gamma )\le \frac{1}{p_{\delta }}\cdot {\rm Y}\left(\frac{\gamma }{c}+b_{\delta }\right).$$

$$\mathbb{P}\left(d^{\prime }(\mathbb{E}[F],\overline{\mathbf{y}})\le {\gamma }^{\prime }|H_0\right)={\rm Y}\left(\frac{\gamma }{c}+b_{\delta }\right)\ge p_{\delta }(1-\alpha ).$$

$$\mathbb{P}\left(d(\mu [Q],\mu [\mathbf{g}])\le \gamma |H_0\right)=\Psi (\gamma )\ge p_{\delta }(1-{\alpha }^{\prime }).$$

$$R=\underset{R\in \{R^{(i)}\}}{\mathrm{argmin}}\left\{\underset{r_m\in R}{\max }\left[\underset{c\in C_m}{\max }d(r_m,c)\right]\right\},$$

$$s_w:=d^{\prime }({\overline{\mathbf{y}}}_0,{\overline{\mathbf{y}}}_w).$$

$$d_{\Sigma }({\overline{\mathbf{y}}}_0,{\overline{\mathbf{y}}}_w)=\sqrt{({\overline{\mathbf{y}}}_0-{\overline{\mathbf{y}}}_w{)}^{\top }{\Sigma }^{-1}({\overline{\mathbf{y}}}_0-{\overline{\mathbf{y}}}_w)},$$

$$\{\begin{array}{l}{S}_w=\max \left\{0,S_{w-1}+(s_w-q)\right\}\\ S_0=0.\end{array}$$

$$d(g,f)\ge \sqrt{\frac{{\lambda }_M}{M}}{d}_{\Sigma }(\zeta (g),\zeta (f)),$$

$$c{d}_{\Sigma }(u_1,u_2)\le d_F(g_1,g_2)\le C{d}_{\Sigma }(u_1,u_2),$$

$$\{O(\rho ,g)|\rho >0,g\in \mathcal{G}\}\cup \{{\zeta }^{-1}(B)|B\in \mathcal{B})\}$$

$$\mathbb{E}[{\mathrm{V}}_{\mathrm{f}}[\mathbf{x}]]=\left(1-\frac{1}{n}\right){\mathrm{V}}_{\mathrm{f}}[P].$$

$${\mathcal{F}}_P(z)={\mathcal{F}}_P(\mathbb{E}[P])+{\Vert \mathbb{E}[P]-z\Vert }_2^2.$$

$${\mathcal{F}}_{\mathbf{x}}(z)={\mathcal{F}}_{\mathbf{x}}(\overline{\mathbf{x}})+{\Vert \overline{\mathbf{x}}-z\Vert }_2^2,$$

$$\bullet$$

$${\mathrm{V}}_{\mathrm{f}}[P]$$

$$=\mathbb{E}\left[{\Vert x\Vert }_2^2-2y^{\top }\mathbb{E}[P]+{\Vert \mu [P]\Vert }_2^2\right]$$

$$=\mathbb{E}\left[{\Vert x\Vert }_2^2\right]-{\Vert \mu [P]\Vert }_2^2.$$

$${\mathrm{V}}_{\mathrm{f}}[\mathbf{x}]$$

$$=\frac{1}{n}\sum _t{x}_t^{\top }{x}_t-\frac{1}{n^2}\sum _{i,j}{x}_i^{\top }{x}_j$$

Full text

Nomenclature

[TABLE]

I Introduction

Learning in non-stationary environments is becoming a hot research topic, as proven by the increasing body of literature on the subject, e.g., see [1, 2] for a survey. Within this learning framework, it is of particular relevance the detection of changes in stationarity of the data generating process. This can be achieved by means of either passive approaches [3], which follow a pure on-line adaptation strategy, or active ones [4, 5], enabling learning only as a proactive reaction to a detected change in stationarity. In this paper, we follow this last learning strategy, though many results are general and can be suitably integrated in passive learning approaches as well.

Most change detection mechanisms have been proposed for numeric i.i.d. sequences and either rely on change point methods or change detection tests. Both change point methods and change detection tests are statistical tests; the former works off-line over a finite number of samples [6] while the latter employs a sequential analysis of incoming observations [7] to detect possible changes. These techniques were originally designed for univariate normal distributed variables, and only later developments extended the methodology to non-Gaussian distributions [8, 9] and multivariate datastreams [10, 11].

A somehow related field to change detection tests is one-class classification (e.g., see [12] and references therein). There, the idea is to model only the nominal state of a given system and detect non-nominal conditions (e.g., outliers, anomalies, or faults) by means of inference mechanisms. However, one-class classifiers typically process data batches with no specific presentation order, while change detection problems are sequential in nature.

The important role played, nowadays, by graphs as description of dynamic systems is boosting, also thanks to recent discoveries of theoretical frameworks for performing signal processing on graphs [13, 14] and for analysing temporal (complex) networks [15, 16, 17]. However, very few works address the problem of detecting changes in stationarity in streams (i.e., sequences) of graphs [18, 19, 20], and, to the best of our knowledge, none of them tackles the problem by considering a generic family of graphs (e.g., graphs with a varying number of vertices/edges, and arbitrary attributes on them). In our opinion, the reason behind such a lack of research in this direction lies in the difficulty of defining a sound statistical framework for generic graphs that, once solved, would permit to also detect changes in time variance in time-dependent graphs. In fact, statistics grounds on concepts like average and expectation, which are not standard concepts in graph domains. Fortunately, recent studies [21, 22] have provided some basic mathematical tools that allow us to move forward in this direction, hence, addressing the problem of detecting changes in stationarity in sequences of graphs.

A key problem in analysing generic graphs refers to assessing their dissimilarity, which is a well-known hard problem [23]. The literature proposes two main approaches for designing such a measure of dissimilarity [24, 25, 26]. In the first case, graphs are analysed in their original domain $\mathcal{G}$, whereas the second approach consists of mapping (either explicitly or implicitly) graphs to numeric vectors. A well-known family of algorithms used to assess dissimilarity between graphs relies on the Graph Edit Distance (GED) approach [27]. More specifically, GED algorithms count and weight the edit operations that are needed in order to make two input graphs equal. Differently, other techniques take advantage of kernel functions [28], spectral graph theory [29, 30], or assess graph similarity by looking for recurring motifs [31]. The computational complexity associated with the graph matching problem inspired researchers to develop heuristics and approximations, e.g., see [32, 33, 34] and references therein.

I-A Problem formulation

In this paper, we consider sequences of attributed graphs, i.e., directed or undirected labelled graphs characterized by a variable number of vertices and edges [21]. Attributed graphs associate vertices and edges with generic labels, e.g., scalars, vectors, categorical, and user-defined data structures. In addition, multiple attributes can be associated to the same vertex/edge, whenever requested by the application. By considering attributed graphs, we position ourselves on a very general framework covering most of application scenarios. However, generality requires a new operational framework, since all assumptions made in the literature to make the mathematics amenable, e.g., graphs with a fixed number of vertices and/or scalar attributes, cannot be accepted anymore. In order to cover all applications modellable through attributed graphs, we propose the following general problem formulation for change detection.

Given a generic pre-metric distance $d(\cdot,\cdot)$ on $\mathcal{G}$, we construct a $\sigma$-algebra $\mathcal{S}$ containing at least all open balls of $(\mathcal{G},d)$ and associate a generic probability measure $Q$ to $(\mathcal{G},\mathcal{S})$. The generated probability space $(\mathcal{G},\mathcal{S},Q)$ allows us to consider graphs as a realisation of a structured random variable $g$ on $(\mathcal{G},\mathcal{S},Q)$. Define $\mathcal{P}$ to be the process generating a generic graph $g_{t}\in\mathcal{G}$ at time $t$ according to a stationary probability distribution $Q$ (nominal distribution). We say that a change in stationarity occurs at (unknown) time $\tau$ when, from time $\tau$ on, $\mathcal{P}$ starts generating graphs according to a non-nominal distribution $\widetilde{Q}\neq Q$, i.e.,

[TABLE]

In this paper, we focus on persistent (abrupt) changes in stationarity affecting the population mean. However, our methodology is general and potentially can detect other types of change, including drifts and transient anomalies lasting for a reasonable lapse of time.

I-B Contribution and paper organization

A schematic description of the proposed methodology to design change detection tests for attributed graphs is shown in Figure 1, and consists of two steps: (i) mapping each graph $g_{t}$ to a numeric vector $y_{t}$ through a prototype-based embedding, and (ii) using a multivariate change detection test operating on the $y$-stream for detecting changes in stationarity.

The novelty content of the paper can be summarized as:

• •

A methodology to detect changes in stationarity in streams of attributed graphs. To the best of our knowledge, this is the first research contribution tackling change detection problems in streams of varying-size graphs with non-identified vertices and user-defined vertex/edge attributes;

• •

A method derived from the methodology to detect changes in stationarity in attributed graphs. We stress that the user can design his/her own change detection method by taking advantage of the proposed methodology;

• •

A set of theoretic results grounding the proposed methodology on a firm basis.

The proposed methodology is general and advances the few existing approaches for change detection in graph sequences mostly relying on the extraction and processing of topological features of fixed-size graphs, e.g., see [35].

It is worth emphasising that the proposed approach assumes neither one-to-one nor partial correspondence between vertices/edges across time steps (i.e., vertices do not need to be uniquely identified). This fact has important practical implications in several applications. As a very relevant example, we refer to the identification problem of neurons in extra-cellular brain recordings based on their activity [36]. In fact, each electrode usually records the activity of a neuron cluster, and single neurons need to be disentangled by a procedure called spike sorting. Hence, a precise identification of neurons (vertices) is virtually impossible in such an experimental setting, stressing the importance of methods that do not require one-to-one correspondence between vertices over time.

The remainder of this paper is structured as follows. Section II contextualizes our contribution and discusses related works. Section III-A presents the proposed methodology for change detection in generic streams of graphs. Theoretical results are sketched in Section III-B; related proofs are given in Appendix C. A specific implementation of the methodology is presented in Sections IV and related proofs in Appendix D. Section V shows experimental results conducted on datasets of attributed graphs. Finally, Section VI draws conclusions and offers future directions. Appendices A and B provide further technical details regarding the problem formulation.

II Related work

The relatively new field of temporal networks deals with graph-like structures that undergo events across time [16, 17]. Such events mostly realise in instantaneous or persistent contacts between pairs of vertices. With such structures one can study dynamics taking place on a network, like epidemic and information spreading, and/or dynamics of the network itself, i.e., structural changes affecting vertices and edges over time. Further relevant directions in temporal networks include understanding the (hidden) driving mechanisms and generative models [15].

The literature in statistical inference on time-varying graphs (or networks) is rather limited [16, 22], especially when dealing with attributed graphs and non-identified vertices. Among the many, anomaly detection in graphs emerged as a problem of particular relevance, as a consequence of the ever growing possibility to monitor and collect data coming from natural and man-made systems of various size. An overview of proposed approaches for anomaly and change detection on time-variant graphs is reported in [35, 37], where the authors distinguish the level of influence of a change. They identify changes affecting vertices and edges, or involving entire sub-networks of different size; this type of change usually concerns static networks, where the topology is often fixed. Other changes have a global influence, or might not be ascribed to specific vertices or edges.

We report that there are several applications in which the vertices are labelled in such a way that, from a time step to another, we are always able to create a partial one-to-one correspondence (identified vertices). This case arises, e.g., when the identity of vertices plays a crucial role and must be preserved over time. Here, we put ourself in the more general scenario where vertices are not necessarily one-to-one identifiable through time.

Within the anomaly detection context, only few works tackle the problem in a classical change detection framework. Among the already published works in detecting changes in stationarity, we mention Barnett et al. [19], whose paper deals with the problem of monitoring correlation networks by means of a change point method. In particular, at every time step $t$, the authors construct the covariance matrix computed from the signals up to time $t$ and the covariance matrix of the remaining data. As statistic for the change point model, they adopt the Frobenius norm between the covariance matrices. The authors evaluate their method on functional magnetic resonance imaging and stock returns. A different way to approach the problem consists in modeling the network-generating process within a probabilistic framework. Graphs with $N$ vertices and disjoint communities can be described by the degree corrected stochastic block model, where some parameters represent the tendency of single vertices to be connected and communities to interact. This model has been adopted by Wilson et al. for monitoring the U.S. Senate co-voting network [18]. As monitoring strategy, they consider the standard deviation of each community, and then apply exponential weighted moving average control chart. A further example of change point method for fixed-size graphs combines a generative hierarchical random graph model with a Bayesian hypothesis test [38].

III Change detection in a stream of graphs

The structure of the section is as follows. Section III-A describes the proposed methodology at a high level to ease the understanding. In Section III-B, we present theoretical results grounding our proposal; their proofs are given in the appendices.

III-A The proposed methodology

The methodology operates on an input sequence of attributed graphs $g_{1},g_{2},\dots,g_{t},\dots\,\in\mathcal{G}$ and, as sketched in Figure 1, it performs two steps: (i) Map (embed) input graphs to a vector domain $\mathcal{D}=\mathbb{R}^{M}$. Embedding is carried out by means of the dissimilarity representation $\zeta:\mathcal{G}\rightarrow\mathcal{D}$, which embeds a generic graph $g_{t}\in\mathcal{G}$ to a vector $y_{t}\in\mathcal{D}$. (ii) Once a multivariate i.i.d. vector stream $y_{1},y_{2},\dots,y_{t},\dots$ is formed, change detection is carried out by inspecting such a numerical sequence with the designer favourite method. The two phases are detailed in the sequel.

III-A1 Dissimilarity representation

The embedding of a generic graph $g\in\mathcal{G}$ is achieved by computing the dissimilarity between $g$ and the prototype graphs in $R=\{r_{1},\dots,r_{M}\}\subset\mathcal{G}$,

[TABLE]

The vector $y$ is referred to as the dissimilarity representation of $g$. Set $R$ has to be suitably chosen to induce informative embedding vectors. For a detailed discussion about dissimilarity representations and possible ways to select prototypes, we suggest [39].

In order to make the mathematics more amenable, here we assume $d(\cdot,\cdot)$ to be a metric distance; nevertheless, in practical applications, one can choose more general dissimilarity measures.

III-A2 Multivariate vector stream

At time step $t$ the process $\mathcal{P}$ generates graph $g_{t}$, and the map $\zeta(\cdot)$ embeds $g_{t}$ onto vector $y_{t}=\zeta(g_{t})\in\mathcal{D}$, inducing a multivariate stream $y_{1},y_{2},\dots,y_{t},\dots$ whose elements lie in $\mathcal{D}$. Figure 1 depicts the continuous embedding of graph process.

Under the nominal condition for process $\mathcal{P}$, graphs $\{g_{t}\}_{t<\tau}$ are i.i.d. and drawn from probability space $(\mathcal{G},\mathcal{S},Q)$. Consequently, also vectors $y_{t}\in\mathcal{D}$ are i.i.d.. We now define a second probability space $(\mathcal{D},\mathcal{B},F)$ associated with embedded vectors $y_{t}$; in particular, here we propose to consider for $\mathcal{B}$ the Borel’s $\sigma$-algebra generated by all open sets in $\mathcal{D}$. $F$ is the push forward probability function of $Q$ by means of $\zeta(\cdot)$, namely

[TABLE]

With such a choice of $F$, we demonstrate in Appendix A that $F$ is a probability measure on $(\mathcal{D},\mathcal{B})$.

III-A3 Change detection test

By observing the i.i.d. vector stream $y_{1},y_{2},\dots,y_{t},\dots$ over time we propose a multivariate change detection procedure to infer whether a change has occurred in the vector stream and, in turn, in the graph stream.

The change detection test is the statistical hypothesis test

[TABLE]

where [math] is the expected value during the nominal – stationary – regime and $\mathbb{E}[\cdot]$ is the expectation operator. Statistic $S_{t}$, which is applied to windows of the vector stream, is user-defined and is requested to increase when the process $\mathcal{P}$ becomes non stationary.

Often, the test comes with a threshold $h_{t}$ so that if

[TABLE]

and the estimated change time is

[TABLE]

Whenever the distribution of $S_{t}$ under hypothesis $H_{0}$ is available – or can be estimated – the threshold can be related to an user-defined significance level $\alpha_{t}$ so that $\alpha_{t}=\mathbb{P}(S_{t}>h_{t}|H_{0})$.

III-B Theoretical results

In this section, we show some theoretical results related to the methodology presented in Section III-A. In particular, we prove the following claims:

1. (C1)

Once a change is detected in the dissimilarity space according to a significance level $\alpha^{\prime}$, then a change occurs in probability with significance level $\alpha$ also in the graph domain; 2. (C2)

If a change occurs in the graph domain having set a significance level $\alpha$, then, with a significance level $\alpha^{\prime}$, a change occurs also in the dissimilarity space.

Figure 2 depicts the central idea behind the methodology. Through transformation $\zeta(\cdot)$, we map graphs to vectors. In the transformed space, we consider the expectation $\mathbb{E}[F]$ and the sample mean $\overline{\mathbf{\boldsymbol{y}}}$ associated with set $\mathbf{\boldsymbol{y}}=\{y_{t_{1}},\dots,y_{t_{n}}\}$ obtained by embedding the graph set $\mathbf{\boldsymbol{g}}=\{g_{t_{1}},\dots,g_{t_{n}}\}$, and design a hypothesis test of the form

[TABLE]

where $\delta$ is a positive threshold and hypothesis $H_{0}$ is associated with a nominal, change-free condition. In this paper, we relate the $y$-test of (4) to a correspondent test $d(\mu[Q],\mu[\mathbf{\boldsymbol{g}}])>\gamma$ operating in the graph domain, where $\mu[Q],\mu[\mathbf{\boldsymbol{g}}]$ are, respectively, the population and sample mean defined according to Fréchet [40]; for further details about Fréchet statistics refer to Appendix B.

Define $\alpha$ and $\alpha^{\prime}$ to be two significance levels, such that:

[TABLE]

In the sequel we relate the threshold $\gamma$ to $\gamma^{\prime}$, so that also the significance levels $\alpha$ and $\alpha^{\prime}$ are in turn related to each other.

In order to address our problem, we introduce mild assumptions to obtain closed-form expressions. Such assumptions are satisfied in most applications. More specifically,

• (A1)

We assume that the attributed graph space $(\mathcal{G},d)$ and the dissimilarity space $(\mathcal{D},d^{\prime})$ are metric spaces; in particular, $(\mathcal{G},d)$ is chosen as a graph alignment space [22] – i.e., a general metric space of attributed graphs – and $d^{\prime}(\cdot,\cdot)$ has to be induced by a norm.

• (A2)

We put ourselves in the conditions of [22] in order to take advantage of results therein; specifically, we assume that the Fréchet function $\mathcal{F}_{Q}(g)=\int_{\mathcal{G}}d^{2}(g,f)\,dQ(f)$ is finite for any $g\in\mathcal{G}$, and there exists a sufficiently asymmetric graph $f$ such that the support of $Q$ is contained in a cone around $f$. In this way, we are under the hypotheses of Theorems 4.1 and 4.2 of [22], which ensure the existence and uniqueness of the Fréchet population and sample mean in $\mathcal{G}$.

• (A3)

The embedding function $\zeta:(\mathcal{G},d)\rightarrow(\mathcal{D},d^{\prime})$ is bilipschitz, i.e., there exist two constants $c,C>0$, such that for any pair $g,f\in\mathcal{G}$

[TABLE]

Let $\Psi(\cdot)$ and $\Upsilon(\cdot)$ be the cumulative density functions (CDFs) of $d(\mu[\mathbf{\boldsymbol{g}}],\mu[Q])$ and $d^{\prime}(\overline{\mathbf{\boldsymbol{y}}},\mathbb{E}[F])$, respectively. Proposition 1 bounds the distribution $\Psi(\cdot)$ in terms of $\Upsilon(\cdot)$. This fact yields the possibility to derive significance levels $\alpha,\alpha^{\prime}$ and thresholds $\gamma,\gamma^{\prime}$ in Equation (5) that are related.

In order to prove Proposition 1, we make use of two auxiliary results, Lemmas 1 and 2. At first we need to comment that, although in general $\overline{\mathbf{\boldsymbol{y}}}\neq\zeta(\mu[\mathbf{\boldsymbol{g}}])$ and $\mathbb{E}[F]\neq\zeta(\mu[Q])$, differences are bounded in practice, as shown by Lemma 1.

Lemma 1**.**

Considering a set $\mathbf{\boldsymbol{g}}$ of i.i.d. random graphs and the associated embedded set $\mathbf{\boldsymbol{y}}$, there exists a constant $v_{2}$ such that

[TABLE]

Lemma 2 is used to derive bounds on the marginal distributions from bounds on the joint distributions that are useful in Proposition 1 to relate the threshold and the significance level in the graph space with the ones in the dissimilarity space.

Lemma 2**.**

Consider a random variable $x\in\mathcal{X}$ and two statistics $d_{1}(\cdot),d_{2}(\cdot):\mathcal{X}\rightarrow\mathbb{R}_{+}$ with associated CDFs $\Phi_{1}(\cdot)$, $\Phi_{2}(\cdot)$, respectively. If function $u:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ is increasing and bijective and $p$ is a constant in $[0,1]$, then

[TABLE]

Claims (C1) and (C2) now follow from the relation between the CDFs $\Psi(\cdot)$ and $\Upsilon(\cdot)$ proven by the subsequent Proposition 1. In particular, regarding Claim (C1) Proposition 1 provides a criterion for setting a specific threshold $\gamma^{\prime}$ for the test (4) for which we can state that $d(\mu[Q],\mu[\mathbf{\boldsymbol{g}}])$ is unexpectedly large with a significance level at most $\alpha$; similarly, we obtain Claim (C2).

Proposition 1**.**

Consider a sample $\mathbf{\boldsymbol{g}}$ of i.i.d. graphs. Under assumptions (A1)–(A3), if $\Psi(\cdot)$ and $\Upsilon(\cdot)$ are the CDFs of statistics $d(\mu[Q],\mu[\mathbf{\boldsymbol{g}}])$ and $d^{\prime}(\mathbb{E}[F],\overline{\mathbf{\boldsymbol{y}}})$, respectively, then for every $\delta>0$ there exist two values $b_{\delta}$ and $p_{\delta}$, depending on $\delta$ but independent of $\mathbf{\boldsymbol{g}}$, such that

[TABLE]

The proofs are given in Appendix C. Results of Proposition 1 allows us to state the major claim (C1): given a sample $\mathbf{\boldsymbol{g}}$, for any significance level $\alpha$ and threshold $\gamma$, as in (5), we can set up a $\gamma^{\prime}$ such that the confidence level of detecting a change in $\mathcal{D}$ is at least $p_{\delta}(1-\alpha)$. Specifically, for $\gamma^{\prime}=\frac{\gamma}{c}+b_{\delta}$ we have that

[TABLE]

Equation (8) states that if no change occurs in $\mathcal{G}$ with confidence level $1-\alpha$, then a change will not be detected in $\mathcal{D}$ according to threshold $\gamma^{\prime}$ at least with confidence level $p_{\delta}(1-\alpha)$. Indeed, this proves (C1) by contraposition. Similarly, Proposition 1 allows to prove Claim (C2). In fact, for any $\alpha^{\prime}$ and $\gamma^{\prime}$ as in (5), we can set $\gamma$ so that $\gamma^{\prime}=\frac{\gamma}{C}-b_{\delta}$ and obtain

[TABLE]

IV Implementations of the methodology

This section describes two examples showing how to implement the proposed methodology and derive actual change detection tests. In Section IV-A we present a specific method for generic families of graphs, whereas in Section IV-B we further specialize the methodological results to a special case considering graphs with identified vertices.

IV-A Change detection based on central limit theorem

Here, we consider specific techniques for prototype selection and change detection test. For the sake of clarity, we keep the subsection structure of Section III-A. Both the prototype selection and change detection test require a training phase. For this reason, the first observed graphs of the stream will serve as training set $T$, that we assume to be all drawn under nominal conditions.

IV-A1 Dissimilarity representation

Since the change detection method operates in the dissimilarity space, we need to define the embedding $\zeta(\cdot)$ that, at each time step, maps a generic graph $g_{t}$ to a vector $y_{t}$.

We comment that the embedding $\zeta(\cdot)$ is completely defined once the graph distance metric $d(\cdot,\cdot)$ and prototype set $R$ are chosen. Here, we adopt a metric GED as graph distance, since it meets Assumption (A1) of $(\mathcal{G},d)$ being a graph alignment space.

Many approaches have been proposed to deal with the prototype selection problem, e.g., see [41] and references therein. While the proposed methodology is general and one can choose any solution to this problem, here we adopt the k-Centres method [42]. The method selects prototypes so as to cover training data with balls of equal radius. Specifically, the algorithm operates as follows:

select $M$ random prototypes $R:=\{r_{1},\dots,r_{M}\}\subseteq T$; 2. 2.

for each $r_{m}\in R$, consider the set $C_{m}$ of all $g\in T$ such that $d(r_{m},g)=\min_{r\in R}d(r,g)$; 3. 3.

for $m=1,\dots,M$ update the prototypes $r_{m}$ with a graph ${c\in T}$ minimising $\max_{g\in C_{m}}d(c,g)$; 4. 4.

if the prototype set $R$ did not change from the previous iteration step then exit, otherwise go to step 2.

In order to improve the robustness of the k-Centres algorithm, we repeat steps 1–4 by randomizing initial conditions and select the final prototype set $R$ to be

[TABLE]

where $\{R^{(i)}\}$ is the collection of prototype sets found at each repetition.

IV-A2 Multivariate vector stream

Every time the process $\mathcal{P}$ generates a graph $g_{t}$, we embed it as $y_{t}=\zeta(g_{t})$ by using the prototype set $R$ identified with the k-Centres approach. This operation results in a multivariate vector stream $y_{1},y_{2},\dots,y_{t},\dots$ on which we apply the change detection test.

IV-A3 Change detection test

We consider here a variation of the cumulative sum (CUSUM) test [8] to design the change detection test. CUSUM is based on the cumulative sum chart [43], it has been proven to be asymptotically optimal [7] and allows for a simple graphical interpretation of results [44]. Here, the CUSUM is adapted to the multivariate case. However, we remind that, in principle, any change detection test can be used on the embedded sequence.

We batch the observed embedding vectors into non-overlapping samples $\mathbf{\boldsymbol{y}}_{w}:=\{y_{(w-1)n+1},\dots,y_{wn}\}$ of length $n$, where index $w$ represents the $w$-th data window. For each window, we compare the sample mean $\overline{\mathbf{\boldsymbol{y}}}_{0}$ estimated in the training set $T$ with that estimated in the $w$-th window, i.e., $\overline{\mathbf{\boldsymbol{y}}}_{w}$ and compute the discrepancy

[TABLE]

By assuming that $y_{1},y_{2},\dots,y_{t},\dots$ are i.i.d., and given sufficiently large $|T|$ and $n$, the central limit theorem grants that $\overline{\mathbf{\boldsymbol{y}}}_{0}$ and $\overline{\mathbf{\boldsymbol{y}}}_{w}$ are normally distributed. In particular, $\overline{\mathbf{\boldsymbol{y}}}_{0}$ and $\overline{\mathbf{\boldsymbol{y}}}_{w}$ share the same expectation $\mathbb{E}[F]$, while covariance matrices are $\frac{1}{|T|}\operatorname{Var}[F]$ and $\frac{1}{{n}}\operatorname{Var}[F]$, respectively.

As a specific choice of $d^{\prime}(\cdot,\cdot)$, we adopt the Mahalanobis’ distance, i.e., $d^{\prime}(\overline{\mathbf{\boldsymbol{y}}}_{0},\overline{\mathbf{\boldsymbol{y}}}_{w}):=d_{\Sigma}(\overline{\mathbf{\boldsymbol{y}}}_{0},\overline{\mathbf{\boldsymbol{y}}}_{w})$ where

[TABLE]

with matrix $\Sigma=\left(\frac{1}{|T|}+\frac{1}{n}\right)\operatorname{Var}[F]$, i.e., the covariance matrix of $\overline{\mathbf{\boldsymbol{y}}}_{0}-\overline{\mathbf{\boldsymbol{y}}}_{w}$. In our implementation, we consider as covariance matrix $\operatorname{Var}[F]$ the unbiased estimator $\frac{1}{|T|-1}\sum_{y\in T}(y-\overline{\mathbf{\boldsymbol{y}}}_{0})\cdot(y-\overline{\mathbf{\boldsymbol{y}}}_{0})^{\top}$.

For each stationary window $w$, the squared Mahalanobis’ distance $s_{w}^{2}$ is distributed as a $\chi^{2}_{M}$.

The final statistic $S_{w}$ inspired by the CUSUM test is defined as

[TABLE]

The difference $s_{w}-q$ increases with the increase of the discrepancy in (9) and positive values support the hypothesis that a change occurred, whereas negative values suggest that the system is still operating under nominal conditions. This behaviour resembles that of the original CUSUM increment associated with the log-likelihood ratio. In particular, the parameter $q$ can be used to tune the sensitivity of the change detection; in fact, if $q^{2}$ is the $\beta$-quantile $\chi^{2}_{M}(\beta)$, then $s_{w}-q$ produces a negative increment with probability $\beta$ and a positive one with probability $1-\beta$.

The last parameter to be defined is the threshold $h_{w}$ as requested in (3). Since the non-nominal distribution is unknown, a common criterion suggests to control the rate of false alarms, $\alpha$. In sequential detection tests, a related criterion requires a specific average run length under the null hypothesis (ARL0) [45]. ARL0 is connected to the false positive rate in the sense that setting a false alarm rate to $\alpha$ yields an ARL0 of $\alpha^{-1}$. Since we propose a sequential test, statistic $S_{w}$ depends on statistics at preceding times. As a consequence, since we wish to keep $\alpha$ fixed, we end up with a time-dependent threshold $h_{w}$. As done in [10, 11], we numerically determine thresholds through Monte Carlo sampling by simulating a large number of processes of repeated independent $\chi^{2}_{M}$ realisations. Threshold $h_{w}$ is then the $1-\alpha$ quantile of the estimated distribution of $S_{w}$.

We point out that, when setting a significance level for the random variable $S_{w}$, we are implicitly conditioning to the event $S_{i}\leq h_{i},\ \forall i<w$; in fact, when $S_{w}$ exceeds $h_{w}$, we raise an alarm and reconfigure the detection procedure.

IV-A4 Theoretical results

The choice of Mahalanobis’ distance ensures that almost all assumptions in Section III-B are met. In particular, the Mahalanobis’ distance meets the requirements of Assumption (A1). Then, the following Lemma 3 provides a lower bound of the form (6) in Assumption (A3); specifically, the lemma shows that, up to a positive factor, the distance between two graphs is larger then the one between the associated embedding vectors.

From these properties, we can apply Proposition 1 to state that Claim (C1) holds. Hence for any $\alpha$, we can set a specific threshold $\gamma^{\prime}$ yielding a confidence level at least $p_{\delta}(1-\alpha)$.

Lemma 3**.**

For any two graphs $g,f\in\mathcal{G}$, we have that

[TABLE]

where $\lambda_{M}$ is the smallest eigenvalue of $\Sigma$.

The proof is reported in Appendix D. Distance $d_{\Sigma}(\cdot,\cdot)$ is well-defined only when $\Sigma$ is positive definite, a condition implying that selected prototypes are not redundant.

IV-B A special case: graphs with identified vertices

Here we take into account the particular scenario where the attribute function $a(\cdot)$ of $g=(V,E,a)\in\mathcal{G}$ assigns numerical attributes in $[0,1]$ to vertices and edges of $g$; the vertex set $V$ is a subset of a predefined finite set $\mathcal{V}$, with $|\mathcal{V}|=N$. The peculiarity of this space $\mathcal{G}$ resides in the fact that any vertex permutation of a graph leads to a different graph. Many real-world applications fall in this setting, for instance correlation graphs obtained by signals coming from sensors or graphs generated by transportation networks.

We show an example of method for this setting, which complies with the methodology and satisfies Assumption (A3). This fact follows from the existence of an injective map $\omega$ from $\mathcal{G}$ to the $[0,1]^{N\times N}$ matrix set. Indeed, we represent each graph with its weighted adjacency matrix whose row/column indices univocally correspond to the vertices in $\mathcal{V}$. By endowing $\mathcal{G}$ with the Frobenius distance111In $\mathcal{G}$, $d_{F}(\cdot,\cdot)$ is a graph alignment distance, as formally shown in Appendix D-B. $d_{F}(g_{1},g_{2}):=\left\|\omega(g_{1})-\omega(g_{2})\right\|_{{F}}$, map $\omega:(\mathcal{G},d_{F})\rightarrow([0,1]^{N\times N},\left\|\cdot\right\|_{{F}})$ is an isometry.

Being the co-domain of $\omega$ an Euclidean space, we compute a matrix $X\in\mathbb{R}^{k\times M}$ whose columns $x_{i}$ constitute a $k$-dimensional vector configuration related to the prototype graphs; this is done via Classical Scaling [39], that is, $\left\|x_{i}-x_{j}\right\|_{{2}}=d(r_{i},r_{j})$ for all pairs of prototypes $r_{i},r_{j}\in R$. As usual, we consider the smallest possible $k$ that preserves the data structure as much as possible. Successively, for any dissimilarity vector $y\in\mathcal{D}$, we define $u:=XJ\,y^{*2}$ to be a linear transformation of $y^{*2}$, obtained by squaring all the components of $y$; the matrix $J$ is the centering matrix $I-\frac{1}{M}\mathbf{\boldsymbol{1}}\mathbf{\boldsymbol{1}}^{\top}$. We apply the same procedure of Section IV, considering $u$ instead of $y$: matrix $\Sigma$ is derived from the non-singular covariance matrix222$\Sigma$ is non-singular, since there are no isolated points in $\mathcal{G}$ and as a consequence of the selection of $k$., and the statistic $s_{w}$ is the Mahalanobis distance $d_{\Sigma}\left(\overline{\mathbf{\boldsymbol{u}}}_{0},\overline{\mathbf{\boldsymbol{u}}}_{w}\right)$.

Considering the space $(\mathcal{G},d_{F})$ and the above transform, we claim that the following lemma holds. Lemma 4 proves the fulfilment of Assumption (A3).

Lemma 4**.**

For any positive definite matrix $\Sigma\in\mathbb{R}^{k\times k}$,

[TABLE]

where $c=\sqrt{\frac{\lambda_{k}(\Sigma)}{4\lambda_{1}(XX^{\top})}}$, $C=\sqrt{\frac{\lambda_{1}(\Sigma)}{4\lambda_{k}(XX^{\top})}}$. $\lambda_{i}(\cdot)$ is the $i$-th eigenvalue in descending order of magnitude.

The proof is reported in Appendix D.

V Experiments

The proposed methodology can operate on very general families of graphs. Besides the theoretical foundations around which the paper is built, we provide some experimental results showing the effectiveness of what proposed in real change detection problems in streams of graphs. In particular, we consider the method introduced in Section IV-A as an instance of our methodology. Source code for replicating the experiments is available in [46].

V-A Experiment description

V-A1 Data

The experimental evaluation is performed on the well-known IAM benchmarking databases [47]. The IAM datasets contain attributed graphs representing drawings and biochemical compounds. Here we consider the Letters, Mutagenicity, and AIDS datasets.

The Letters dataset contains 2-dimensional geometric graphs. As such, each vertex of the graphs is characterized by a real vector representing its location on a plane. The edges define lines such that the graphical planar representation of the graphs resembles a latin-script letter. The dataset is composed of 15 classes (one for each letter333The IAM letters database [47] considers only non-curved letters; hence, e.g., letters A and E are considered, whereas B and C are excluded.) containing 150 instances each.

Conversely, the Mutagenicity and AIDS datasets contain biological molecules. Molecules are represented as graphs by considering each atom as a vertex and each chemical link as an edge. Each vertex is attributed with a chemical symbol, whereas the edges are labelled according to the valence of the link. Both datasets contain two classes of graphs: mutagenic/non-mutagenic for Mutagenicity and active/inactive for AIDS. The two datasets are imbalanced in terms of size of each classes, in particular Mutagenicity has 2401 mutagenic and 1963 non-mutagenic molecules; AIDS contains 400 active and 1600 inactive molecules.

We considered these datasets because they contain different types of graphs with variegated attributes (numerical and categorical). We refer the reader to [47] and references therein for a more in-depth discussion about the datasets.

V-A2 Simulating the generating process $\mathcal{P}$

For each experiment in Table I, we consider two collections of graphs containing all possible observations in the nominal and non-nominal regimes, respectively. Each collection is composed by graphs present in one or more predefined classes of the dataset under investigation. The collections have to be different, but they do not need to be disjoint; as such, some graphs can belong to both collections.

Next, we simulate the process $\mathcal{P}$ by bootstrapping graphs from the first collection up to the predefined change time $\tau$; this is the nominal regime. After $\tau$, we bootstrap objects from the second collection hence modelling a change.

Regarding molecular datasets, we considered two distinct experiments. For the Mutagenicity (AIDS) experiment, we set the non-mutagenic (inactive) class as nominal collection and mutagenic (active) class as non-nominal one. On the other side, for the Letter dataset we design four different experiments depending on which classes will populate the collections. Table I reports the settings of all the experiments and next section describes the relevant parameters.

V-A3 Parameters setting

For all experiments, the offset $q$ is set to the third quartile of the $\chi^{2}(M)$ distribution. The time-dependent threshold $h_{w}$ has been numerically estimated by Monte Carlo sampling. We drew one million processes of i.i.d. random variables $s_{w}$ by taking the square root of i.i.d. random variables distributed as $\chi^{2}(M)$. For each obtained process (stream), we computed the sequence of cumulative sums $S_{w}$ like in (10) and estimated the threshold $h_{w}$ as the quantile of order $\alpha_{w}=$1/ARL0, with ARL0 $=200$ (windows).

We divided the training set $T$ in two disjoint subsets, $T_{c},T_{p}$, used during the prototype selection and change detection learning phases, respectively. We set $|T_{p}|=300$ and $|T_{c}|=1000$, afterwards we generated a stream of graphs containing $20\cdot n\cdot$ARL0 observations associated with the operational phase. The change is forced at time $\tau=12\cdot n\cdot\text{ARL0}$. As for the distance $d(\cdot,\cdot)$, we considered the bipartite GED implemented in [48], where we selected the Volgenant and Jonker assignment algorithm. The other GED parameters are set according to the type of graphs under analysis, i.e., for geometric graphs we consider the Euclidean distance between numerical attributes, and a binary (0-1 distance) for categorical attributes. The k-Centres procedure is repeated 20 times.

We believe that the selected parameter settings are reasonable. Nevertheless, a proper investigation of their impact w.r.t. performance metrics is performed in a companion paper [49], hence it is outside the focus of the present one.

V-A4 Figures of merit

We assess the performance of the proposed methodology by means of the figures of merit here described. Such measurements are obtained by replicating each experiment one hundred times; we report the average of the observed measures with their estimated 95% confidence interval (95CI) or standard deviation (std).

First of all, we consider the observed ARL0 introduced in Section III-A3 and the delay of detection (DoD). Both of them are computed as the average time lapses between consecutive alarms, but limiting to those alarms raised before and after time $\tau$, respectively. From them we estimate the rate of detected changes (DCR), by assessing the rate of simulations in which the DoD is less then the observed ARL0. Finally, we consider also the estimated rate of false anomalies within 1000 samples (FA1000). This is computed as the ratio between the count of raised false alarms and the total number of thousands of time steps under the nominal condition.

We point out that the measures ARL0, DoD, and FA1000 are computed with the window as unitary time step.

V-A5 Baseline methods

As previously mentioned, state-of-the-art change detection methods for graph streams usually assume a given topology with a fixed number of vertices and/or simple numeric features for vertices/edges. As reported in [35], considering a variable topology, a common methodology for anomaly detection on graphs consists of extracting some topological features at each time step and then applying a more general anomaly detector on the resulting numeric sequence. Accordingly, in addition to the method proposed in Section IV, we consider two baseline methods for comparison. More precisely, we considered two topological features: the density of edges $\phi_{1}(g)=\frac{|E|}{|V|(|V|-1)}$ and the spectral gap $\phi_{2}(g)=|\lambda_{1}(L(g))|-|\lambda_{2}(L(g))|$ of the Laplacian matrix $L(g)$ [50]. The particular choices of $\phi_{1}$ and $\phi_{2}$ can be justified by considering that both features are suitable for describing graphs with a variable number of vertices and edges. We implemented two CUSUM-like change detection tests as in Section IV where, for $i=1,2$, the statistic $s_{t}$ is now given by $\left|\phi_{i}(g_{t})-\mathbb{E}[\phi_{i}(g_{t})]\right|$, and $\mathbb{E}[\phi_{i}(g_{t})]$ is numerically estimated in the training phase.

In addition, we consider a further baseline implemented as a degenerate case of our method by selecting only $M=1$ prototype and window size $n=25$. This last baseline is introduced to show that the strength of our methodology resides also in the embedding procedure, and not only in the graph distance $d(\cdot,\cdot)$.

V-B Results on IAM graph database

For the sake of readability, we show results for our method and baselines in two different tables, that is, Table II and III, respectively.

In all experiments shown in Table I, there is a parameter setting achieving a detection rate (DCR) statistically larger than $0.99$. Indeed, in Table II the 95% confidence interval (95CI) of the DCR is above $0.99$ (see symbol $*\;$in the table). Looking at Table II more in detail, we notice that both the window size and the number of prototypes yield higher DCR. In particular, this can be seen in L-O experiments, where all DCR estimates have disjoint 95CIs (i.e., differences are statistically significant). The same phenomenon appears also for L-D5, L-S, and MUT, e.g., see symbols $\circ\;$and $\bullet\;$in Table II. As far as other figures of merit are concerned, we do not observe statistical evidence of any trend. Still, with the exception of a few cases in MUT, all 95CIs related to ARL0 contain the target value ARL0=200; hence, we may say that the threshold estimation described in Section V-A3 completed as expected.

Here, we limit the analysis to the proposed parameter settings for $n$ or $M$, since we already reach the highest possible DCR, achieving one hundred detections out of one hundred. We believe that, by increasing the window size $n$ the false alarms will decrease, as our method relies on the central limit theorem. Concerning the number $M$ of prototypes, we point out that, in the current implementation of the methodology, the number of parameters to be estimated scales as $M^{2}$; accordingly, we need to increase the number of samples.

The second experimental analysis addresses the performance assessment of the three baselines of Section V-A5 on the experiments of Table I. The results reported in Table III show that, in some cases, the considered baselines achieve sound performance, which is comparable to the one shown in Table II. Comparing Table III with Table II, by intersecting the 95% confidence intervals, we notice that there is always one of the proposed methods which attains DCR that is statistically equivalent or better than the baselines, see symbol $*\;$in Table II. In particular, the proposed method performs significantly better than the baselines on the MUT dataset.

Finally, Table III shows that the method based on edge density $\phi_{1}$ is significantly more accurate in terms of DCR than the one based on spectral gap $\phi_{2}$ in almost all experiments; confidence intervals at level 95% do not intersect. The first method performed better than the degenerate one (M1) with only one prototype: see L-O and L-S. Conversely, M1 outperforms $\phi_{1}$ on L-D5.

VI Conclusions

In this paper, we proposed a methodology for detecting changes in stationarity in streams of graphs. Our approach allows to handle general families of attributed graphs and is not limited to graphs with fixed number of vertices or graphs with (partial) one-to-one correspondence between vertices at different time steps (uniquely identified vertices). The proposed methodology consists of an embedding step, which maps input graphs onto numeric vectors, bringing the change detection problem back to a more manageable setting.

We designed and tested a specific method as an instance of the proposed methodology. The embedding has been implemented here as a dissimilarity space representation, hence relying on a suitable set of prototype graphs, which, in our case, provide also a characterization of the nominal condition, and a dissimilarity measure between graphs, here implemented as a graph edit distance. The method then computes the Mahalanobis’ distance between the mean values in two windows of embedded graphs and adopts a CUSUM-like procedure to identify changes in stationarity.

We provided theoretical results proving that, under suitable assumptions, the methodology is able to relate the significance level with which we detect changes in the dissimilarity space with a significance level that changes also occurred in the graph domain; also the vice versa has been proven. We also showed that our methodology can handle more basic, yet relevant scenarios with uniquely identified vertices.

Finally, we performed experiments on IAM graph datasets showing that the methodology can be applied both to geometric graphs (2-dimensional drawings) and graphs with categorical attributes (molecules), as instances of possible data encountered in real-world applications. Results show that the proposed method attains at least comparable (and often better) results w.r.t. other change detectors for graph streams.

In conclusion, we believe that the proposed methodology opens the way to designing sound change detection methods for sequences of attributed graphs with possibly time-varying topology and non-identified vertices. In future studies, we plan to work on real-world applications and focus on the automatic optimization of relevant parameters affecting the performance.

Appendix A Correspondence between probability spaces

Consider the measurable space $(\mathcal{D},\mathcal{B})$ introduced in Section III-A2, and the probability space $(\mathcal{G},\mathcal{S},Q)$ of Section I-A.

Let us define the preimage function $\zeta^{-1}(B):=\{g\in\mathcal{G}:\zeta(g)\in B\}$, with $\zeta^{-1}(\varnothing)=\varnothing$ and consider the smallest $\sigma$-algebra $\mathcal{S}$ containing all open balls444A ball $O(\rho,g)$ is defined as a set $\{f\in\mathcal{G}:d(g,f)<\rho\}$ of all graphs $f\in\mathcal{G}$ having distance $d(f,g)$ w.r.t. a reference graph $g\in\mathcal{G}$ smaller than the radius $\rho>0$. $O(\rho,g)$ and preimage sets w.r.t. any $B\in\mathcal{B}$

[TABLE]

and a generic probability density function $Q:\mathcal{S}\rightarrow[0,1]$ on $\mathcal{S}$. Then, we can define the function $F:\mathcal{B}\rightarrow[0,1]$ as in Equation (2).

The triple $(\mathcal{D},\mathcal{B},F)$ is a probability space. The following three properties provide a proof that $F$ is a probability measure on $(\mathcal{D},\mathcal{B})$:

• •

$F(\mathcal{D})=Q(\zeta^{-1}(\mathcal{D}))=Q(\mathcal{G})=1$,

• •

$F(\varnothing)=Q(\varnothing)=0$, and

• •

for any countable collection $\{B_{i}\}\in\mathcal{B}$ of pairwise disjoint sets, $\zeta^{-1}(\cup B_{i})=\cup\zeta^{-1}(B_{i})$, hence the sets $\zeta^{-1}(B_{i})$ are pairwise disjoint and $F(\cup B_{i})=\cup Q(\zeta^{-1}(B_{i}))=\cup F(B_{i})$.

Notice also that, indicating with $\operatorname{Im}(\zeta)$ the image set $\{\zeta(g)|g\in\mathcal{G}\}\subseteq\mathcal{D}$ of $\zeta(\cdot)$, we have $F(B)=0$ for any $B\in\mathcal{B}$ such that $B\cap\operatorname{Im}(\zeta)=\varnothing$.

Appendix B Fréchet mean

Given a probability space $(\mathcal{X},\mathcal{S},P)$ defined on a metric space $(\mathcal{X},d)$, we consider a random sample $\mathbf{\boldsymbol{x}}=\{x_{t_{1}},\dots,x_{t_{n}}\}$.

B-A Definition of Fréchet mean and variation

For any object $x\in\mathcal{X}$, let us define the functions $\mathcal{F}_{\mathbf{\boldsymbol{x}}}(x):=\frac{1}{n}\sum_{t=t_{1}}^{t_{n}}d(x,x_{t})^{2}$ and $\mathcal{F}_{P}(x):=\int_{\mathcal{G}}d(x,x^{\prime})^{2}\,dP(x^{\prime})$. A Fréchet sample (population) mean is any object $x$ attaining the minimum of the function $\mathcal{F}_{\mathbf{\boldsymbol{x}}}(\cdot)$ ($\mathcal{F}_{P}(\cdot)$). We point out that the minimum might not exist in $\mathcal{X}$ and, if it does, it can be attained at multiple objects. Whenever the minimum exists and is unique, we refer to it as $\mu[\mathbf{\boldsymbol{x}}]$ ($\mu[P]$). In addition, we define Fréchet sample (population) variation as the infimum $\operatorname{V_{f}}[\mathbf{\boldsymbol{x}}]:=\inf_{x}\mathcal{F}_{\mathbf{\boldsymbol{x}}}(x)$ ($\operatorname{V_{f}}[P]:=\inf_{x}\mathcal{F}_{P}(x)$).

B-B Fréchet mean in Euclidean spaces

In the case of a set $\mathcal{X}\subseteq\mathbb{R}^{d}$ and distance $d(\cdot,\cdot)=\left\|\cdot-\cdot\right\|_{{2}}$, the space is Euclidean. Firstly, we show that $\mu[P]=\mathbb{E}[P]$ and $\mu[\mathbf{\boldsymbol{x}}]=\overline{\mathbf{\boldsymbol{x}}}$, then we show that

[TABLE]

1. The following equality holds

[TABLE]

as we will show below, in Part 3. This result proves that the minimum is attained in the expectation $z=\mathbb{E}[P]$.

2. Similarly, in the ‘sample’ case

[TABLE]

proving that $\mu[\mathbf{\boldsymbol{x}}]=\overline{\mathbf{\boldsymbol{x}}}$.

3. The results of previous Parts 1 and 2 are derived from the following three equalities: for any $z\in\mathbb{R}^{d}$

[TABLE]

4. Let us move on proving the second result. Notice that

[TABLE]

5. Then,

[TABLE]

Now, thanks to the independence of the observations

[TABLE]

which proves the thesis.

6. Finally, in this part, we prove a result that holds in general. Notice that $\frac{1}{n}\sum_{t}d(x_{t},\mu[P])^{2}\geq\operatorname{V_{f}}[\mathbf{\boldsymbol{x}}]$. Then, since $\mathbf{\boldsymbol{x}}$ are i.i.d. so also are the $\{d(x_{t},\mu[P])^{2}\}$. Thanks to the monotonicity of the expectation, we have

[TABLE]

Appendix C Proofs of Section III-A

C-A Proof of Lemma 1

Notice that the means $\mu[\mathbf{\boldsymbol{y}}]=\overline{\mathbf{\boldsymbol{y}}}$ and $\mu[F]=\mathbb{E}[F]$ are computed w.r.t. the Euclidean metric, the and $d^{\prime}(\cdot,\cdot)$ is deployed as statistic for the change detection test.

1. From Equation (12), for any $z\in\mathcal{D}$ we have

[TABLE]

2. We provide a second inequality that will be useful later. Given three graphs $g,f,r\in\mathcal{G}$, we have $d(g,f)\geq d(g,r)-d(r,f)$ from the triangular inequality (Assumption (A1)); since this holds for any prototype $r\in R$, it proves that

[TABLE]

3. Exploiting the inequality (14) in Part 2, and taking $z=\zeta(\mu[Q])$, we obtain

[TABLE]

In Appendix A we observed that $F(\mathcal{D}\setminus\operatorname{Im}(\zeta))=0$, then

[TABLE]

Eventually, combining with Part 1, we obtain the first part of the thesis for $v_{0}:=M\,\operatorname{V_{f}}[Q]-\operatorname{V_{f}}[F]\geq 0$, in fact

[TABLE]

4. Similarly, we have

[TABLE]

and, for $v_{\mathbf{\boldsymbol{g}}}:=M\,\operatorname{V_{f}}[\mathbf{\boldsymbol{g}}]-\operatorname{V_{f}}[\mathbf{\boldsymbol{y}}]$,

[TABLE]

By exploiting (11) and (13), and the monotonicity of the expected value

[TABLE]

As final remarks, we point out that here we left the dependency from $n$ on purpose, but an independent bound can be easily found, e.g., $v_{2}=M\,\operatorname{V_{f}}[Q]-\tfrac{1}{2}\operatorname{V_{f}}[F]$.

5. The random variable $\left\|\overline{\mathbf{\boldsymbol{y}}}-\zeta(\mu[\mathbf{\boldsymbol{g}}])\right\|_{{2}}^{2}$ is non-negative. As such we can apply Theorem 1, Section 5.4 in [51] (sometimes called Markov’s inequality) and obtain, for each $\delta>0$,

[TABLE]

Now, from the previous Part 4 we conclude that

[TABLE]

C-B Proof of the Lemma 2

For convenience, let us define the quantities

[TABLE]

1. By the law of total probability, and for any $\gamma\geq 0$,

[TABLE]

Lower-bounding the second addendum with zero and by hypothesis,

[TABLE]

2. Notice that $\mathbb{P}(d_{1}(x)\leq u(\gamma)|d_{2}(x)=\gamma,x\in A)=1$ for all $\gamma\geq 0$, thanks to the fact that $x\in A$; hence, we have $\rho(-|-)=1$. Applying again the law of total probabilities,

[TABLE]

Combining with the above Part 1, we have

[TABLE]

3. A final remark is that Lemma 2 proves also that, if $\ell:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ is bijective and increasing providing a lower-bound, then for any $\gamma>0$ and $p>0$

[TABLE]

In fact $d_{1}\geq\ell(d_{2})$ if and only if $d_{2}\leq\ell^{-1}(d_{1})$, therefore, we obtain the result by applying the lemma to $u(\cdot)=\ell^{-1}(\cdot)$ and inverting the roles of $d_{1}$ and $d_{2}$.

C-C Proof of Proposition 1

By Assumption (A2), we can consider the following function of $\mathbf{\boldsymbol{g}}$ and their respective CDFs

[TABLE]

recalling that $\mathbf{\boldsymbol{y}}=(\dots,\zeta(g_{i}),\dots)^{\top}$.

1. From triangular inequality (Assumption (A1)), we have

[TABLE]

By the equivalence of the 1- and 2-norm in $\mathbb{R}^{2}$,

[TABLE]

Now, since $d^{\prime}(\cdot,\cdot)$ is induced by a norm (Assumption (A1)), we can deploy the equivalence of any pair of norms in $\mathbb{R}^{M}$. Letting $m$ be a constant for relating the distance $d^{\prime}(\cdot,\cdot)$ with the Euclidean one, and applying Lemma 1, we obtain

[TABLE]

By exploiting again Lemma 1, for any $\delta>0$, we have

[TABLE]

So, with $b(\delta)=\sqrt{\frac{m^{2}}{2}(v_{2}+\delta)}$ and $p(\delta)=1-\frac{v_{2}}{\delta}$, we have the following estimates

[TABLE]

2. With Assumption (A3), we can exploit Lemma 2 by setting $u(x)=C\,x$, $d_{1}(\mathbf{\boldsymbol{g}})=d(\mathbf{\boldsymbol{g}})$ and $d_{2}(\mathbf{\boldsymbol{g}})=d^{\prime}(\mathbf{\boldsymbol{g}})$. Here $\mathbb{P}(d_{1}\leq u(d_{2}))=1$, hence we obtain that $\Psi(\gamma)\geq\Phi\left(\frac{\gamma}{C}\right)$. Accordingly, employing also (15), we have

[TABLE]

3. Consider this time (16) with $d_{1}(\mathbf{\boldsymbol{g}})=d^{\prime}_{0}(\mathbf{\boldsymbol{g}})$ and $d_{2}(\mathbf{\boldsymbol{g}})=d^{\prime}(\mathbf{\boldsymbol{g}})$. By applying Lemma 2 twice, with $u(x)=x+b(\delta)$ and $\ell(x)=x-b(\delta)$, we obtain

[TABLE]

4. Combining previous parts, Parts 2 and 3, we obtain

[TABLE]

Appendix D Proofs of Sections IV and IV-B

D-A Proof of Lemma 3

Recall that $\left\|x\right\|_{{2}}^{2}\geq\lambda_{M}x^{\top}\Sigma^{-1}x$, where $\lambda_{M}$ is the smallest eigenvalue of the covariance matrix $\Sigma$. By the triangular inequality we have (14), which leads to $c=\sqrt{\frac{\lambda_{M}}{M}}$.

D-B An instance of graph alignment space

A graph alignment space is a pair $(\mathcal{G},d)$. According to what defined in Section IV-B, we explicit an attribute kernel determining a graph alignment distance $d$ acting like $d_{F}$.

We consider $\mathcal{A}=\mathcal{V}\times[0,1]$ as attribute set. Following [21, Ex.3.4], for any pair $(x,v)\in\mathcal{A}$ we define the feature map $\mathbf{\boldsymbol{\Phi}}:(x,v)\mapsto(x,e_{v})$, where $e_{v}\in\mathbb{R}^{N}$ is zero everywhere except for the position $v$ where it assumes the fixed value $\nu$. The attribute kernel $k:\mathcal{A}\times\mathcal{A}\rightarrow\mathbb{R}$ is then defined accordingly. We consider a vertex disabled if there are no incoming or outgoing edges, but it has always a label of the form $(0,v)$.

It might happen that the “cost” of matching vertices is overcome by the improvement of a better topology alignment. This is something which is not contemplated in the distance $d_{F}$. For avoiding this, it suffices to set the constant $\nu$ sufficiently large, so that the trivial identity alignment is always the optimal alignment (i.e., vertices are uniquely identified). Setting $\nu=N^{2}$ ensures this behaviour and yields an alignment distance $d(\cdot,\cdot)$ acting like $d_{F}(\cdot,\cdot)$.

D-C Proof of Lemma 4

1. For every prototye set $R$, the dissimilarity matrix $D(R,R):=[d_{F}(r_{i},r_{j})]_{i,j}$ is Euclidean, since it has been generated by using the Euclidean distance. In general, the dimension of $\mathbb{R}^{k}$ is $k=N^{2}$, although it can be lower. We remind that $N$ is the number of vertices assumed for the input graphs. We use the notation $x_{1},\dots,x_{M}\in\mathbb{R}^{k}$ to denote the prototypes $r_{1},\dots,r_{M}\in R$ w.r.t. the classical scaling process and $X$ is obtained by stacking $[x_{1}|\dots|x_{M}]$ as columns (see [39, Sec. 3.5.1]). The minimal number of prototypes for representing $\mathbb{R}^{k}$ is $M=k+1$, and this holds true only when the matrix $X$ is full rank. For obtaining the embedding $z$ of a generic graph $g\in\mathcal{G}\setminus R$, we compute $y^{*2}:=[d_{F}(g,r_{1})^{2},\dots,d_{F}(g,r_{M})^{2}]^{\top}$ and we can solve the following linear system with respect to $z$

[TABLE]

where $J=I-\frac{1}{n}\mathbf{\boldsymbol{1}}\mathbf{\boldsymbol{1}}^{\top}$ is the centering matrix, $I$ is the identity matrix, and $D^{*2}$ is the component-wise square of $D(R,R)$ (see [39, Cor. 3.7]). The solution is unique, provided that the rank of $X$ is $k$.

2. The differences $\delta z=z_{1}-z_{2}$ and $\delta y^{*2}=y^{*2}_{1}-y^{*2}_{2}$ are related by the equation $2\,XX^{\top}\delta z=-XJ\,\delta y^{*2}$.

Being $\left\|a\right\|_{{A}}^{2}=a^{\top}\,A\,a$ we have that $4{\left\|\delta z\right\|_{{(XX^{\top})^{2}}}}=4\left\|XX^{\top}z\right\|_{{I}}={\left\|XJ\,\delta y^{*2}\right\|_{{I}}}$ and

[TABLE]

We want to show that $c\leq{\left\|\delta z\right\|_{{I}}}/{\left\|XJ\,\delta y^{*2}\right\|_{{\Sigma^{-1}}}}\leq C$, in fact $\left\|\delta z\right\|_{{I}}=\left\|\delta z\right\|_{{2}}=d_{F}(g_{1},g_{2})$, and $\left\|XJ\,\delta y^{*2}\right\|_{{\Sigma^{-1}}}=d_{\Sigma}(XJ\,y_{1}^{*2},XJ\,y_{2}^{*2})$.

3. Applying [52, Theo.1], we can bound ${\left\|a\right\|_{{A}}}/{\left\|a\right\|_{{B}}}$ in terms of the values $\beta$ for which there exists non-null $q\in\mathbb{R}^{M}$ such that $(A-\beta B)\,q=\mathbf{\boldsymbol{0}}$. For the pair $(I,(XX^{\top})^{2})$ the values $\beta$ are provided by the square of the inverse eigenvalues $\lambda_{1}(XX^{\top}),\dots,\lambda_{k}(XX^{\top})$ of $XX^{\top}$. The eigenvalues are reported in descending order. For what concerns $(I,\Sigma^{-1})$, instead, we see that the values $\beta$ corresponds to the eigenvalues $\lambda_{1}(\Sigma),\dots,\lambda_{k}(\Sigma)$ of $\Sigma$. In the end, the positive constants $c,C$ are

[TABLE]

The non-singularity of $XX^{\top}$ makes the above fractions feasible.

A final comment on how to exploit Lemma 4 in Proposition 1 is reported in Appendix D-D.

D-D Proposition 1 for the special case of Section IV-B

When considering identified vertices, while Lemma 2 is still valid, Lemma 1 has to be adapted. Despite the new embedding $\zeta_{0}(g)=XJ\zeta(g)^{*2}$ is slightly different from the original one in (1), $\zeta_{0}(\cdot)$ and $\mathcal{D}_{0}$ can be treated in a similar way. In particular, a result equivalent to Lemma 1 can be proved for $\zeta_{0}(\cdot)$ by solely adapting Part 2 of the proof shown in Appendix C-A,

[TABLE]

Paying attention in substituting the map $\zeta(\cdot)$ with $\zeta_{0}(\cdot)$, the rest of the proof holds with $M=\lambda_{1}(XX^{\top})^{2}$ thus obtaining the Claims (C1) and (C2), with a final result similar to the one of Proposition 1.

Acknowledgements

This research is funded by the Swiss National Science Foundation project 200021_172671: “ALPSFORT: A Learning graPh-baSed framework FOr cybeR-physical sysTems”.