Comparing temporal graphs using dynamic time warping

Within many real-world networks, the links between pairs of nodes change over time. Thus, there has been a recent boom in studying temporal graphs. Recognizing patterns in temporal graphs requires a proximity measure to compare different temporal graphs. To this end, we propose to study dynamic time warping on temporal graphs. We define the dynamic temporal graph warping (dtgw) distance to determine the dissimilarity of two temporal graphs. Our novel measure is flexible and can be applied in various application domains. We show that computing the dtgw-distance is a challenging (in general) NP-hard optimization problem and identify some polynomial-time solvable special cases. Moreover, we develop a quadratic programming formulation and an efficient heuristic. In experiments on real-world data, we show that the heuristic performs very well and that our dtgw-distance performs favorably in de-anonymizing networks compared to other approaches.


Introduction
A fundamental concept for pattern recognition is the notion of proximity (i.e., (dis)similarity) between objects. For objects that are represented by numerical feature vectors, there exist a lot of well-known proximity measures such as p-norms or positive semi-definite kernels. In structural pattern recognition, objects are often more naturally represented by complex (discrete) data structures such as graphs, strings, or time series. For these representations, one can often not simply use vector-based proximity measures. Instead, one needs to define suitable domain-specific proximity measures such as the edit distance on graphs or strings or the dynamic time warping distance 1 on time series.
The majority of graph proximity measures focus on static graphs. This includes the graph edit distance (Riesen 2015), graph kernels (Kriege et al. 2020), and geometric graph distances (Jain 2016). However, many complex systems are not static as the links between entities dynamically change over time. Thus, there is a steadily growing research interest in analyzing temporal graphs (we also use the term temporal network interchangeably) (Rozenshtein and Gionis 2019;Saramäki 2013, 2019). Such temporal graphs can be represented by a series of temporal edges between a fixed set of vertices. Examples are social contact networks, disease spreading networks, traffic networks, attack networks in computer security, or protein-protein interaction networks in biology (Kostakos 2009;Saramäki 2013, 2019;Li et al. 2017;Vijayan et al. 2017). Examples of data mining problems on temporal social networks include community detection (Dakiche et al. 2019), epidemics analysis (Rozenshtein et al. 2016), and influence spreading (Gayraud et al. 2015).
An extended abstract of this article appeared in the Proceedings of the 8th International Conference on Complex Networks and their Applications, Springer, SCI, vol 882, pp. 469-480, 2019(Froese et al. 2019. This article contains all proofs in full detail and presents extended experimental results. Many processes described by temporal graphs naturally vary in duration and temporal dynamics (for example, chemical reactions or the spread of a disease might proceed at different speeds), which makes data mining tasks such as classification challenging. Hence, one needs to find suitable proximity measures, which has seemingly not been done so far.
Our paper proposes such a measure. We introduce a novel proximity measure on temporal graphs based on vertex signature graph distance and dynamic time warping, called dynamic temporal graph warping (dtgw). Dynamic time warping allows to cope with variations in temporal dynamics. Thus, by combining established methods from graphbased pattern recognition and time series data mining in a nontrivial way, we obtain a suitable tool to analyze temporal network data. We study the computational complexity of the dtgw-distance, develop efficient algorithms and study their behavior on real-world data, the latter indicating the strong potential for future applications.
Related Work. Graph distance based on vertex mappings using local vertex signatures was introduced by Jouili and Tabbone (2009). The idea of using vertex mappings can also be found in optimal assignment kernels (Fröhlich et al. 2005;Kriege et al. 2016;Bento and Ioannidis 2018).
Regarding proximity measures on temporal graphs, seemingly little work has been done so far. In fact, we are not aware of other approaches for numerically measuring the proximity of two temporal graphs. Related concepts, however, have been investigated for temporal graphs. For example, one approach is based on network embeddings where nodes are mapped onto certain feature spaces incorporating the temporal behavior (Ahmed and Karypis 2015;Zuo et al. 2018;Nguyen et al. 2018). Another approach is based on network alignments (Vijayan et al. 2017;Elhesha et al. 2019) where a vertex mapping that optimizes some criteria is computed. However, dynamic time warping has not been used in this context so far.
Dynamic time warping (Sakoe and Chiba 1978) is an established measure for mining time series data (Rakthanmanon et al. 2012;Wang et al. 2013) which is specifically designed to cope with temporal distortion in the data via nonlinear alignment of time series. It can be applied to time series of different lengths (in contrast to the Euclidean distance for example) which is a relevant aspect in time series averaging. We lift this approved concept to the domain of temporal graphs.
Our Contributions. We define the dynamic temporal graph warping (dtgw) distance as a twofold discrete minimization problem involving computation of an optimal vertex mapping and an optimal "warping path" (see Sect. 3). As a by-product, our approach does not only yield a distance measure but also yields an interpretable mapping between vertices of the two temporal graphs which can, for example, be used for de-anonymization of individuals in social networks.
We show that the dtgw-distance is NP-hard to compute in general (Theorem 1). In contrast, we point out several polynomial-time solvable special cases. This includes the case when either a vertex mapping or a temporal alignment is fixed (Observation 1), the case of deciding whether the dtgw-distance is zero (Theorem 2), and the case when the lifetimes of the two temporal graphs differ only by a constant and the warping path length is restricted (Proposition 3). Moreover, we give a quadratic programming formulation (Sect. 5.2) and propose an efficient and effective heuristic approach (Sect. 5.3).
We empirically evaluate the heuristic 2 in preliminary experiments on real-world temporal social networks (faceto-face contact and spatial proximity networks) against the quadratic program and some simple baseline methods to show its efficiency and solution quality. Moreover, we demonstrate that our concept can successfully be used for deanonymization of real-world temporal social networks and is faster than other existing methods such as DynaMAGNA++ (Vijayan et al. 2017) or HTNE (Zuo et al. 2018) (Sect. 6).
Compared to the conference version (Froese et al. 2019), this version contains the NP-hardness proof of Theorem 1, the proof of Proposition 3, and the quadratic programming formulation (Sect. 5.2). Moreover, we present additional experimental results including a comparison of the heuristic with the quadratic program (Sect. 6.2) and a clustering experiment to demonstrate robustness with respect to noisy data (Sect. 6.3).
Organization. Section 2 contains basic definitions. Section 3 presents our main definition of the dtgw-distance followed by NP-hardness results in Sect. 4 and positive algorithmic results in Sect. 5. Section 6 presents experimental results on some real-world data. We conclude in Sect. 7 with an outlook on future applications and challenges.

Preliminaries
Te m p o r a l G r a p h s . A t e m p o r a l g r a p h G = (V, E 1 , E 2 , … , E T ) consists of a vertex set V and a sequence of T ≥ 1 edge sets E i , each being a set of unordered vertex pairs. By G i = (V, E i ) , we denote the i th layer of G and we call T the lifetime of G . The underlying graph of G is the graph (V, We remark that all definitions and results in this work can easily be extended to labeled temporal graphs (with vertex and/or edge labels).

Page 3 of 16 50
Vertex Mapping. A vertex mapping between two vertex sets V and W is a set M ⊆ V × W containing min(|V|, |W|) tuples such that each x ∈ V ∪ W is contained in at most one tuple of M. We denote the set of all vertex mappings between V and W by M(V, W) . Let V M ⊆ V be the subset of vertices in V that are contained in some tuple of M ( W M ⊆ W is defined analogously). Note that V M = V or W M = W holds since |M| = min(|V|, |W|).
Assignment Problem. Computing optimal vertex mappings between two temporal graphs can be solved via the Assignment Problem which is a fundamental problem in combinatorial optimization. Given two sets A and B of equal size and a cost function c ∶ A × B → ℚ , the goal is to find a bijection ∶ A → B such that ∑ a∈A c(a, (a)) is minimized. It is well known that the Assignment Problem is solvable in O(|A| 3 ) time (Ahuja et al. 1993, Theorem 12.2).
Dynamic Time Warping. The dynamic time warping distance (Sakoe and Chiba 1978) is a distance between time series. It is based on the concept of a warping path. A warping path of order n × m is a set We denote the set of all warping paths of order n × m by P n,m . For two temporal graphs G = (V, E 1 , … , E T ) , H = (W, F 1 , … , F U ) , every order-(T × U) warping path p defines a warping between G and H , that is, a pair (i, j) ∈ p warps layer G i to layer H j .
Parameterized Complexity. We assume the reader to be familiar with basic concepts of computational complexity theory such as NP-completeness (Garey and Johnson 1979). In parameterized complexity theory (Downey and Fellows 2013;Cygan et al. 2015), one considers running times with respect to two dimensions. One dimension is the size of the input instance x, and the other dimension is a parameter k (usually a numerical value). An instance of a parameterized problem is a pair (x, k). The class XP contains all parameterized problems that can be solved in polynomial time for every constant parameter value, that is, in |x| f (k) time for some function f only depending on k.

Dynamic temporal graph warping (DTGW)
In this section, we define our temporal graph distance based on dynamic time warping using a vertex-signaturebased graph distance as cost function. We choose this graph distance for the following reasons. First, in contrast to the NP-hard edit distance, it is polynomial-time computable. Second, it is based on a mapping between the two vertex sets which allows to enforce a consistency over time. This consistency assumption is useful in many temporal network applications where the vertices in both networks correspond to the same set of objects over time. This implicitly allows to identify vertices within the two networks. Third, vertex signatures allow for a high flexibility since they can be chosen arbitrarily in order to incorporate essential information (local or global) for the application at hand (e.g., one might use feature vectors obtained via network embedding).
Graph Distance Based on Vertex Signatures.
The following approach is due to Jouili and Tabbone (2009) For two (static) graphs G = (V, E) and H = (W, F) with vertex signatures f G ∶ V → ℚ k and f H ∶ W → ℚ k and a given vertex mapping M between V and W, we define the cost of M as where G (v) ∈ ℚ is the (predefined) cost of "deleting" vertex v from G since it is not mapped by M to any vertex in the other vertex set. The value G (v) might, for example, depend on the vertex signature of v. Note that "deleting" a vertex does not affect the signatures of other vertices. Note also that one of the last two sums on the right-hand side above is always zero.
The vertex-signature-based distance between G and H is then defined as Depending on the application, one might normalize the distance D by some appropriate factor (typically depending on |V| and |W|; e.g., Jouili and Tabbone (2009) normalized by min(|V|, |W|) −1 ).
Throughout this work, we assume that vertex signature functions f G are computable in polynomial time in the size of G and we assume all metrics d to be polynomial-time computable. In the rest of the paper, we neglect the running times for computing the values of f G and d because we assume that all vertex signatures are precomputed once in polynomial time.
Dynamic Time Warping Distance for Temporal Graphs. We transfer the concept of dynamic time warping to temporal graphs in the following way. Let G = (V, E 1 , … , E T ) and H = (W, F 1 , … , F U ) be two temporal graphs, and let We then define the vertex-signature-based dynamic temporal graph warping distance (dtgw-distance) between G and H as Figure 1 depicts an example illustrating the dtgwdistance of two temporal graphs. Intuitively, the vertex mapping identifies vertices with similar behavior over time and the warping path identifies the time layers with similar vertex behavior. Note that (for T = U) if one fixes p = {(1, 1), (2, 2), … , (T, T)} , then we get a temporal graph distance without time warping (similar to the Euclidean distance).
The following results are easily observed and play a central role for our subsequent algorithms.
O b s e r v a t i o n 1 L e t G = (V, E 1 , … , E T ) a n d H = (W, F 1 , … , F U ) be two temporal graphs with |V| ≤ |W| =∶ n .

For a fixed vertex mapping
Proof i) For a given vertex mapping M, an optimal warping path can be computed by a well-known dynamic program for dynamic time warping in O(T ⋅ U ⋅ n) time (Sakoe and Chiba 1978). Here, O(n) is the time for computing the costs C(G i , H j , M) . Note that faster dynamic time warping algorithms for special cases are known (Abboud et al. 2015;Gold and Sharir 2018;Kuszmaul 2019;Froese et al. 2020).
Then, we need to find a vertex mapping M ∈ M(V � , W) that minimizes ∑ (u,v)∈M (u, v) . Note that M defines a bijection between V ′ and W. Hence, computing M is an Assignment Problem instance solvable in O(n 3 ) time (Ahuja et al. 1993, Theorem 12.2). Computing all (u, v) values can be done in O(n 2 ⋅ |p|) time. ◻ Note that Observation 1 i) implies that if we already know the vertex mapping up to a constant number of vertices, then dtgw can be computed in polynomial time (since we can try out all polynomially many possible vertex mappings). Furthermore, Observation 1 ii) implies that dtgw is polynomial-time computable if the optimal temporal alignment between G and H is known beforehand. In particular, dtgw can be computed in polynomial time if one temporal graph has a constant lifetime or a constant number of vertices since there are only polynomially many possible warping paths or polynomially many vertex mappings.

