Time Series Prediction for Graphs in Kernel and Dissimilarity Spaces

TL;DR

This paper introduces a novel approach for predicting the future state of entire graphs over time using kernel and dissimilarity-based regression techniques, enabling better understanding of graph dynamics in various fields.

Contribution

It proposes framing graph evolution prediction as a regression problem with kernel and dissimilarity methods, including strategies to accelerate Gaussian process regression.

Findings

Kernel regression captures theoretical graph dynamics effectively.

Gaussian process regression improves accuracy on real-world data.

Speed-up strategies make Gaussian processes computationally feasible.

Abstract

Graph models are relevant in many fields, such as distributed computing, intelligent tutoring systems or social network analysis. In many cases, such models need to take changes in the graph structure into account, i.e. a varying number of nodes or edges. Predicting such changes within graphs can be expected to yield important insight with respect to the underlying dynamics, e.g. with respect to user behaviour. However, predictive techniques in the past have almost exclusively focused on single edges or nodes. In this contribution, we attempt to predict the future state of a graph as a whole. We propose to phrase time series prediction as a regression problem and apply dissimilarity- or kernel-based regression techniques, such as 1-nearest neighbor, kernel regression and Gaussian process regression, which can be applied to graphs via graph kernels. The output of the regression is a point embedded in a pseudo-Euclidean space, which can be analyzed using subsequent dissimilarity- or kernel-based processing methods. We discuss strategies to speed up Gaussian Processes regression from cubic to linear time and evaluate our approach on two well-established theoretical models of graph evolution as well as two real data sets from the domain of intelligent tutoring systems. We find that simple regression methods, such as kernel regression, are sufficient to capture the dynamics in the theoretical models, but that Gaussian process regression significantly improves the prediction error for real-world data.

Tables

Table 1. Table 1: The mean RMSE and runtime across cross validation trials for both theoretical data sets (x-axis) and all methods (y-axis). The standard deviation is shown in brackets. Runtime entries with 0.000 0.000 0.000 had a shorter runtime (and standard deviation) than 10 − 3 superscript 10 3 10^{-3} milliseconds. The best (lowest) value in each column is highlighted by bold print.

	Barabási-Albert		Game of Life
method	RMSE	runtime [ms]	RMSE	runtime [ms]
identity	$0.137$ ($0.005$)	$0.000$ ($0.000$)	$1.199$ ($0.455$)	$0.000$ ($0.000$)
1-nn	$0.073$ ($0.034$)	$0.111$ ($0.017$)	$1.191$ ($0.442$)	$0.112$ ($0.025$)
KR	$0.095$ ($0.039$)	$0.122$ ($0.016$)	$0.986$ ($0.398$)	$0.120$ ($0.040$)
GPR	$0.064$ ($0.028$)	$0.148$ ($0.022$)	$0.965$ ($0.442$)	$0.127$ ($0.026$)
rBCM	$0.062$ ($0.015$)	$0.312$ ($0.083$)	$0.967$ ($0.461$)	$0.267$ ($0.077$)

Table 2. Table 2: The mean RMSE and runtime across cross validation trials for both Java data sets (x-axis) and all methods (y-axis). The standard deviation is shown in brackets. Runtime entries with 0.000 0.000 0.000 had a shorter runtime (and standard deviation) than 10 − 3 superscript 10 3 10^{-3} seconds. The best (lowest) value in each column is highlighted by bold print.

	MiniPalindrome		Sorting
method	RMSE	runtime [s]	RMSE	runtime [s]
identity	$0.295$ ($0.036$)	$0.000$ ($0.000$)	$0.391$ ($0.029$)	$0.000$ ($0.000$)
1-NN	$0.076$ ($0.047$)	$0.000$ ($0.000$)	$0.090$ ($0.042$)	$0.000$ ($0.000$)
KR	$0.115$ ($0.031$)	$1.308$ ($0.171$)	$0.112$ ($0.027$)	$1.979$ ($0.231$)
GPR	$0.075$ ($0.064$)	$111.417$ ($0.304$)	$0.020$ ($0.034$)	$114.394$ ($0.301$)
rBCM	$0.044$ ($0.052$)	$11.698$ ($0.085$)	$0.010$ ($0.025$)	$18.5709$ ($0.121$)

Equations

$$f(\vec{x}):={\vec{y}}_{i^{+}}\text{ where }{i}^{+}=\underset{i\in \{1,\dots ,N\}}{\mathrm{argmin}}d(\vec{x},{\vec{x}}_i)$$

$$f(\vec{x}):=\frac{{\sum }_{i=1}^{N}k(\vec{x},{\vec{x}}_i)\cdot {\vec{y}}_i}{{\sum }_{i=1}^{N}k(\vec{x},{\vec{x}}_i)}$$

$$\vec{k}$$

$$\mathbf{K}$$

$$p(\vec{y}_{1},\ldots,\vec{y}_{N},\vec{y}|\vec{x}_{1},\ldots,\vec{x}_{N},\vec{x})=\mathcal{N}\Big{(}\vec{\theta}_{1},\dots,\vec{\theta}_{N},\vec{\theta},\begin{pmatrix}\bm{K}+\tilde{\sigma}^{2}\cdot\bm{I}^{N}&{\vec{k}}^{T}\\ \vec{k}&k(\vec{x},\vec{x})\end{pmatrix}\Big{)}$$

$$p(\vec{y}$$

$$\vec{\mu }$$

$${\sigma }^2$$

$$s_d(x,x^{\prime }):=\exp \left(-\frac{1}{2}\frac{d(x,x^{\prime }{)}^2}{{\psi }^2}\right)$$

$$d(x_{i},x_{i^{\prime}})=\sqrt{{\Big{(}\phi(x_{i})-\phi(x_{i^{\prime}})\Big{)}}^{T}\cdot\Lambda\cdot\Big{(}\phi(x_{i})-\phi(x_{i^{\prime}})\Big{)}}$$

$$k(x_i,x_{i^{\prime }})={\varphi (x_i)}^T\cdot \varphi (x_{i^{\prime }})$$

$$\vec{\mu }={\vec{x}}^T+\vec{k}\cdot (\mathbf{K}+{\tilde{\sigma }}^2\cdot {\mathbf{I}}^N{)}^{-1}\cdot (\mathbf{Y}-\mathbf{X})$$

$$\frac{k(x,x_t^j)}{{\sum }_{j^{\prime }=1}^M{\sum }_{t^{\prime }=1}^{T_{j^{\prime }}-1}k(x,x_{t^{\prime }}^{j^{\prime }})}$$

$$\vec{\mu }$$

$$\displaystyle=\vec{x}+\sum_{j=1}^{M}-\gamma^{j}_{1}\cdot\vec{x}^{j}_{1}+\bigg{(}\sum_{t=2}^{T_{j}-1}(\gamma^{j}_{t-1}-\gamma^{j}_{t})\cdot\vec{x}^{j}_{t}\bigg{)}+\gamma^{j}_{T_{j}-1}\cdot\vec{x}^{j}_{T_{j}}$$

$${\Big{(}\phi(x_{i})-\sum_{i^{\prime}=1}^{N}\alpha_{i^{\prime}}\phi(x_{i^{\prime}})\Big{)}}^{T}\cdot\Lambda\cdot\Big{(}\phi(x_{i})-\sum_{i^{\prime}=1}^{N}\alpha_{i^{\prime}}\phi(x_{i^{\prime}})\Big{)}=\sum_{i^{\prime}=1}^{N}\alpha_{i^{\prime}}d(x_{i},x_{i^{\prime}})^{2}-\frac{1}{2}\vec{\alpha}\cdot\mathrm{D}^{2}\cdot{\vec{\alpha}}^{T}$$

$${\varphi (x_i)}^T\cdot \left(\sum _{i^{\prime }=1}^N{\alpha }_{i^{\prime }}\varphi (x_{i^{\prime }})\right)=\sum _{i^{\prime }=1}^N{\alpha }_{i^{\prime }}k(x_i,x_{i^{\prime }})$$

$${\sigma }_{\mathrm{rBCM}}^{-2}$$

$${\vec{\mu }}_{\mathrm{rBCM}}$$

$${\sigma }_{\mathrm{rBCM}}^2\cdot \left(\sum _{c=1}^C{\alpha }_c+{\alpha }_{\mathrm{prior}}\right)=1$$

$$\sum _{c=1}^C{\sigma }_{\mathrm{rBCM}}^2\cdot {\alpha }_c\cdot \left(\sum _{i=1}^N{\gamma }_i^c\right)=\sum _{c=1}^C{\sigma }_{\mathrm{rBCM}}^2\cdot {\alpha }_c\cdot 0=0$$

$$E=\sqrt{\frac{1}{T-1}\sum _{t^{\prime }=1}^{T-1}\sum _{j=1}^M\sum _{t=1}^{T_j}{\alpha }_{t^{\prime },t}^{j}d(x_t^j,x_{t^{\prime }+1}^{\prime }{)}^2+{\alpha }_{t^{\prime }}^{\prime }d(x_{t^{\prime }}^{\prime },x_{t^{\prime }+1}^{\prime }{)}^2-\frac{1}{2}{{\vec{\alpha }}_{t^{\prime }}}^T{\mathrm{D}}^2{\vec{\alpha }}_{t^{\prime }}}$$

$$\psi (v,t)$$

$$\rho ((u,v),t)$$

Full text

1 Introduction

To model connections between entities, graphs are oftentimes the method of choice, e.g. to model traffic connections between cities [47], data lines between computing nodes [13], communication between people in social networks [38], or the structure of a student’s solution to a learning task in an intelligent tutoring system [41, 45]. In all these examples, nodes as well as connections change significantly over time. For example, in traffic graphs, the traffic load changes significantly over the course of a day, making optimal routing a time-dependent problem [47]; in distributed computing, the distribution of computing load and communication between machines crucially depends on the availability and speed of connections and the current load of the machines, which changes over time [13]; in social networks or communication networks new users may enter the network, old users may leave and the interactions between users may change rapidly [38]; in intelligent tutoring systems, students change their solution over time to get closer to a correct solution [35, 41]. In all these cases it would be beneficial to predict the next state of the graph in question, because it provides the opportunity to intervene before negative outcomes occur, for example by re-routing traffic, providing additional bandwidth where required or to provide helpful hints to students.

Traditionally, predicting the future development based on knowledge of the past is the topic of time series prediction, which has wide-ranging applications in physics, sociology, medicine, engineering, finance and other fields [54, 58]. However, classic models in time series prediction, such as ARIMA, NARX, Kalman filters, recurrent networks or reservoir models focus on vectorial data representations, and they are not equipped to handle time series of graphs [58]. Accordingly, past work on predicting changes in graphs has focused on simpler sub-problems which can be phrased as classic problems, e.g. predicting the overall load in an energy network [2] or predicting the appearance of single edges in a social network [38].

In this contribution, we develop an approach to address the time series prediction problem for graphs, which we frame as a regression problem with structured data as input and as output. Our approach has two key steps: First, we represent graphs via pairwise kernel values, which are well-researched in the scientific literature [3, 12, 16, 17, 19, 57]. This representation embeds the discrete set of graphs in a smooth kernel space. Second, within this space, we can apply similarity- and kernel-based regression methods, such as nearest neighbor regression, kernel regression [43] or Gaussian processes [49] to predict the next position in the kernel space given the current position. Note that this does not provide us with the graph that corresponds to the predicted point in the kernel space. Indeed, identifying the corresponding graph in the primal space is a kernel pre-image problem which is in general hard to solve [4, 5, 36]. However, we will show that this data point can still be analyzed with subsequent kernel- or dissimilarity based methods.

If the underlying dynamics in the kernel space can not be captured by a simple regression scheme, such as 1-nearest neighbor or kernel regression, more complex models, such as Gaussian processes (GPs), are required. However, GPs suffer from a relatively high computational complexity due to a kernel matrix inversion. Fortunately, Deisenroth and Ng have suggested a simple strategy to permit predictions in linear time, namely distributing the prediction to multiple Gaussian processes, each of which handles only a constant-sized subset of the data [18]. Further, pre- and post-processing of kernel data - e.g. for eigenvalue correction - usually requires quadratic or cubic time. This added complexity can be avoided using the well-known Nyström approximation as investigated by [26].

The key contributions of our work are: First, we provide an integrative overview of seemingly disparate threads of research which are related to time-varying graphs. Second, we provide a scheme for time series prediction in dissimilarity and kernel spaces. This scheme is compatible with explicit vectorial embeddings, as are provided by some graph kernels [12, 16], but does not require such a representation. Third, we discuss how the predictive result, which is a point in an implicit kernel feature space, can be analyzed using subsequent kernel- or dissimilarity-based methods. Fourth, we provide an efficient realization of our prediction pipeline for Gaussian processes in linear time. Finally, we evaluate our proposed approaches on two theoretical and two practical data sets.

We start our investigation by covering related work and introducing notation to describe dynamics on graph data. Based on this work, we develop our proposed approach of time series prediction in kernel and dissimilarity spaces. We also discuss speedup techniques to provide predictions in linear time. Finally, we evaluate our approach empirically on four data sets: two well-established theoretical models of graph dynamics, namely the Barabasi-Albert model [6] and Conway’s Game of Life [23], as well as two real-world data sets of Java programs, where we try to predict the next step in program development [41, 42]. We find that for the theoretical models, even simple regression schemes may be sufficient to capture the underlying dynamics, but for the Java data sets, applying Gaussian processes significantly reduces the prediction error.

2 Background and Related Work

Dynamically changing graphs are relevant in many different fields, such as traffic [47], distributed computing [13], social networks [38] or intelligent tutoring systems [35, 41]. Due to the breadth of the field, we focus here on relatively general concepts which can be applied to a wide variety of domains. We begin with two formalisms to model dynamics in graphs, namely time-varying graphs [13], and sequential dynamical systems [8]. We then turn towards our research question: predicting the dynamics in graphs. This has mainly been addressed in the domain of social networks under the umbrella of link prediction [38, 60], as well as in models of graph growth [27]. Finally, as preparation for our own approach, we discuss graph kernels and dissimilarities as well as prior work on kernel-based approaches for time series prediction.

2.1 Models of Graph Dynamics

Time-Varying Graphs:

Time-varying graphs have been introduced by Casteigts and colleagues in an effort to integrate different notations found in the fields of delay-tolerant networks, opportunistic-mobility networks or social networks [13]. The authors note that, in all these domains, graph topology changes over time and the influence of such changes is too severe as to model them only in terms of system anomalies. Rather, dynamics have to be regarded as an “integral part of the nature of the system” [13]. We revisit a slightly varied version of the notation developed in their work.

Definition 1** (Time-Varying Graph [13]).**

A time-varying graph is defined as a five-tuple $\mathcal{G}=(V,E,\mathcal{T},\psi,\rho)$ where

• •

$V$ is an arbitrary set called nodes,

• •

$E\subseteq V\times V$ is a set of node tuples called edges,

• •

$\mathcal{T}=\{t:t_{0}\leq t\leq T,t\in\mathbb{N}\}$ is called lifetime of the graph,

• •

$\psi:V\times\mathcal{T}\to\{0,1\}$ is called node presence function, and node $v$ is called present at time $t$ if and only if $\psi(v,t)=1$, and

• •

$\rho:E\times\mathcal{T}\to\{0,1\}$ is called edge presence function, and edge $e$ is called present at time $t$ if and only if $\rho(e,t)=1$.

In figure 1, we show an example of a time-varying graph modeling the connectivity in simple public transportation graph over the course of a day. In this example, nodes model stations and edges model train connections between stations. In the night (left), all nodes may be present but no edges, because no lines are active yet. During the early morning (middle), some lines become active while others remain inactive. Finally, in mid-day (right), all lines are scheduled to be active, but due to a disturbance - e.g. construction work - a station is closed and all adjacent connections become unavailable. Note that in this example, all nodes, or stations, are assumed to be known in advance. If not all nodes are known in advance, e.g. a long-term public transportation model where new stations may be built in the future, the definition of time-varying graphs can be extended to allow an infinite node set $V$ with the restriction that only a finite number of nodes may be present at each time. This is also relevant for modelling social networks or student solutions in intelligent tutoring systems, where the introduction of new nodes happens frequently.

Using the notion of a presence function, many interesting concepts from static graph theory can be generalized to a dynamic version. For example, the neighborhood of a node $u$ at time $t$ can be defined as the set of all nodes $v$ with $\psi(v,t)=1$ for which an edge $(u,v)\in E$ exists such that $\rho((u,v),t)=1$. Similarly, a path at time $t$ between two nodes $u$ and $v$ can be defined as a sequence of edges $(u_{1},v_{1}),\ldots,(u_{K},v_{K})$, such that $u=u_{1}$, $v=v_{K}$ and for all $k$ $(u_{k},v_{k})\in E$, $u_{k+1}=v_{k}$, $\psi(u_{k},t)=1$, and $\rho((u_{k},v_{k}),t)=1$. Two nodes $u$ and $v$ can be called connected at time $t$ if a path between them exists at time $t$. Further, we can define the temporal subgraph $\mathcal{G}_{t}$ of graph $\mathcal{G}$ at time $t$ as the graph of all nodes and edges of $\mathcal{G}$ which are present at time $t$.

Note that we have assumed discrete time in our definition of a time-varying graph. This is justified by the following consideration: If a graph is embedded in continuous time, changes to the graph take the form of value changes in the node or edge presence function, which happen discretely, because a presence function can only take the values [math] or $1$. We call points in continuous time where such a discrete change occurs events. Assuming that there are only finitely many such events, we can write all events in the lifetime of a graph as an ascending sequence $t_{1},\ldots,t_{K}$. Accordingly, all changes in the graph are fully described by the sequence of temporal subgraphs $\mathcal{G}_{t_{1}},\ldots,\mathcal{G}_{t_{K}}$ [13, 55]. Therefore, even time-varying graphs defined on continuous time can be fully described by considering the discrete lifetime $\{1,\ldots,K\}$.

Sequential Dynamical Systems:

Sequential dynamical systems (SDS) have been introduced by Barret, Reidys and Mortvart to generalize cellular automata to arbitrary graphical structures [8, 10]. In essence, they assign a binary state to each node $v$ in a static graph $\mathcal{G}=(V,E)$. This state is updated according to a function $f_{v}$ which maps the current states of the node itself and all of its neighbors to the next state of the node $v$. This induces a discrete dynamical system on graphs (where edges and neighborhoods stay fixed) [8, 9, 10]. Interestingly, SDSs can be related to time-varying graphs by interpreting the binary state of a node $v$ at time $t$ as the value of its presence function $\psi(v,t)$. As such, sequential dynamical systems provide a compact model for the node presence function of a time-varying graph. Furthermore, if the graph dynamics is governed by a known SDS, its future can be predicted by simply simulating the dynamic system. Unfortunately, if the underlying SDS is unknown, this technique can not be applied, and to our knowledge, no learning schemes exists to date to infer an underlying SDS from data. Therefore, other predictive methods are required.

2.2 Predicting Changes in Graphs

In accordance with classical time series prediction, one can describe time series prediction in graphs as predicting the next temporal subgraph $\mathcal{G}_{t+1}$ given a sequence of temporal subgraphs $\mathcal{G}_{0},\ldots,\mathcal{G}_{t}$ in a time-varying graph. To our knowledge, there exists no approach which addresses this problem in this general form. However, more specific subproblems have been addressed in the literature.

Link Prediction:

In the realm of social network analysis, Liben-Nowell and Kleinberg have formulated the link prediction problem, which can be stated as: Given a sequence of temporal subgraphs $\mathcal{G}_{0},\ldots,\mathcal{G}_{t}$ for a time-varying graph $\mathcal{G}$, which edges will be added to the graph in the next time step, i.e. for which edges do we find $\rho(e,t)=0$ but $\rho(e,t+1)=1$ [38, 60]? For example, given all past collaborations in a scientific community, can we predict new collaborations in the future? The simplest approach to address this challenge is to compute a similarity index between nodes, rank all non-existing edges $(u,v)$ according to the similarity index of their nodes $s(u,v)$ and predict all edges which have an index above a certain threshold [38, 39]. Typical similarity indices to this end include the number of common neighbors at time $t$, the Jaccard index at time $t$ or the Adar index at time $t$ [38]. A more recent approach is to train a classifier that predicts the value of the edge presence function $\rho(e,t+1)$ for all edges with $\rho(e,t)=0$ using a vectorial feature representation of the edge $e$ at time $t$, where features include the similarity indices discussed above [39]. In a survey, Lü and Zhou further list maximum-likelihood approaches on stochastic models and probabilistic relational models for link prediction [40].

Growth models:

In a seminal paper, Barabási and Albert described a simple model to incrementally grow an undirected graph node by node from a small, fully connected seed graph [6]. Since then, many other models of graph growth have emerged, most notably stochastic block models and latent space models [14, 27]. Stochastic block models assign each node to a block and model the probability of an edge between two nodes only dependent on their respective blocks [34]. Latent space models embed all nodes in an underlying, latent space and model the probability of an edge depending on the distance in this space [31]. Both classes of models can be used for link prediction as well as graph generation. Further, they can be trained with pre-observed data in order to provide more accurate models of the data. However, none of the discussed models addresses the question of time series prediction in general because the deletion of nodes or edges or label changes are not covered by the models.

2.3 Graph Dissimilarities and Kernels

Instead of predicting the next graph in a sequence directly, one can consider indirect approaches, which base a prediction on pairwise dissimilarities $d(\mathcal{G},\mathcal{G}^{\prime})$ or kernels $k(\mathcal{G},\mathcal{G}^{\prime})$. As a simple example, consider a 1-nearest-neighbor approach: We first aggregate a database of training data consisting of graphs $\mathcal{G}_{t}$ and their respective successors $\mathcal{G}_{t+1}$. If we are confronted with a new graph $\mathcal{G}^{\prime}$, we simply predict the graph $\mathcal{G}_{t+1}$ such that $d(\mathcal{G}_{t},\mathcal{G}^{\prime})$ is minimized or $k(\mathcal{G}_{t},\mathcal{G}^{\prime})$ is maximized. The task of quantifying distance and similarity between graphs has attracted considerable research efforts, which we can divide roughly into two main streams: graph edit distances and graph kernels.

Graph Edit Distance:

The graph edit distance between two graphs $\mathcal{G}$ and $\mathcal{G}^{\prime}$ is traditionally defined as the minimum number of edit operations required to transform $\mathcal{G}$ into $\mathcal{G}^{\prime}$. Permitted edit operations include node insertion, node deletion, edge insertion, edge deletion, and the substitution of labels in nodes or edges [53]. This problem is a generalization of the classic string or tree edit distance, which is defined as the minimum number of operations required to transform a string into another or a tree into another respectively [37, 62]. Unfortunately, while the string edit distance and the tree edit distance can be efficiently computed in $\mathcal{O}(n^{2})$ and $\mathcal{O}(n^{4})$ respectively, computing the exact graph edit distance is NP-hard [61]. However, many approximation schemes exist, e.g. relying on self-organizing maps, Gaussian mixture models, graph kernels or binary linear programming [22]. A particularly simple approximation scheme is to order the nodes of a graph in a sequence and then apply a standard string edit distance measure on these sequences (refer e.g. to [45, 51, 52]). For our experiments on real-world Java data we will rely on a kernel over such an approximated graph distance as suggested in [45].

Note that for labeled graphs one may desire to assign non-uniform costs for substituting labels. Indeed, some labels may be semantically more related than others, implying a lower edit cost. Learning such edit costs from data is the topic of structure metric learning, which has mainly been investigated for string edit distances [11]. However, if string edit distance is applied as a substitute for graph edits, the results apply to graphs as well [45].

Graph Kernels:

Complementary to the view of graph edit distances, graph kernels represent the similarity between two graphs instead of their dissimilarity. The only formal requirement on a graph kernel is that it implicitly or explicitly maps a graph $\mathcal{G}$ to a vectorial feature representation $\phi(\mathcal{G})$ and computes the pairwise kernel values of two graphs $\mathcal{G}$ and $\mathcal{G}^{\prime}$ as the dot-product of their feature vectors, that is, $k(\mathcal{G},\mathcal{G}^{\prime})={\phi(\mathcal{G})}^{T}\cdot\phi(\mathcal{G}^{\prime})$ [17]. If each graph has a unique representation (that is, $\phi$ is injective) computing such a kernel is at least as hard as the graph isomorphy problem [24]. Thus, efficient graph kernels rely on a non-injective feature embedding, which is still expressive enough to capture important differences between graphs [17]. A particularly popular class of graph kernels are walk kernels which decompose the kernel between two graphs into kernels between paths which can be taken in the graphs [12, 17, 19]. More recently, advances have been made using decompositions of graphs into constituent parts instead of parts [16, 57]. The overall kernel is constructed via a sum of kernel values between the parts, which is guaranteed to preserve the kernel property according to the convolutional kernel principle [30]. In our experiments on artificial data, we apply the shortest-path-length kernel suggested by Borgwardt and colleagues, which compares the lengths of shortest paths in both graphs to construct an overall graph kernel, which is sufficiently expressive for these simple cases [12].

Any of these quantitative measures of distance or similarity supports time series prediction in a 1-nearest neighbor fashion as stated above. However, we would like to apply more sophisticated prediction methods as well. To this end, we turn to kernel-based methods.

2.4 Kernel-Based approaches for Vectorial Time Series Prediction

The idea to apply kernel-based methods for time series prediction as such is not new. One popular example is the use of support vector regression with wide-ranging applications in finance, business, environmental research and engineering [54]. Another example are gaussian processes to predict chemical processes [25], motion data [59] and physics data [50]. If a graph kernel is used that provides an explicit vectorial embedding of graphs, then these methods can be readily applied for time series prediction on graph data. However, in this work, we generalize this approach to kernels where only an implicit vectorial embedding is available. This enables us to use a broader range of graph distances and kernels, such as the graph edit distance approximation for our Java program data.

3 Time Series Prediction for Graphs

3.1 Time Series Prediction as a Regression Problem

Relying on the notation of time-varying graphs as introduced above, we can describe a time series of graphs as a sequence of temporal subgraphs $\mathcal{G}_{0},\ldots,\mathcal{G}_{t}$. Time series prediction is the problem of predicting the next graph in the series, $\mathcal{G}_{t+1}$. We first transform this problem into a regression problem as suggested by Sapankevych and colleagues, that is, we try to learn a function $f$ which maps the past $K$ states of a time series to a successor state [54]. Relying on distance- and kernel-based approaches we can address this problem by computing pairwise distance- or kernel values on graph sequences of length $K$, that is, $d((\mathcal{G}_{t-K},\ldots,\mathcal{G}_{t}),(\mathcal{G}^{\prime}_{t^{\prime}-K},\ldots,{\mathcal{G}^{\prime}}_{t^{\prime}}))$ or $k((\mathcal{G}_{t-K},\ldots,\mathcal{G}_{t}),(\mathcal{G}^{\prime}_{t^{\prime}-K},\ldots,{\mathcal{G}^{\prime}}_{t^{\prime}}))$ respectively. Using the sequential edit distance framework [45] or the convolutional kernel principle [17, 24], one can easily construct distances and kernels on graph sequences from distances and kernels from graphs. However, for simplicity, we will assume $K=1$ in the remainder of this paper, such that simple graph distances or kernels are sufficient. This is equivalent to the Markov assumption, that is, we assume that the next temporal subgraph $\mathcal{G}_{t+1}$ only depends on the temporal subgraph $\mathcal{G}_{t}$ and is conditionally independent from all temporal subgraphs $\mathcal{G}_{0},\ldots,\mathcal{G}_{t-1}$. Note that the Markov assumption is quite common in the related literature. For example, sequential dynamical systems make the very same assumption [8], as do most link prediction methods [40] and graph growth models [14, 27]. Further note that our framework generalizes to the non-Markovian case by replacing the graph distance or kernel to a sequence-of-graph distance or kernel and using everything else as-is.

Under the Markov assumption we can express our training data as tuples $\{(x_{i},y_{i})\}$ for $i\in\{1,\ldots,N\}$ where $y_{i}$ is the successor of $x_{i}$ in some time series, that is, $y_{i}=\mathcal{G}_{t+1}$ and $x_{i}=\mathcal{G}_{t}$. We denote the unseen test data point, for which a prediction is desired, as $x$.

We proceed by introducing three non-parametric regression techniques, namely $1$-nearest neighbor, kernel regression and Gaussian process regression. We then revisit some basic research on dissimilarities, similarities and kernels and leverage these insights to derive time series prediction purely in latent spaces without referring to any explicit vectorial form. Finally, we discuss how to speed up Gaussian process prediction from cubic to linear time.

3.2 Non-Parametric Regression techniques

In this section, we introduce three non-parametric regression techniques in their standard form, assuming vectorial input and output data. Further below we explain how these methods can be applied to structured input and output data. We assume that a training data set of tuples $\{(\vec{x}_{i},\vec{y}_{i})\}$ is available, where $\vec{y}_{i}$ is the desired output for input $\vec{x}_{i}$. Further, we assume that a dissimilarity $d$ or a kernel $k$ on the input space is available wherever needed.

1-Nearest Neighbor (1-NN):

We define the predictive function for one-nearest neighbor regression (1-NN) as follows:

[TABLE]

Note that the predictive function of 1-NN is not smooth, because the $\operatorname*{argmin}$ result is ill-defined if there exist two $i,i^{\prime}$ with $i\neq i^{\prime}$ and $\vec{y}_{i}\neq\vec{y}_{i^{\prime}}$ such that $d(\vec{x},\vec{x}_{i})=d(\vec{x},\vec{x}_{i^{\prime}})$. At these points, the predictive function is discontinuous and jumps between $\vec{y}_{i}$ and $\vec{y}_{i^{\prime}}$. A simple way to address this issue is kernel regression.

Kernel Regression (KR):

Kernel regression was first proposed by Nadaraya and Watson and can be seen as a generalization of 1-nearest neighbor to a smooth predictive function $f$ using a kernel $k$ instead of a dissimilarity [43]. The predictive function is given as:

[TABLE]

Note that kernel regression generally assumes kernel values to be positive and requires at least one $i$ for which $k(\vec{x},\vec{x}_{i})>0$, that is, if the test data point is not similar to any training data point, the prediction degenerates. Another limitation of kernel regression is that an exact reproduction of the training data is not possible, i.e. $f(\vec{x}_{i})\neq\vec{y}_{i}$. To achieve both, training data reproduction as well as a smooth predictive function, we turn to Gaussian process regression.

Gaussian Process Regression (GPR):

In Gaussian process regression (GPR) we assume that the output points (training as well as test) are a realization of a multivariate random variable with a Gaussian distribution [49]. The model extends KR in several ways. First, we can encode prior knowledge regarding the output points via the mean of our prior distribution, denoted as $\vec{\theta}_{i}$ and $\vec{\theta}$ for $\vec{y}_{i}$ and $\vec{y}$ respectively. Second, we can cover Gaussian noise on our training output points within our model. For this noise, we assume mean [math] and standard deviation $\tilde{\sigma}$.

Let now $k$ be a kernel on $\mathcal{X}$. We define

[TABLE]

Then under the GP model, the conditional density function of the output points given the input points is

[TABLE]

where $\bm{I}^{N}$ is the $N$-dimensional identity matrix and $\mathcal{N}(\vec{\mu},\bm{C})$ is the multivariate Gaussian probability density function for mean $\vec{\mu}$ and covariance matrix $\bm{C}$. Note that our assumed distribution takes all outputs $\vec{y}_{1},\ldots,\vec{y}_{N},\vec{y}$ as argument, not just a single point. The posterior distribution for just $\vec{y}$ can be obtained by marginalization:

Theorem 1** (Gaussian Process Posterior Distribution).**

Let $\bm{Y}$ be the matrix ${(\vec{y}_{1},\ldots,\vec{y}_{N})}^{T}$ and $\bm{\Theta}:={(\vec{\theta}_{1},\ldots,\vec{\theta}_{N})}^{T}$. Then the posterior density function for Gaussian process regression is given as:

[TABLE]

We call $\vec{\mu}$ the predictive mean and $\sigma^{2}$ the predictive variance.

Proof.

Refer e.g. to [49, p. 27]. ∎

Note that the posterior distribution is, again, a Gaussian distribution. For such a distribution, the mean corresponds to the point of maximum probability density, such that we can define our predictive function as $f(\vec{x}):={\vec{\mu}}^{T}$ where $\vec{\mu}$ is the predictive mean of the posterior distribution for point $\vec{x}$. Further note that the predictive mean becomes the prior mean if $\vec{k}$ is the zero vector, i.e. if the test data point is not similar to any training data point.

The main drawback of GPs is their high computational complexity: For training, the inversion of the matrix $(\bm{K}+\tilde{\sigma}^{2}\cdot\bm{I}^{N})^{-1}$ requires cubic time. We will address this problem in section 3.5.

3.3 Dissimilarities, Similarities and Kernels

Each of the methods described above relies either on a dissimilarity (in the case of 1-NN), or on a kernel (in the case of KR and GPR). Note that we have not strictly defined these concepts until now. Indeed, dissimilarities as well as similarities are inherently ill-defined concepts. A dissimilarity for the set $\mathcal{X}$ is any function of the form $d:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$, such that $d$ decreases if the two input arguments are in some sense more related to each other [48]. Conversely, a similarity for the set $\mathcal{X}$ is any function $s:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ which increases if the two input arguments are in some sense more related to each other [48]. While in many cases, more rigorous criteria apply as well (such as non-negativity, symmetry, or that any element is most similar to itself) any or all of them may be violated in practice [56]. On the other hand, kernels are strictly defined as inner products. If we require a kernel to apply our prediction method but only have a dissimilarity or a similarity available, we require a conversion method to obtain a corresponding kernel. In a first step we note that it is simple to convert from a dissimilarity $d$ to a similarity $s$. We can simply apply any monotonously decreasing function. In the context of this work, we will apply the radial basis function transformation

[TABLE]

where $\psi$ is a positive real number called bandwidth.

Further, for any finite set of data points $x_{1},\ldots,x_{N}$, we can transform a symmetric similarity into a kernel. Let $\bm{S}$ be a symmetric matrix of pairwise similarities for these data points with the entries $\bm{S}_{ii^{\prime}}=s(x_{i},x_{i^{\prime}})$. Such a symmetric similarity matrix $\bm{S}$ is a kernel matrix if and only if it is positive semi-definite [48]. Positive semi-definiteness can be enforced by applying an eigenvalue decomposition $\bm{S}=\bm{U}\bm{\Lambda}\bm{V}$ and removing the negative eigenvalues in the matrix $\bm{\Lambda}$, e.g. by setting them to zero or by taking their absolute value. The resulting matrix $\bm{\tilde{\Lambda}}$ can then be applied to obtain a kernel matrix $\bm{K}=\bm{U}\bm{\tilde{\Lambda}}\bm{V}$ [26, 48, 56]. Note that the eigenvalue decomposition of a similarity matrix has cubic time-complexity in the number of data points. However, linear-time approximations have recently been discovered based on the Nyström method [26]. Therefore, in our further work, we can assume that a kernel can be obtained as needed.

3.4 Prediction in Kernel and Dissimilarity Spaces

In introducing the three regression techniques above we have assumed that each data point is a vector, such that algebraic operations like scalar multiplication (for kernel regression), as well as matrix-vector multiplication and vector addition (for Gaussian process regression) are permitted. In the case of graphs, such operations are not well-defined. However, we can rely on prior work regarding vectorial embeddings of dissimilarities and similarities to define our regression in a latent space. In particular:

Theorem 2** (Pseudo-Euclidean Embeddings [48]).**

For dissimilarities: Let $x_{1},\ldots,x_{N}$ be some points from a set $\mathcal{X}$, and let $d:\mathcal{X}\times\mathcal{X}\to\mathbb{R}^{+}$ be a function such that for all $x_{i},x_{i^{\prime}}$ it holds: $d(x_{i},x_{i^{\prime}})\geq 0$, $d(x_{i},x_{i})=0$ and $d(x_{i},x_{i^{\prime}})=d(x_{i^{\prime}},x_{i})$. Then, there exists a vector space $\mathcal{Y}\subset\mathbb{R}^{R}$ and a function $\phi:\mathcal{X}\to\mathcal{Y}$ such that for all $i,i^{\prime}\in\{1,\ldots,N\}$ it holds:

[TABLE]

for some diagonal matrix $\Lambda\in\{-1,0,1\}^{R\times R}$.

For kernels: Let $k$ be a kernel on $\mathcal{X}$. Then there exists an Euclidean space $\mathcal{Y}$ and a function $\phi:\mathcal{X}\to\mathcal{Y}$ such that for all $i,i^{\prime}\in\{1,\ldots,N\}$ it holds:

[TABLE]

Proof.

Refer to [48]. ∎

The key idea of our approach to time series prediction for structured data (i.e. graphs) is to apply time series prediction in the implicit (pseudo-)Euclidean space $\mathcal{Y}$, without any need for an explicit embedding in $\mathcal{Y}$. This is obvious for 1-NN: Our prediction is just the successor to the closest training data point. This is a graph itself, providing us with a function mapping directly from graphs as inputs to graphs as outputs. The situation is less clear for KR or GPR. Here, our approach is to express the predictive output as a linear combination of known data points, which permits further processing as we will show later. In particular, we can proof that the resulting linear combinations are convex or at least affine. Linear, convex and affine combinations are defined as follows:

Definition 2** (Linear, Affine, and Convex Combinations).**

Let $\vec{x}_{1},\ldots,\vec{x}_{N}$ be vectors from some vector space $\mathcal{X}$ and let $\alpha_{1},\ldots,\alpha_{N}$ be real numbers. The sum $\sum_{i=1}^{N}\alpha_{i}\cdot\vec{x}_{i}$ is called a linear combination and the numbers $\alpha_{i}$ linear coefficients. If it holds $\sum_{i=1}^{N}\alpha_{i}=1$ the linear combination is called affine. If it additionally holds that for all $i$ $\alpha_{i}\geq 0$ the affine combination is called convex.

Our proof that the prediction provided by GPR is an affine combination requires two additional assumptions. First, we require a certain, natural prior, namely: In the absence of any further knowledge, our best prediction for the point $x$ is the identity, i.e. staying where we are. In this case, we obtain $\vec{\theta}=\vec{x}$ and $\vec{\theta}_{i}=\vec{x}_{i}$ and the predictive mean

[TABLE]

where $\bm{X}={(\vec{x}_{1},\ldots,\vec{x}_{N})}^{T}$. The predictive variance is still the same as in equation 8. We will assume this prior in the remainder of this paper.

Further, we observe that our training inputs and outputs for time series prediction have a special form: Training data is presented in the form of sequences $x^{j}_{1},\ldots,x^{j}_{T_{j}}$ for $j\in\{1,\ldots,M\}$ where $x^{j}_{t+1}$ is the successor of $x^{j}_{t}$. This implies that each point that is a predecessor of another point is also a successor of another point (except for the first element in each sequence) and vice versa (except for the last element in each sequence).

Theorem 3** (Predictive results as combinations).**

If training data is provided in terms of sequences $x^{j}_{1},\ldots,x^{j}_{T_{j}}$ for $j\in\{1,\ldots,M\}$ where $x^{j}_{t+1}$ is the successor of $x^{j}_{t}$, then it holds:

The predictive result for 1-NN is a single training data point (i.e. a convex combination of only a single data point). 2. 2.

The predictive result for KR is an affine combination of training data points. 3. 3.

If the prior $\theta=x$ and $\theta^{j}_{t}=x^{j}_{t}$ is assumed, the predictive result for GPR is an affine combination of training data points and the test data point, where the test data point has the coefficient $1$.

Proof.

This follows directly from the form of the predictive function in equation 1 2. 2.

The coefficients are given as

[TABLE]

for any point $x^{j}_{t+1}$ with $j\in\{1,\ldots,M\}$ and $t\in\{1,\ldots,T_{j}-1\}$. This combination is affine due to the normalization by $\sum_{j^{\prime}=1}^{M}\sum_{t^{\prime}=1}^{T_{j^{\prime}}-1}k(x,x^{j^{\prime}}_{t^{\prime}})$. If the kernel is non-negative, the combination is convex. 3. 3.

We define $\vec{\gamma}=(\gamma^{1}_{1},\ldots,\gamma^{1}_{T_{1}-1},\ldots,\gamma^{M}_{1},\ldots,\gamma^{M}_{T_{M}-1}):=\vec{k}\cdot(\bm{K}+\tilde{\sigma}^{2}\cdot\bm{I}^{N})^{-1}$. The predictive mean is given as $\vec{\mu}=\vec{x}+\vec{\gamma}\cdot(\bm{Y}-\bm{X})$ which yields

[TABLE]

Thus, the coefficients within each sequence are $-\gamma^{j}_{1}$, $\gamma^{j}_{1}-\gamma^{j}_{2}$, $\ldots$, $\gamma^{j}_{T_{j}-2}-\gamma^{j}_{T_{j}-1}$, $\gamma^{j}_{T_{j}-1}$. These coefficients add up to zero. Finally, we have a coefficient of $1$ for $\vec{x}$, which is our prior for the prediction. Therefore, the overall sum of all coefficients is $1$ and the combination is affine.

∎

In effect, we can reformulate the predictive function as $f:\mathcal{X}\to\mathbb{R}^{N+1}$ mapping a test data point to the coefficients of an affine combination, which represent our actual predictive result in the (pseudo-)Euclidean space $\mathcal{Y}$. Note that these coefficients do not provide us with a primal space representation of our data, i.e. we do not know what the graph which corresponds to a particular affine combination looks like. Indeed, as graph kernels are generally not injective [24], there might be multiple graphs that correspond to the same affine combination. Finding such a graph is called the kernel pre-image problem and is hard to solve even for vectorial data [4, 5, 36]. This poses a challenge with respect to further processing: How do we interpret a data point for which we have no explicit representation [33]?

Fortunately, we can still address many classical questions of data analysis for such a representation relying on our already existing implicit embedding in $\mathcal{Y}$. We can simply extend this embedding for the predicted point via the affine combination as follows:

Theorem 4** (Distances and Kernels for Affine Combinations).**

For dissimilarities: Let $x_{1},\ldots,x_{N}$, $\phi$, $d$ and $\Lambda$ be as in theorem 2. Let $\mathrm{D}^{2}$ be the matrix of squared pairwise dissimilarities between the points $x_{1},\ldots,x_{N}$, and let $\vec{\alpha}$ be a $1\times N$ vector of affine coefficients. Then it holds for any $i\in\{1,\ldots,N\}$:

[TABLE]

For kernels: Let $x_{1},\ldots,x_{N}$, $\phi$ and $k$ be as in theorem 2 and let $\bm{K}$ be the matrix of pairwise kernel values between the points $x_{1},\ldots,x_{N}$. Let $\vec{\alpha}$ be a $1\times N$ vector of linear coefficients. Then it holds for any $i\in\{1,\ldots,N\}$:

[TABLE]

Proof.

Refer to [28] for a proof of 16. 17 follows by simple linear algebra. ∎

Using this extended embedding, we can simply apply any dissimilarity- or kernel-based method on the predicted point. For dissimilarities, this includes relational learning vector quantization for classification [29] or relational neural gas for clustering [28]. For kernels, this includes the non-parametric regression techniques discussed here, but also techniques like kernel vector quantization [32] or support vector machines for classification [15] and kernel variants of $k$-means, SOM and neural gas for clustering [20]. Therefore, we have achieved a full methodological pipeline for preprocessing, prediction and post-processing, which can be summarized as follows:

If we intend to use a dissimilarity measure on graphs, we start off by computing the matrix of pairwise dissimilarities $\mathrm{D}$ on our training data. If required we symmetrize this matrix by setting $\mathrm{D}\leftarrow\frac{1}{2}\cdot(\mathrm{D}+{\mathrm{D}}^{T})$ and set the diagonal to zero. Implicitly, this step embeds our training data in a pseudo-Euclidean space $\mathcal{Y}$ where $\mathrm{D}$ are pairwise pseudo-Euclidean distances. We transform this matrix into a similarity matrix $\bm{S}$ using, for example, the radial basis function transformation. 2. 2.

If we intend to use a similarity measure on graphs, we start off by computing the matrix of pairwise similarities $\bm{S}$ on our training data. Otherwise we use the transformed matrix from the previous step. 3. 3.

We transform $\bm{S}$ to a kernel matrix $\bm{K}$ via Eigenvalue correction [26, 33, 48]. Implicitly, this step embeds our data in an Euclidean space $\mathcal{Y}$ where our pairwise kernel values represent inner products of our training data. 4. 4.

For any test data point $x$ we compute the vector of kernel values $\vec{k}$ to all training data points. 5. 5.

We apply 1-NN, KR or GPR as needed to infer a prediction for our test data point $x$ in form of an affine coefficient vector $\vec{\alpha}$. 6. 6.

We extend our dissimilarity matrix and/or kernel matrix for the predicted point via equations 16 and/or 17 respectively. 7. 7.

We apply any further dissimilarity- or kernel-based method on the predicted point as desired.

The only challenge left is to speed up predictions in GPR to reduce the cubic time complexity to linear time complexity.

3.5 Linear Time Predictions

GP regression involves the inversion of the matrix $(\bm{K}+\tilde{\sigma}^{2}\cdot\bm{I}^{N})$, resulting in $\mathcal{O}(N^{3})$ complexity. A variety of efficient approximation schemes exist [49]. Recently, the robust Bayesian Committee Machine (rBCM) has been introduced as a particularly fast and accurate approximation [18]. The rBCM approach is to distribute the examples into $C$ disjoint sets, based e.g. on clustering in the input data space. For each of these sets, a separate GP regression is used, yielding the predictive distributions $\mathcal{N}(\vec{\mu}_{c},\sigma^{2}_{c})$ for $c\in\{1,\ldots,C\}$. These distributions are combined to the final predictive distribution $\mathcal{N}(\vec{\mu}_{\mathrm{rBCM}},\sigma^{2}_{\mathrm{rBCM}})$ with

[TABLE]

Here, $\sigma_{\mathrm{prior}}^{2}$ is the variance of the prior for the prediction, which is a new meta-parameter introduced in the model. The weights $\beta_{c}$ can be seen as a measure for the predictive power of the single GP experts. As suggested by the authors, we use the differential entropy, given as \beta_{c}=\frac{1}{2}\cdot\Big{(}\log(\sigma_{\mathrm{prior}}^{2})-\log(\sigma^{2}_{c})\Big{)} [18]. Note that $\vec{\theta}=\vec{x}$ in our case. This approach results in linear-time complexity if the size of any single cluster is considered to be constant (i.e. the number of clusters is proportional to $N$), such that only a matrix of constant size has to be inverted.

Two challenges remain to apply rBCM productively within our proposed pipeline. First, it remains to show that the predictive mean of rBCM still has the form of an affine combination. Second, we require a dissimilarity- or kernel-based clustering scheme which runs in linear time, as to ensure overall linear-time complexity. Regarding the first issue we show:

Theorem 5** (rBCM prediction as affine combination).**

If training data is provided in terms of sequences $x^{j}_{1},\ldots,x^{j}_{T_{j}}$ for $j\in\{1,\ldots,M\}$ where $x^{j}_{t+1}$ is the successor of $x^{j}_{t}$, the predictive mean of an rBCM is an affine combination of training data points and the test data point.

Proof.

We define $\alpha_{c}:=\frac{\beta_{c}}{\sigma_{c}^{2}}$ and \alpha_{\mathrm{prior}}:=\bigg{(}1-\sum_{c=1}^{C}\beta_{c}\bigg{)}\cdot\frac{1}{\sigma_{\mathrm{prior}}^{2}}. For every cluster, the respective predictive mean has the shape $\vec{\mu}_{c}=\vec{x}+\sum_{i=1}^{N}\gamma^{c}_{i}\cdot\vec{x}_{i}$ where $\sum_{i=1}^{N}\gamma^{c}_{i}=0$ (as demonstrated in the proof to theorem 3). So the overall coefficient assigned to $\vec{x}$ is

[TABLE]

and all other coefficients add up to

[TABLE]

Therefore, we obtain an affine combination. ∎

With regards to the second issue, we apply relational neural gas (RNG) [28]. RNG clusters data points using a Voronoi-tesselation with respect to prototypes. To ensure roughly constant-sized clusters we have to supply the method with a number of prototypes proportional to the data set size. Further, we have to ensure that the clustering itself takes only linear time. As such, relational neural gas requires quadratic time for training. However, it can be sped up by training the prototype positions only on a constant-sized subset of the data and assigning all remaining data points to the closest prototype. Computing the distance of all data points to all prototypes requires only the pairwise distances of all data points to the data points in the training subset, which is in $\mathcal{O}(N)$.

The approach is illustrated in figure 2. Relational neural gas places prototypes (large circle and square) into the data (small circles and squares) and thereby distributes data points into disjoint clusters (distinguished by shape). Successor relations are depicted as arrows. For each cluster, a separate Gaussian process is trained. For a test data point (diamond shape), each of the GPs provides a separate predictive Gaussian distribution, which are given in terms of their means (half-tansparent circle and square) and their variance (dashed, half-transparent circles). The predictive distributions are merged to an overall predictive distribution with the mean from equation 19 (solid diamond shape) and the variance from equation 18 (dashed circle). Note that the overall predictive distribution is more similar to the prediction of the circle-cluster because the test data point is closer to this cluster and thus the predictive variance for the circle-cluster is lower, giving it a higher weight in the merge process.

4 Experiments

In our experimental evaluation, we apply the pipeline introduced in the previous section to four data sets, two of which are theoretical models and two of which are real-world data sets of Java programs. In all cases, we evaluate the root mean square error (RMSE) of the prediction for each method in a leave-one-out-crossvalidation over the sequences in our data set. We denote the current test trajectory as $x^{\prime}_{1},\ldots,x^{\prime}_{T}$, the training trajectories as $\{x^{j}_{1},\ldots,x^{j}_{T_{j}}\}_{j=1,\ldots,M}$, the predicted affine coefficients for point $x^{\prime}_{t^{\prime}}$ as $\vec{\alpha}_{t^{\prime}}=(\alpha^{1}_{t^{\prime},1},\ldots,\alpha^{M}_{t^{\prime},T_{M}},\alpha^{\prime}_{t^{\prime}})$ and the matrix of squared pairwise dissimilarities (including the test data points) as $\mathrm{D}^{2}$. Accordingly, the RMSE for each fold has the following form (refer to equation 16).

[TABLE]

We evaluate our four regression models, namely 1-nearest neighbor (1-NN), kernel regression (KR), Gaussian process regression (GPR) and the robust Bayesian committee machine (rBCM), as well as the identity function as baseline, i.e. we predict the current point as next point. We optimized the hyper parameters for all methods (i.e. the radial basis function bandwidth $\psi$ and the noise standard deviation $\tilde{\sigma}$ for GPR and rBCM) using a random search with $10$ random trials111Our implementation of time series prediction is available online at http://doi.org/10.4119/unibi/2913104. In each trial, we evaluated the RMSE in a nested leave-one-out-crossvalidation over the training sequences and chose the parameters which corresponded to the lowest RMSE. Let $\bar{d}$ be the average dissimilarity over the data set. We drew $\psi$ from a uniform distribution in the range $[0.05\cdot\bar{d},\bar{d}]$ for the theoretical data sets and fixed it to $0.3\cdot\bar{d}$ for the Java data sets to avoid the need for a new eigenvalue correction in each random trial. We drew $\tilde{\sigma}$ from an exponential distribution in the range $[10^{-3}\cdot\bar{d},\bar{d}]$ for the theoretical and $[10^{-2}\cdot\bar{d},\bar{d}]$ for the Java data sets. We fixed the prior standard deviation $\sigma_{\mathrm{prior}}=\bar{d}$ for all data sets. For rBCM we preprocessed the data via relational neural gas (RNG) clustering with $\left\lfloor\frac{N}{100}\right\rfloor$ clusters for all data sets. As this pre-processing could be applied before hyper-parameter selection the runtime overhead of clustering was negligible and we did not need to rely on the linear-time speedup described above but could compute the clustering on the whole training data set.

Our experimental hypotheses are that all prediction methods should yield lower RMSE compared to the baseline (H1), that rBCM should outperform 1-NN and KR (H2) and that rBCM should be not significantly worse compared to GPR (H3). To evaluate significance we use a Wilcoxon signed rank sum test.

4.1 Theoretical Data Sets

We investigate the following theoretical data sets:

Barabási-Albert model:

This is a simple stochastic model of graph growth in undirected graphs with hyper-parameters $m_{0}$, $k$ and $m$ [6]. The growth process starts with a fully connected initial graph of $m_{0}$ nodes and adds $m-m_{0}$ nodes one by one. Each newly added node is connected to $k$ of the existing nodes which are randomly selected with the probability $P(u)=\deg_{t}(u)/(\sum_{v}\deg_{t}(v))$ where $\deg_{t}$ is the node degree at time $t$, i.e. $\deg_{t}(v)=\sum_{u}\rho((u,v),t)$. It has been shown that the edge distribution resulting from this growth model is scale-free, more specifically the probability of a certain degree $k$ is $P(k)\propto k^{-3}$ [6]. Our data set consists of $20$ graphs with $m=27$ nodes each, grown from an initial graph of size $m_{0}=3$ and $k=2$ new edges per node. This resulted in $500$ graphs overall.

Conway’s Game of Life:

John Conway’s Game of Life [23] is a simple, 2-dimensional cellular automaton model. Nodes are ordered in a regular, 2-dimensional grid and connected to their eight neighbors in the grid. Let $\mathcal{N}(v)$ denote this eight- neighborhood in the grid. Then we can describe Conway’s Game of Life with the following equations for the node presence function $\psi$ and the edge presence function $\rho$ respectively:

[TABLE]

Note that Conway’s Game of Life is turing-complete and its evolution is, in general, unpredictable without computing every single step according to the rules [1]. We created $30$ trajectories by initializing a $20\times 20$ grid with one of six standard patterns at a random position, namely blinker, beacon, toad, block, glider, and block and glider (see figure 3). The first four patterns are simple oscillators with a period of two, the glider is an infinitely moving structure with a period of two (up to rotation) and the block and glider is a chaotic structure which converges to a block of four and a glider after 105 steps. 222Also refer to the Life Wiki http://conwaylife.com/wiki/ for more information on the patterns. We let the system run for $T=10$ time steps resulting in $300$ graphs overall. In every step, we further activated $5\%$ of the cells at random, simulating observational noise.

As data representation for the theoretical data sets we use an explicit feature embedding inspired by the shortest-path-kernel of Borgwardt and colleagues [12]. Using the Floyd-Warshall algorithm [21] we compute all shortest paths in the graph and then compute a histogram over the lengths of these shortest paths as a feature map (see figure 4 for an example). As dissimilarity we use the Euclidean distance on these feature vectors, which we normalize by the average distance over the data set. We transformed the distance to a kernel via the radial basis function transformation. Eigenvalue correction was not required.

The RMSE and runtimes for all three data sets are shown in table 1. As expected, KR, GPR and rBCM outperform the identity-baseline ($p<10^{-3}$ for both data sets), supporting H1. 1-NN outperforms the baseline only in the Barabási-Albert data set ($p<10^{-3}$). Also, our results lend support to H2 as rBCM outperforms 1-NN in both data sets ($p<0.05$ for Barabási-Albert, and $p<0.01$ for Conway’s Game of Life). However, rBCM is significantly better than KR only for the Barabási-Albert data set ($p<0.001$), indicating that for simple data sets such as our theoretical ones, KR might already provide a sufficient predictive quality. Finally, we do not observe a significant difference between rBCM and GPR, as expected in H3. Interestingly, for these data sets, rBCM is slower compared to GP, which is probably due to constant overhead for maintaining multiple models.

4.2 Java Programs

Our two real-world java data sets consist of programs for two different problems from beginner’s programming courses. The motivation for time series prediction on such data is to help students achieve a correct solution in an intelligent tutoring system (ITS). In such an ITS, students incrementally work on their program until they might get stuck and do not know how to proceed. Then, we would like to predict the most likely next state of their program, given the trajectories of other students who have already correctly solved the problem. This scheme is inspired by prior work in intelligent tutoring systems, such as the hint factory by Barnes and colleagues [7].

As our data sets consist only of final, working versions of programs, we have to simulate the incremental growth. To that end we first represented the programs as graphs via their abstract syntax tree, and then recursively removed the last semantically important node (where we regarded nodes as important which introduce a visibility scope in the Java program, such as class declarations, method declarations, and loops), until the program was empty. Then, we reversed the order of the resulting sequence, as to achieve a growing program. In detail, we have the following two data sets:

MiniPalindrome:

This data set consists of $48$ Java programs, each realizing one of eight different strategies to recognize palindromic input (see figure 5) [41] 333The data set is available online at http://doi.org/10.4119/unibi/2900666. The abstract syntax trees of these programs contain $135$ nodes on average. The programs come in eight different variations described in [41]. Our simulation resulted in $834$ data points.

Sorting:

This is a benchmark data set of $64$ Java sorting programs taken from the web, implementing one of two sorting algorithms, namely BubbleSort or InsertionSort (see figure 5) [42] 444The data set is available online at http://doi.org/10.4119/unibi/2900684. The abstract syntax trees contain $94$ nodes on average. Our simulation resulted in $800$ data points.

To achieve a dissimilarity representation we first ordered the nodes of the abstract syntax trees in order of their appearance in the original program code. Then, we computed a sequence alignment distance on the resulting node sequences, similar to the method described by Robles-Kelly and Hancock [52]. In particular, we used an affine sequence alignment with learned node dissimilarity as suggested in [45]. We transformed the dissimilarity to a similarity via the radial basis function transformation and obtained a kernel via clip Eigenvalue correction [26].

We show the RMSEs and runtimes for both data sets in table 2. Contrary to the theoretical data sets before, we observe strong differences both between the predictive models and the baseline, as well as between rBCM and 1-NN as well as KR ($p<0.01$ in all cases), which supports both H1 and H2. Interestingly, rBCM appears to achieve better results compared to GPR, which might be the case due to additional smoothing provided by the averaging operation over all cluster-wise GPR results. This result supports H3. Finally, we observe that rBCM is about 10 times faster compared to GPR.

5 Discussion and Conclusion

Our results indicate that it is possible to achieve time series prediction in kernel and dissimilarity spaces, in particular for graphs. In all our experiments, even simple predictive models (1-nearest neighbor and kernel regression) did outperform the baseline of staying where you are. For real-world data we could further improve the predictive error significantly by applying a more complex predictive model, namely the robust Bayesian Committee Machine (rBCM). This indicates a trade-off in terms of model choice: Simpler models are faster and in the case of 1-nearest neighbor the result is easier to interpret, because it is given as a graph. However, in the case of real-world data, it is likely that a more complex predictive model is required to accurately describe the underlying dynamics in the kernel space. Fortunately, the runtime overhead is only a constant factor as rBCM can be applied in linear time.

Our key idea is to apply time series prediction in a (pseudo-)Euclidean space and representing the graph output as a point in this space. Even though we do not know the graph corresponding to the predicted location in the latent space, we have shown that this point can be analyzed using subsequent dissimilarity- or kernel-based methods as the dissimilarities and kernel values with respect to the predicted point can still be calculated.

However, a few challenges for further research remain. First, usual hyperparameter optimization techniques for Gaussian processes depend on a vectorial data representation [18] and thus are not necessarily applicable in our proposed pipeline. Therefore, alternative hyperparameter selection techniques are required. Second, theoretic or empirical results regarding the number of data required to make accurate predictions are still lacking for this novel approach. Finally, for some applications, the predicted point in the primal space may be required, that is, we need a prediction in form of a graph, for example for feedback provision in intelligent tutoring systems or predicting the precise structure of a sensor network, e.g. for predictive maintenance in an oil refinery. For such cases, we have to solve a pre-image problem: Finding the original point that maps to the affine combination in the pseudo-Euclidean space. This problem has proved particularly challenging in the literature up until now and could profit from further consideration [4, 5, 36].