Corollary 1
The dtgw-distance between two temporal graphs can be computed in polynomial time if at least one of the following applies: 1. The vertex mapping is known up to a constant number of vertices. 2. The warping path is known. Fig. 1 Example of the dtgw-distance between two temporal graphs G and H on four vertices with lifetimes six and five. The vertex coloring indicates an optimal vertex mapping M.

At least one of the temporal graphs has a constant lifetime or a constant number of vertices.
For given vertex signature function and metric, we refer to the decision problem of testing whether two temporal graphs have dynamic temporal graph warping distance at most some given value c by DtgW .

Computational hardness
Even though the dynamic time warping distance and the vertex-signature-based graph distance are both computable in polynomial time, their combined application to temporal graphs yields a distance measure that is generally NP-hard to compute. Intuitively, this is due to the fact that the vertex mapping has to be consistent for all layers. This introduces non-trivial dependencies between the time warping and the vertex mapping which render the problem computationally hard. Indeed, this is not a singular case for temporal graph problems where for many cases the temporal counterparts of problems solvable in polynomial time turn NP-hard; examples include the computation of matchings in graphs (Heeger et al. 2019;Baste et al. 2020;Mertzios et al. 2020), short path computations (Casteigts et al. 2019;Fluschnik et al. 2020), or the computation of separators .

Theorem 1 DtgW is NP -complete for every metric when the vertex signatures are vertex degrees.
Proof DtgW is clearly contained in NP since for a given vertex mapping and warping path (both having polynomial size), one can check in polynomial time whether the dtgw -distance is at most c (also see Observation 1).
To show NP-hardness, we give a polynomial-time manyone reduction from the NP-complete 3 -sAt problem. Let d ∶ ℚ × ℚ → ℚ be any metric, and let = C 1 ∧ … ∧ C m be an instance of 3 -sAt over the variables x 1 , … , x n . Each clause C j is then a disjunction of three literals The idea is to represent each literal by a vertex which can be mapped to either ⊤ (true) or ⊥ (false). We then build, for each clause, a clause box gadget consisting of three consecutive layers. The choice of a warping path will then, for each clause, implicitly select one of its literals, and the costs caused by each clause box will attain their minimum value if and only if that particular literal is mapped to ⊤. Now, a detailed description of the reduction follows. Let D and D ′ be two copies of the graph ∐ 22m i=1 K 2 (consisting of 22m disjoint edges), where for each vertex v ∈ V(D) we denote its copy in V(D � ) by v ′ . We construct two temporal graphs G and H . Their vertex sets each contain 2n + 47m + 8 vertices as follows: Both temporal graphs have 2n + 26m layers defined as follows. For each i ∈ [n] , we set For j ∈ [m] , we set and Finally, for j ∈ [22m] , we set We call the layers containing |E(D)| edges separation layers. Furthermore, for each j ∈ [m] we say that the layers 2n + 4j − 3 , 2n + 4j − 2 , and 2n + 4j − 1 form the clause block corresponding to C j (see Fig. 2 for an example). Let c ∶= 42m ⋅ d(0, 1).
We claim that dtgw (G, H) ≤ c if and only if has a satisfying truth assignment. "⇐ " : G i v e n a s a t i s f y i n g a s s i g n m e n t ∶ {x 1 , … , x n } → { , } of , we define the following vertex mapping: Only relevant vertices are shown in each layer The three possible warpings between layers of a clause block. Each edge is labeled with the minimal cost it causes under the assumption that the set To construct a warping path, we begin by defining, for each j ∈ [m] , the following three sub-paths (see also Fig. 3): is true. We then build the warping path p as the union of all k j j , using the trivial warping path for all remaining layers: It is then not difficult to calculate that each clause block adds cost of exactly 42 ⋅ d(0, 1) and there are no other costs. Thus, . Note that any non-separation layer contains at most eight edges. So if p warps any separation layer to any non-separation layer, then the resulting cost would be at least (44m − 16) ⋅ d(0, 1) > c . Thus, we may assume that every separation layer i of G is only warped to layer i of H and vice versa. Since the last 22m layers of each temporal graph are all identical and M and p are chosen to have minimal cost, we can conclude that } , then the 22m layers 2n + 4m + 1, … , 2n + 4m + 22m each would cause cost of at least 2 ⋅ d(0, 1) , thus exceeding c in total. Hence, M has to contain a bijection from }. Now, consider the clause block corresponding to C j = 1 j ∨ 2 j ∨ 3 j . From the arguments above, it follows that G 2n+4j−3 and G 2n+4j−1 are warped to H 2n+4j−3 and H 2n+4j−1 , respectively. This already costs 32 ⋅ d(0, 1) . We distinguish three cases (corresponding to 1 j through 3 j above): 1. G 2n+4j−2 is warped to H 2n+4j−3 . This causes costs of at least 2 ⋅ d(0, 1) . Then, H 2n+4j−2 must be warped to G 2n+4j−1 or p would not have minimal cost. Thus, there are additional costs of at least 8 ⋅ d(0, 1) . This is the situation illustrated in Fig. 3a.
In summary, the costs contributed by each clause block are at least 42 ⋅ d(0, 1) . Hence, to meet the bound of c, all layers outside of clause blocks must not cause any additional cost.
, the clause block corresponding to C j must have cost of exactly 42 ⋅ d(0, 1) . If we are in Case (1) as above, then this is only possible if M maps each degree-1 vertex of G 2n+4j−2 to some degree-1 vertex of H 2n+4j−3 . Thus, 2 j , ⊤ (j,2) ∈ M . Otherwise, if we are in Case (2), respectively, Case (3), then analogous arguments yield that 1 j , ⊤ (j,1) ∈ M , respectively, 3 j , ⊤ (j,3) ∈ M . Hence, in any case there is some Consequently, is a satisfying assignment for . ◻ Let us take a closer look at the reduction in the proof of Theorem 1. Note that the corresponding optimal warping path is always close to the diagonal (that is, |i − j| ≤ 1 holds for every pair (i, j)). Hence, it lies within the so-called Sakoe-Chiba band (Sakoe and Chiba 1978) of width one. Moreover, the maximum degree in each layer is one. Finally, the number of vertices and the number of layers of both temporal graphs and the target cost c are all upper-bounded linearly in the size of the 3 -sAt formula, which allows to conclude a running time lower bound based on the Exponential Time Hypothesis 3  [together with the Sparsification Lemma ]. These observations are summarized in the following corollary. (Recall that V, W are the vertex sets and T, U are the lifetimes of G , H.) Corollary 2 DtgW is NP -complete for every metric and vertex degrees as vertex signatures even when the maximum degree of each layer is one and the warping path is Furthermore, if the dtgw-distance is normalized (e.g., divided by the number of vertices), then we obtain NPhardness for a constant value of c (by the reduction in the proof of Theorem 1).

Corollary 3 DtgW with normalized vertex-signaturebased distance is NP-complete for a constant value of c.
Due to the proven worst-case hardness of DtgW , there is little hope to solve the problem efficiently in general. Nevertheless, in the next section we provide some methods to cope with this intractability.

Algorithms
In this section, we first point out polynomial-time solvable special cases. Then, we develop a mathematical programming formulation as well as a heuristic approach to approximate the dtgw-distance in practice.

Exact polynomial-time algorithms for special cases
Our first algorithmic result is to show that one can determine in polynomial time whether two temporal graphs with the same number of vertices have dtgw-distance zero. This basic case occurs when checking for duplicates within a data set. In contrast, determining whether two (static) graphs have graph edit distance zero is not known to be polynomial-time solvable (as this is equivalent to the famous grAPh isomor-Phism problem). Proof We will show that for distance zero, an optimal warping path can easily be determined. Polynomial-time solvability then follows from Observation 1. Let G = (V, E 1 , … , E T ) and H = (W, F 1 , … , F U ) be two t e m p o r a l g r a p h s w i t h V =∶ {v 1 , … , v n } and W =∶ {w 1 , … , w n } . For each i ∈ [T] , we define the i th layer signature of G as f (

Theorem 2 Let
). Assum-ing dtgw (G, H) = 0 , it follows that there exists a vertex mapping M ⊆ V × W and a warping path p ∈ P T,U such that holds for every (i, j) ∈ p . Since d is a metric, this implies that and layer i is warped to layer j and layer i ′ is warped to layer j ′ , then f (H j ) ≠ f (H j � ) since otherwise the cost will not be zero. By the definition of a warping path, it follows that the layers 1, … , i 1 of G can only be warped to layers 1, … , j 1 of H and the layers i 1 + 1, … , i 2 of G can only be warped to layers j 1 + 1, … , j 2 of H and so on. Note that this is only possible if q = r . If this is the case, then we can assume that the warping path p has the following form: By Observation 1, we can now check whether there exists a vertex mapping that yields distance zero for the warping path p in O(n 2 ⋅ (T + U) + n 3 ) time. Computing p can be done in O(n(T + U)) time. ◻ We remark that if the vertex signatures and the metric satisfy the property that every pair of different vertex signatures has distance at least for some constant > 0 , then DtgW parameterized by the resulting cost c is in XP . For example, this is the case when the vertex signatures contain only integers and d is any p -norm (for p ≥ 1 ). Then, every pair of different signatures has distance at least = 1 . The idea of the algorithm is to "guess" the tuples of a warping path which cause nonzero cost (at most c∕ many) and to check whether it is possible to complete the warping path without further costs. The latter can be done in polynomial time using similar arguments as for the case c = 0 (Theorem 2).

Corollary 4 DtgW is in XP with respect to c if the vertex signatures have a constant minimum pairwise distance
In contrast, if the dtgw-distance is normalized, then the differences between vertex signatures can be arbitrarily small, in which case DtgW is NP-hard for constant c (Corollary 3).
To overcome this hardness, in the following, we consider special cases based on parameters regarding the warping path length. We assume that the lifetimes of the inputs differ by at most a constant, that is, T = U + t for some t ≥ 0 (which might often be the case in practice). Note that, by definition, every warping path of order T × U has length at least T. We define the parameter to be the difference between the warping path length and the lower bound T, that is, we consider only order-(T × U) warping paths of length at most T + . (In practice, overly long warping paths might be considered unnatural.) We prove that DtgW is in XP with respect to the combined parameter ( , t).

Theorem 3 For all vertex signatures and all metrics,
DtgW is solvable in time if n = max(|V|, |W|) , T = U + t, and the warping paths have length at most T + .
t e m p o r a l g r a p h s , a n d let p = {p 1 = (i 1 , j 1 ), … , p L = (i L , j L )} ∈ P T,U be a warping path. The warping path p contains L − 1 steps p +1 − p = (i +1 − i , j +1 − j ) ∈ {(1, 0), (0, 1), (1, 1)} for 1 ≤ < L . We call a step horizontal if p +1 − p = (1, 0) , and we call it vertical if p +1 − p = (0, 1) , and otherwise we call it diagonal. Let ≤ U − 1 denote the number of vertical steps in p. Then, p contains also + t horizontal and U − 1 − d i a g o n a l s t e p s , t h a t i s , O(max(|V|, |W|) 2 ⋅ (T + ) + max(|V|, |W|) 3 ) Note that Proposition 3 implies polynomial-time solvability of DtgW if t and are constants. For unbounded t, however, we conjecture that DtgW is NP-hard even if the warping paths are restricted to have length max(T, U) , which is the minimum possible length (that is, = 0).

Quadratic programming
We give a formalization of DtgW as a quadratic minimization problem with linear constraints (QP). This is NP-hard to solve in general but can be used to solve small instances exactly with the state-of-the-art QP-solvers such as Gurobi 4 .
We use the following variables: Computing dtgw (G, H) is the following quadratic 5 minimization problem.
4 www.gurob i.com 5 It is also possible to convert our formulation into a linear problem by introducing further variables and constraints for replacing the product w s,t ⋅ m i,j in the objective function. However, in our experiments we found that the quadratic formulation can be solved faster.

Page 10 of 16
The constraints 1b and 1c ensure that the vertex mapping variables define a correct vertex mapping, that is, every vertex is mapped to exactly one other vertex (or is deleted). Constraints 1d to 1g ensure that the warping variables define a valid warping path. Here, the constraints 1e to 1g imply that if the warping path contains a pair (s, t), then it also contains at least one of the pairs (s + 1, t) , (s, t + 1) , or (s + 1, t + 1). (Since the objective is minimized, any solution will actually select only one of these pairs.) The number of variables is in O(|V| ⋅ |W| + T ⋅ U) , and the number of constraints is in O(|V| + |W| + T ⋅ U).

Heuristic approaches
In this section, we present a heuristic to compute the dtgwdistance, which typically yields good (not necessarily optimal) solutions in practice.
The approach is to simply start with an arbitrary initial vertex mapping (or warping path) and to compute an optimal warping path (vertex mapping) based on Observation 1 in polynomial time. This process is then repeated by alternating between optimal warping path and optimal vertex mapping computation until the solution converges to a local minimum (or some other criterion is reached).
Note that it is a convenient feature of our heuristic to be able to stop the process after any number of iterations to obtain some approximate solution (a so-called anytime algorithm). It is further possible to incorporate prior knowledge, for example, by fixing the mapping for some vertices. Note also that convergence is guaranteed since we decrease the objective in each alternation and the search space is finite. We propose several initialization options.
Initial Warping Path. A first idea for initialization is to choose a shortest warping path (that is, of length max(T, U) ). Note that for T ≠ U several such paths exist. Without further knowledge about the instances, choosing a path within the Sakoe-Chiba band of small width is a reasonable default. This initialization is very simple and only requires O(T + U) time.
Another idea is to compute a warping path using D(G i , H j ) as a cost for warping layer i to layer j. This is of course an optimistic estimate since it allows to use a different vertex mapping for each pair of layers. Then, a vertex mapping can be computed by Observation 1. This initialization takes O(T ⋅ U ⋅ n 3 ) time where n ∶= max(|V|, |W|).
Initial Vertex Mapping. The idea is to compute a vertex mapping by solving an Assignment Problem instance for approximate costs. Let (u, v) be some approximate cost for mapping vertex u ∈ V to v ∈ W . For example, one could use the following estimations: The first option * estimates the cost of mapping u to v over all possible warpings between any two layers. (This is usually more than any warping path will incur.) The definition of opt only considers for each layer of the first temporal graph the minimal cost over all layers of the other temporal graph. (This estimate might be too low.) Both of these options require O(T ⋅ U ⋅ |V| ⋅ |W|) time.
Based on the estimated costs, one computes a vertex mapping by solving an Assignment Problem instance and then computes an optimal warping path for this vertex mapping based on Observation 1.
The running time of one iteration (that is, computing a vertex mapping and an optimal warping path) is O (T + U) ⋅ n 2 + n 3 + T ⋅ U ⋅ n . While the number of iterations might depend on the choice of initialization, in our experiments the heuristic always converged after very few iterations. Regarding the solution quality, while it is possible to construct adversarial examples where the heuristic performs poorly, our experiments in Sect. 6.2 indicate that it performs well in practice.

Experiments
We conducted several experiments 6 to demonstrate the merit of our dtgw-distance in applications and to evaluate the performance of the alternating minimization heuristic (Am) we described in Sect. 5.3. For computations, we used a 4.0 GHz i7-6700K processor (single-threaded).

Data sets
We used three data sets from the SocioPatterns project (Génois and Barrat 2018). Each of these consists of two temporal social networks, both recorded simultaneously with the same individuals as vertices. The first network is a faceto-face contact network, whereas the second one is a copresence network where edges represent spatial proximity. All six networks have a temporal resolution of 20 seconds.
For our experiments, only the first day of each network was used and vertices without any edges were discarded. The three different data sets were recorded at a primary school ("LyonSchool," 237 vertices, 1700 layers), a scientific conference ("SFHH," 403 vertices, 2300 layers), and a workplace ("InVS15," 180 vertices, 2100 layers).

Comparison of heuristic and exact solutions
We compared the solutions of our Am heuristic under different initialization schemes against the optimal solutions obtained from the QP formulation. Due to long running times of the QP-solver, we were restricted to very small temporal networks. We randomly selected 10 children from class 1A of the primary school face-to-face network and extracted 225 consecutive layers (during a high contact period) which we split into 15 temporal subnetworks with 15 consecutive layers each. We used vertex degrees as signatures (with absolute value metric) and computed all pairwise dtgw-distances between these 15 networks with the following algorithms: • QP: exact QP-solver (Gurobi 8.0.1), • Am * : Am with * initialization, • Am opt : Am with opt initialization, • Am swp : Am with shortest warping path initialization, • Am owp : Am with optimistic warping path initialization.
We implemented the Am heuristic in Python, using a C++ implementation 7 of the Jonker-Volgenant algorithm (Jonker and Volgenant 1987) to solve the Assignment Problem. Figure 4 shows for each initialization variant the estimated cumulative distribution function (ecdf) of the error percentage = 100 ⋅ (d AM − d QP )∕d QP , where d AM is the approximate dtgw-distance obtained by an Am heuristic and d QP is the exact dtgw-distance obtained by the QP-solver. A point ( , P) on an ecdf-curve of an Am heuristic means that the error percentage of Am is at most with estimated probability P.
All Am variants found the correct solution for a majority of samples (P 0 > 0.5 ). The average error percentages are rather small and vary between 3.0 by Am owp and 5.5 by Am opt . The Am owp heuristic performed the best, having P 0 ≈ 0.71 and maximum error percentage max ≈ 36.4 . These findings indicate that for small instances the approximations of the four heuristics are close to the optimal solution on average but may fail considerably in some cases with up to a maximum error of 63.6% . We remark that based on our experimental experience the relative error becomes smaller on larger instances.
Regarding running times, the Am heuristic took less than 0.01 seconds per instance, usually converging after at most three iterations (independently of the chosen initialization). In comparison, the QP was slower by a factor of more than 10 000, requiring 8 minutes on average (median 2 minutes) with some instances approaching 2 hours.

Sensitivity of DTGW to noise
The goal of this experiment was to assess how sensitive the dtgw-distance is to noise, that is, how well can original data be reconstructed from noisy data. We compared our dtgwdistance approach to the following two baseline methods.  • Non-consistent: Instead of using one consistent vertex mapping for all layers, one can allow a different mapping for each pair of layers. Note that the resulting distance can be computed in O(T ⋅ U ⋅ n 3 ) time, thus being faster than an exact computation of the dtgw-distance but much slower than a single iteration of the Am heuristic. • Non-temporal: A naive approach is to ignore the time information and solely compute an optimal vertex mapping between the underlying graphs. This requires O(n 3 ) time plus the (usually linear) time to build the underlying graphs.
We used the primary school face-to-face network from which we extracted five reference temporal networks representing the contacts between children of the same grade, each containing 45-50 vertices and 3100 layers. For each of the five reference networks, we generated nine noisy copies as follows: We used the Am heuristic to approximate the pairwise dtgw-distances (using degrees as vertex signatures) between all reference and noisy temporal networks. In all of these instances, shortest warping path initialization (which is fastest) was used since preliminary tests showed that the other initializations produce very similar results. Figure 5 shows the dendrogram obtained by hierarchical clustering using complete linkage of the approximated pairwise dtgw-distances. Both the dtgw-distance and the nonconsistent baseline were able to partition the instances into five clusters, each of which consists of a reference network and its nine noisy copies. Hence, they successfully recovered the original reference networks from noise. However, the clusters produced by the dtgw-distance are more compact than the ones of the non-consistent baseline. In contrast, the non-temporal baseline was not able to separate the graphs of different grades.
In all instances, the heuristic converged within at most six iterations, taking less than 15 seconds. In comparison, the non-consistent baseline required 4 minutes on average for each instance, while the non-temporal baseline was the fastest (below 1 second).

De-anonymization
Besides measuring a distance between temporal graphs, the dtgw-distance additionally provides a mapping between the vertex sets which implicitly allows to identify vertices. This allows the de-anonymization (Narayanan and Shmatikov 2009) of temporal social networks. Since the data sets used contain the original mapping between the vertex sets, we can employ this as an easy benchmark for the accuracy of the dtgw-distance.
We used the Am heuristic (with shortest warping path initialization) to compute the dtgw-distance (with degrees as vertex signatures 9 ) on the three data sets mentioned in Sect. 6.1. We counted how many vertices were correctly reidentified (that is, mapped to their copies) in the resulting vertex mapping. We compared our results to the following alternative algorithms found in the literature: • DynaMAGNA++ (Vijayan et al. 2017): A search-based evolutionary algorithm computing a vertex mapping that maximizes edge conservation and node conservation over time. • Temporal network embedding (Zuo et al. 2018): Hawkes process-based temporal network embedding (HTNE) computes a low-dimensional embedding of the vertices of a temporal network. From this, we computed a vertex mapping minimizing the Euclidean distances between the vertex feature vectors. • Fixed dtgw: Note that our dtgw-distance allows to fix the warping path beforehand (Observation 1 ii)). Since in our case each pair of temporal graphs was recorded using synchronized clocks, it is natural to use a fixed warping path that aligns layer i of the first graph with layer i of the second graph.
To simulate a situation in which the temporal graphs represent processes which do not run synchronously in time, we created two modified versions of each of the data sets. In the first one, called "shifted," all events of the first graph were delayed by three minutes. In the second version, called "randomized," each layer of each of the graphs was randomly and independently replaced by X layers where X ∈ {1, 2, … } is a random variable with ℙ(X ≥ x) = x −3 . Since we pretend that the nature of these modifications is unknown to the tested algorithms, dtgw with fixed warping path is not applicable to these variants. For DynaMAGNA++, we used a population of size 15 000 and a maximum of 10 000 generations. With HTNE, we computed 128-dimensional vertex embeddings using a batch of size 10 000, a learning rate of 0.1, a history length of 2, and 5 negative samples. Unlike dtgw, both methods utilized all four processor cores.
The results and running times are listed in Table 1. Most notably, the re-identification rate of HTNE was poor on all data sets, suggesting that these embeddings are ill-suited for comparing vertices taken from different networks. Furthermore, all methods failed to re-identify any significant number of vertices on the primary school data set. This might be explained by the fact that the co-presence network is very different from the face-to-face contacts due to a low spatial resolution (as was also noted by Génois and Barrat (2018)).
The overall performance was much better on the other two data sets, especially on the conference data where up to 90% of participants could be re-identified, whereas on the workplace data set the best result was 51%. Unsurprisingly, fixing the correct layer alignment on the unmodified graphs sped up the dtgw computation significantly while also yielding slightly better results. On these instances, dtgw performed comparably to DynaMAGNA++, being better in one case and worse in the other, although requiring much less computational effort. In contrast, on the shifted and randomized data sets dtgw always achieved the best results. (Notably, the re-identification performance using dtgw did not decrease on shifted data.) In all cases, the Am heuristic converged after at most six iterations and DynaMAGNA++ converged within 2 000 generations.

Conclusion
We introduced a new proximity measure for comparing temporal graphs by transferring dynamic time warping from time series to temporal graphs. This yields a challenging computational problem for which we proposed exact algorithms and a heuristic approach to solve it. While exact solutions can only be computed for very small instances, we empirically showed that our heuristic runs fast in practice and yields good approximations of optimal solutions. In our experiments, it was also capable of de-anonymizing social networks.
Our work opens several directions for future research. We believe that the dtgw-distance is a promising tool, for example, in biology and chemistry. Processes like epidemic disease spreading or chemical reactions can naturally be viewed as temporal graphs where the vertices represent individuals or (macro) molecules. [Unfortunately, we could not test this, as there is still a lack of openly available temporal molecular data (Vijayan et al. 2017).] Since the exact timescales of these processes often vary, the ability of dynamic time warping to compensate for such differences would be especially helpful in this context. One might also use the dtgw-distance to understand the learning process of neural networks. The training phases of neural networks yield temporal networks which can be analyzed to gain insight into how different conditions influence the learning process. Another potential application is analyzing team sports data via temporal graphs to reveal similar strategies or roles of individual players. Depending on the application domain, there is a wide range of possibilities to test the performance when using different vertex signatures or even other graph distances.
Besides experimenting with various application domains and further definition variants, already the proven computational worst-case hardness of DtgW may trigger further algorithmic research. A concrete open question is whether DtgW is in XP (or even fixed-parameter tractable) with respect to , when the warping path length is restricted to be at most max(T, U) + where T, U are the respective lifetimes. It is also interesting to further study the influence of graph-specific parameters in the spirit of a multivariate complexity analysis (Niedermeier 2010;Fellows et al. 2013).
provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.