Year: 2020

Abstract:
Representation learning models for graphs are a successful family of techniques that project nodes into feature spaces that can be exploited by other machine learning algorithms. Since many real-world networks are inherently dynamic, with interactions among nodes changing over time, these techniques can be defined both for static and for time-varying graphs. Here, we show how the skip-gram embedding approach can be used to perform implicit tensor factorization on different tensor representations of time-varying graphs. We show that higher-order skip-gram with negative sampling (HOSGNS) is able to disentangle the role of nodes and time, with a small fraction of the number of parameters needed by other approaches. We empirically evaluate our approach using time-resolved face-to-face proximity data, showing that the learned representations outperform state-of-the-art methods when used to solve downstream tasks such as network reconstruction. Good performance on predicting the outcome of dynamical processes such as disease spreading shows the potential of this new method to estimate contagion risk, providing early risk awareness based on contact tracing data. The source code and data are publicly available at [link to anonymized repository].

heading: 1 INTRODUCTION
text: A great variety of natural and artificial systems can be represented as networks of elementary structural entities coupled by relations between them. The abstraction of such systems as networks helps us understand, predict and optimize their behaviour (Newman, 2003; Albert & Barabási, 2002). In this sense, node and graph embeddings have been established as standard feature representations in many learning tasks (Cai et al., 2018; Goyal & Ferrara, 2018). Node embedding methods map nodes into low-dimensional vectors that can be used to solve downstream tasks such as edge prediction, network reconstruction and node classification.
Node embeddings have proven successful in achieving low-dimensional encoding of static network structures, but many real-world networks are inherently dynamic (Holme & Saramäki, 2012). Timeresolved networks are also the support of important dynamical processes, such as epidemic or rumor spreading, cascading failures, consensus formation, etc. (Barrat et al., 2008). Time-resolved node embeddings have been shown to yield improved performance for predicting the outcome of dynamical processes over networks, such as information diffusion and disease spreading (Sato et al., 2019), providing estimation of infection and contagion risk when used with contact tracing data.
Since we expect having more data on proximity networks being used for contact tracing and as proxies for epidemic risk (Alsdurf et al., 2020), learning meaningful representations of time-resolved proximity networks can be of extreme importance when facing events such as epidemic outbreaks (Kapoor et al., 2020; Gao et al., 2020). The manual and automatic collection of time-resolved proximity graphs for contact tracing purposes presents an opportunity for quick identification of possible infection clusters and infection chains. Even before the COVID-19 pandemic, the use of wearable proximity sensors for collecting time-resolved proximity networks has been largely discussed in the literature and many approaches have been used to describe patterns of activity and community structure, and to study spreading patterns of infectious diseases (Sapienza et al., 2015; Gauvin et al., 2014; Génois et al., 2015).
Here we propose a representation learning model that performs implicit tensor factorization on different higher-order representations of time-varying graphs. The main contributions are as follows: Given that the skip-gram embedding approach implicitly performs a factorization of the shifted pointwise mutual information matrix (PMI) (Levy & Goldberg, 2014), we generalize it to perform implicit factorization of a shifted PMI tensor. We then define the steps to achieve this factorization using higher-order skip-gram with negative sampling (HOSGNS) optimization. We show how to apply 3rd-order and 4th-order SGNS on different higher-order representations of time-varying graphs. Finally, we show that time-varying graph representations learned via HOSGNS outperform stateof-the-art methods when used to solve downstream tasks, even using a fraction of the number of embedding parameters.
We report the results of learning embeddings on empirical time-resolved face-to-face proximity data and using such representations as predictors for solving two different tasks: network reconstruction and predicting the outcomes of a SIR spreading process over the time-varying graph. We compare these results with state-of-the art methods for time-varying graph representation learning.

heading: 2 PRELIMINARIES AND RELATED WORK
text: Skip-gram representation learning. The skip-gram model was designed to compute word embeddings in WORD2VEC (Mikolov et al., 2013), and afterwards extended to graph node embeddings (Perozzi et al., 2014; Tang et al., 2015; Grover & Leskovec, 2016). Levy & Goldberg (2014) established the relation between skip-gram trained with negative sampling (SGNS) and traditional low-rank approximation methods (Kolda & Bader, 2009; Anandkumar et al., 2014), showing the equivalence of SGNS optimization to factorizing a shifted PMI matrix (Church & Hanks, 1990). This equivalence was later retrieved from diverse assumptions (Assylbekov & Takhanov, 2019; Allen et al., 2019; Melamud & Goldberger, 2017; Arora et al., 2016; Li et al., 2015), and exploited to compute closed form expressions approximated in different graph embedding models (Qiu et al., 2018). In this work, we refer to the shifted PMI matrix also as SPMIκ = PMI− log κ, where κ is the number of negative samples.
Random walk based graph embeddings. Given an undirected, weighted and connected graph G = (V, E) with nodes i, j ∈ V , edges (i, j) ∈ E and adjacency matrix A, graph embedding methods are unsupervised models designed to map nodes into dense d-dimensional representations (d |V|) (Hamilton et al., 2017). A well known family of approaches based on the skip-gram model consists in sampling random walks from the graph and processing node sequences as textual sentences. In DEEPWALK (Perozzi et al., 2014) and NODE2VEC (Grover & Leskovec, 2016), the skip-gram model is used to obtain node embeddings from co-occurrences in random walk realizations. Although the original implementation of DEEPWALK uses hierarchical softmax to compute embeddings, we will refer to the SGNS formulation given by Qiu et al. (2018). Since SGNS can be interpreted as a factorization of the word-context PMI matrix (Levy & Goldberg, 2014), the asymptotic form of the PMI matrix implicitly decomposed in DEEPWALK can be derived (Qiu et al., 2018). Given the 1-step transition matrix P = D−1A, where D = diag(d1, . . . , d|V|) and di = ∑ j∈V Aij is the (weighted) node degree, the expected PMI for a node-context pair (i, j) occurring in a T -sized window is:
E[ PMIDEEPWALK(i, j) | T ] = log ( 1 2T ∑T r=1 [p ∗(i)(Pr)ij + p ∗(j)(Pr)ji]
p∗(i) p∗(j)
) (2.1)
where p∗(i) = divol(G) is the unique stationary distribution for random walks (Masuda et al., 2017) and vol(G) = ∑i,j∈V Aij . We will use this expression in Section 3.2 to build PMI tensors from higher-order graph representations.
Time-varying graphs and their algebraic representations. Time-varying graphs (Holme & Saramäki, 2012) are defined as triples H = (V, E , T ) , i.e. collections of events (i, j, k) ∈ E , representing undirected pairwise relations among nodes at discrete times (i, j ∈ V , k ∈ T ). H can be seen as a temporal sequence of static graphs {G(k)}k∈T , each of those with adjacency matrix A(k) such that A(k)ij = ω(i, j, k) ∈ R is the weight of the event (i, j, k) ∈ E . We can concatenate the list of time-stamped snapshots [A(1), . . . ,A(|T |)] to obtain a single 3rd-order tensor
Astat(H) ∈ R|V|×|V|×|T | which characterize the evolution of the graph over time. This representation has been used to discover latent community structures of temporal graphs (Gauvin et al., 2014) and to perform temporal link prediction (Dunlavy et al., 2011). Indeed, beyond the above stacked graph representation, more exhaustive representations are possible. In particular, the multi-layer approach (De Domenico et al., 2013) allows to map the topology of a time-varying graph H into a static network GH = (VH, EH) (the supra-adjacency graph) such that vertices of GH correspond to pairs (i, k) ≡ i(k) ∈ V × T of the original time-dependent network. This representation can be stored in a 4th-order tensor Adyn(H) ∈ R|V|×|V|×|T |×|T | equivalent, up to an opportune reshaping, to the adjacency matrix A(GH) ∈ R|V||T |×|V||T | associated to GH. Multi-layer representations for time-varying networks have been used to study time-dependent centrality measures (Taylor et al., 2019) and properties of spreading processes (Valdano et al., 2015).
Time-varying graph representation learning. Given a time-varying graph H = (V, E , T ), we define as temporal network embedding a model that learns from data, implicitly or explicitly, a mapping function: f : (v, t) ∈ V × T 7→ v(t) ∈ Rd (2.2) which project time-stamped nodes into a latent low-rank vector space that encodes structural and temporal properties of the original evolving graph. Many existing methods learn node representations from sequences of static snapshots through incremental updates in a streaming scenario: deep autoencoders (Goyal et al., 2017), SVD (Zhang et al., 2018), skip-gram (Du et al., 2018) and random walk sampling (Béres et al., 2019; Mahdavi et al., 2018; Yu et al., 2018). Another class of models learn dynamic node representations by recurrent/attention mechanisms (Goyal et al., 2020; Li et al., 2018; Sankar et al., 2020; Xu et al., 2020) or by imposing temporal stability among adjacent time intervals (Zhou et al., 2018; Zhu et al., 2016). DYANE (Sato et al., 2019) and WEG2VEC (Torricelli et al., 2020) project the dynamic graph structure into a static graph, in order to compute embeddings with WORD2VEC. Closely related to these ones are Nguyen et al. (2018) and Zhan et al. (2020), which learn node vectors according to time-respecting random walks or spreading trajectory paths. Moreover, Kumar et al. (2019) proposed an embedding framework for user-item temporal interactions, and Malik et al. (2020) suggested a tensor-based convolutional architecture for dynamic graphs.
Methods that perform well for predicting outcomes of spreading processes make use of timerespecting supra-adjacency representations such as the one proposed by Valdano et al. (2015). In this representation, random itineraries correspond to temporal paths of the original time-varying graph. The supra-adjacency representation GH that we refer in Section 3.2, also used in DYANE, with adjacency matrix A(GH), is defined by two rules:
1. For each event (i, j, t0), if i is also active at time t1 > t0 and in no other time-stamp between the two, we add a cross-coupling edge between supra-adjacency nodes j(t0) and i(t1). In addition, if the next interaction of j with other nodes happens at t2 > t0, we add an edge between i(t0) and j(t2). The weights of such edges are set to ω(i, j, t0).
2. For every case as described above, we also add self-coupling edges (i(t0), i(t1)) and (j(t0), j(t2)), with weights set to 1.
Figure 1 shows the differences between a time-varying graph and its time-aware supra-adjacency representation, according to the formulation described above. DYANE computes, given a node i ∈ V , one vector representation for each time-stamped node i(t) ∈ V(T ) = {(i, t) ∈ V×T : ∃ (i, j, t) ∈ E} of this supra-adjacency representation. Similar models that learn time-resolved node representations
require a quantity O(|V| · |T |) of embedding parameters to represent the system in the latent space. Compared with these methods, our approach requires a quantity O(|V|+ |T |) of embedding parameters for disentangled node and time representations.

heading: 3 PROPOSED METHOD
text: Given a time-varying graphH = (V, E , T ), we propose a representation learning method that learns disentangled representations for nodes and time slices. More formally, we learn a function:
f∗ : (v, t) ∈ V × T 7→ v, t ∈ Rd
through a number of parameters proportional to O(|V|+ |T |). This embedding representation can then be reconciled with the definition in Eq. (2.2) by combining v and t in a single v(t) representation using any combination function c : (v, t) ∈ Rd × Rd 7→ v(t) ∈ Rd. Starting from the existing skip-gram framework for node embeddings, we propose a higher-order generalization of skip-gram with negative sampling (HOSGNS) applied to time-varying graphs. We show that it allows to implicitly factorize higher-order relations that characterize tensor representations of time-varying graphs, in the same way that the classical SGNS decomposes dyadic relations associated to a static graph. Similar approaches have been applied in NLP for dynamic word embeddings (Rudolph & Blei, 2018), and higher-order extensions of the skip-gram model have been proposed to learn context-dependent (Liu et al., 2015) and syntactic-aware (Cotterell et al., 2017) word representations. Moreover tensor factorization techniques have been applied to include the temporal dimension in recommender systems (Xiong et al., 2010; Wu et al., 2019), knowledge graphs (Lacroix et al., 2020; Ma et al., 2019) and face-to-face contact networks (Sapienza et al., 2015; Gauvin et al., 2014). But this work is the first to merge SGNS with tensor factorization, and then apply it to learn time-varying graph embeddings.

heading: 3.1 HIGHER-ORDER SKIP-GRAM WITH NEGATIVE SAMPLING AS IMPLICIT TENSOR FACTORIZATION
text: Here we address the problem of generalizing SGNS to learn embedding representations from higherorder co-occurrences. We analyze here the 3rd-order case, giving the description of the general N -order case in the Supplementary Information. Later in this work we will focus 3rd and 4th order representations since these are the most interesting for time-varying graphs.
We consider a set of training samples D = {(i, j, k), i ∈ W, j ∈ C, k ∈ T } obtained by collecting co-occurrences among elements from three setsW , C and T . While SGNS is limited to pairs of nodecontext (i, j), hereD is constructed with three (or more) variables, e.g. sampling random walks over a higher-order data structure. We denote as #(i, j, k) the number of times the triple (i, j, k) appears in D. Similarly we use #i =∑j,k#(i, j, k), #j =∑i,k#(i, j, k) and #k =∑i,j #(i, j, k) as the number of times each distinct element occurs in D, with relative frequencies PD(i, j, k) = #(i,j,k)|D| , PD(i) = #i |D| , PD(j) = #j |D| and PD(k) = #k |D| .
Optimization is performed as a binary classification task, where the objective is to discern occurrences actually coming from D from random occurrences. We define the likelihood for a single observation (i, j, k) by applying a sigmoid (σ(x) = (1 + e−x)−1) to the higher-order inner product [[·]] of corresponding d-dimensional representations:
P [ (i, j, k) ∈ D | wi, cj , tk ] = σ ( [[wi, cj , tk]] ) ≡ σ ( ∑d r=1 WirCjrTkr ) (3.1)
where embedding vectors wi, cj , tk ∈ Rd are respectively rows of W ∈ R|W|×d, C ∈ R|C|×d and T ∈ R|T |×d. In the 4th-order case we will also have a fourth embedding matrix S ∈ R|S|×d related to a fourth set S. For negative sampling we fix an observed (i, j, k) ∈ D and independently sample jN and kN to generate κ negative examples (i, jN , kN ). In this way, for a single occurrence (i, j, k) ∈ D, the expected contribution to the loss is:
`(i, j, k) = log σ ( [[wi, cj , tk]] ) + κ · E
jN ,kN∼PN
[ log σ ( − [[wi, cjN , tkN ]] )] (3.2)
where the noise distribution is the product of independent marginal probabilities PN (j, k) = PD(j) · PD(k). Thus the global objective is the sum of all the quantities of Eq. (3.2) weighted with the corresponding relative frequency PD(i, j, k). The full loss function can be expressed as:
L = − |W|∑ i=1 |C|∑ j=1 |T |∑ k=1 [ PD(i, j, k) log σ ( [[wi, cj , tk]] ) + κ PN (i, j, k) log σ ( − [[wi, cj , tk]] )] (3.3)
In Supplementary Information we show the formal steps to obtain Eq. (3.3) for the N -order case and that it can be optimized with respect to the embedding parameters, satisfying the low-rank tensor approximation of the multivariate shifted PMI tensor into factor matrices W,C,T:∑d
r=1 WirCjrTkr ≈ log
( PD(i, j, k)
PN (i, j, k)
) − log κ ≡ SPMIκ(i, j, k) (3.4)

heading: 3.2 TIME-VARYING GRAPH EMBEDDING VIA HOSGNS
text: While a static graph G = (V, E) is uniquely represented by an adjacency matrix A(G) ∈ R|V|×|V|, a time-varying graphH = (V, E , T ) admits diverse possible higher-order adjacency relations (Section 2). Starting from these higher-order relations, we can either use them directly or use random walk realizations to build a dataset of higher-order co-occurrences. In the same spirit that random walk realizations give place to dyadic co-occurrences used to learn embeddings in SGNS, we use higher-order co-occurrences to learn embeddings via HOSGNS.
As discussed in Section 3.1, the statistics of higher-order relations can be summarized in multivariate PMI tensors, which derive from proper co-occurrence probabilities among elements. Once such PMI tensors are constructed, we can again factorize them via HOSGNS. To show the versatility of this approach, we choose PMI tensors derived from two different types of higher-order relations:
1. A 3rd-order tensor P(stat)(H) ∈ R|V|×|V|×|T | which gather relative frequencies of nodes occurrences in temporal edges:
(P(stat))ijk = ω(i, j, k)
vol(H) (3.5)
where vol(H) = ∑i,j,k ω(i, j, k) is the total weight of interactions occurring in H. These probabilities are associated to the snapshot sequence representation Astat(H) = [A(1), . . . ,A(|T |)] and contain information about the topological structure ofH.
2. A 4th-order tensor P(dyn)(H) ∈ R|V|×|V|×|T |×|T |, which gather occurrence probabilities of time-stamped nodes over random walks of the supra-adjacency graph GH (as used in DYANE). Using the numerator of Eq. (2.1) tensor entries are given by:
(P(dyn))ijkl = 1
2T T∑ r=1 [ d(ik) vol(GH) (Pr)(ik)(jl) + d(jl) vol(GH) (Pr)(jl)(ik) ] (3.6)
where (ik) and (jl) are lexicographic indices of the supra-adjacency matrix A(GH) corresponding to nodes i(k) and node j(l). These probabilities encode causal dependencies among temporal nodes and are correlated with dynamical properties of spreading processes.
We also combined the two representations in a single tensor that is the average of P(stat) and P(dyn)
(P(stat|dyn))ijkl = 1
2
[ (P(stat))ijkδkl + (P(dyn))ijkl ] (3.7)
where δkl = 1[k = l] is the Kronecker delta.
Figure 2 summarizes the differences between graph embedding via classical SGNS and time-varying graph embedding via HOSGNS. Here, indices (i, j, k, l) correspond to (node, context, time, contexttime) in a 4th-order tensor representation ofH. The above tensors gather empirical probabilities PD(i, j, k . . . ) corresponding to positive examples of observable higher-order relations. The probabilities of negative examples PN (i, j, k . . . ) can be
Under review as a conference paper at ICLR 2021
<latexit sha1_base64="Q4CeBecyp8O/0jf6b0ct3gMnVwc=">AAAB+nicbVDLSgMxFL1TX7W+prp0EyxCBSkzUlAXQlFElxXsA9qhZNJMG5p5kGSUMvZT3LhQxK1f4s6/MdPOQlsPBA7n3Ms9OW7EmVSW9W3klpZXVtfy64WNza3tHbO425RhLAhtkJCHou1iSTkLaEMxxWk7EhT7Lqctd3SV+q0HKiQLg3s1jqjj40HAPEaw0lLPLHZ9rIYE8+Tmotw8vj6a9MySVbGmQIvEzkgJMtR75le3H5LYp4EiHEvZsa1IOQkWihFOJ4VuLGmEyQgPaEfTAPtUOsk0+gQdaqWPvFDoFyg0VX9vJNiXcuy7ejINKue9VPzP68TKO3MSFkSxogGZHfJijlSI0h5QnwlKFB9rgolgOisiQywwUbqtgi7Bnv/yImmeVOxq5fyuWqpdZnXkYR8OoAw2nEINbqEODSDwCM/wCm/Gk/FivBsfs9Gcke3swR8Ynz/0w5Mr</latexit>
<latexit sha1_base64="LirSwMR5Yx1+7KqOq1UjwmeB8Fg=">AAAB/HicbVDLSgMxFL3js9bXaJdugkWoUMqMFNSFUBShywp9QTuUTJq2oZkHSUYYhvorblwo4tYPceffmGlnoa0HAodz7uWeHDfkTCrL+jbW1jc2t7ZzO/ndvf2DQ/PouC2DSBDaIgEPRNfFknLm05ZiitNuKCj2XE477vQu9TuPVEgW+E0Vh9Tx8NhnI0aw0tLALPQ9rCYE86R+U2qX78vN89nALFoVaw60SuyMFCFDY2B+9YcBiTzqK8KxlD3bCpWTYKEY4XSW70eShphM8Zj2NPWxR6WTzMPP0JlWhmgUCP18hebq740Ee1LGnqsn06hy2UvF/7xepEZXTsL8MFLUJ4tDo4gjFaC0CTRkghLFY00wEUxnRWSCBSZK95XXJdjLX14l7YuKXa1cP1SLtdusjhycwCmUwIZLqEEdGtACAjE8wyu8GU/Gi/FufCxG14xspwB/YHz+AA4qk8A=</latexit>
obtained as the product of marginal distributions PD(i), PD(j), PD(k) . . . Computing exactly the objective function in Eq. (3.3) (or the 4th-order analogous) is computationally expensive, but it can be approximated by a sampling strategy: picking positive tuples according to the data distribution PD and negative ones according to independent sampling PN , HOSGNS objective can be asymptotically approximated through the optimization of the following weighted cross entropy loss:
L(bce) = − 1 B [ B∑ (ijk... )∼PD log σ ( [[wi, cj , tk . . . ]] ) +κ· B∑ (ijk... )∼PN log σ ( −[[wi, cj , tk . . . ]] )] (3.8)
where B is the number of the samples drawn in a training step and κ is the negative sampling constant. We additionally apply the warm-up steps explained in Supplementary Information to speed-up the training convergence.

heading: 4 EXPERIMENTS
text: For the experiments we use time-varying graphs collected by the SocioPatterns collaboration (http://www.sociopatterns.org) using wearable proximity sensors that sense the face-to-face proximity relations of individuals wearing them. After training the proposed models (HOSGNS applied to P(stat) , P(dyn) or P(stat|dyn)) on each dataset, embedding matrices W,C,T (and S in case of P(stat)) are mapped to embedding vectors wi, cj , tk (and sl) where i, j ∈ V and k, l ∈ T , and we use them to solve different downstream tasks: node classification and temporal event reconstruction.

heading: 4.1 EXPERIMENTAL SETUP
text: Datasets. We performed experiments with both empirical and synthetic datasets describing face-toface proximity of individuals. We used publicly available empirical contact data collected by the SocioPatterns collaboration (Cattuto et al., 2010), with a temporal resolution of 20 seconds, in a variety of contexts: in a school (“LYONSCHOOL”), a conference (“SFHH”), a hospital (“LH10”), a highschool (“THIERS13”), and in offices (“INVS15”) (Génois & Barrat, 2018). This is currently the largest collection of open datasets sensing proximity in the same range and temporal resolution used by modern contact tracing systems. In addition, we used social interactions data generated by the agent-based-model OpenABM-Covid19 (Hinch et al., 2020) to simulate an outbreak of COVID-19 in a urban setting.
We built a time-varying graph from each dataset, and for the empirical data we performed aggregation on 600 seconds time windows, neglecting those snapshots without registered interactions at that time scale. The weight of the link (i, j, k) is the number of events recorded between nodes (i, j) in a certain aggregated window k. For synthetic data we maintained the original temporal resolution and we set links weights to 1. Table 1 shows statistics for each dataset.
Table 1: Summary statistics about empirical and synthetic time-varying graph data. In order: number of single nodes |V|, number of steps |T |, number of events |E|, number of active nodes |V(T )|,
average weight of events 1|E|
∑
e∈E ω(e), nodes density |V(T )| |V||T | and links density 2|E| |V|(|V|−1)|T | .
Dataset |V| |T | |E| |V(T )| Average Weight Nodes Density Links Density
LYONSCHOOL 242 104 44820 17174 2.806 0.6824 0.0148 Empirical SFHH 403 127 17223 10815 4.079 0.2113 0.0017 graphs LH10 76 321 7435 4880 4.448 0.2000 0.0081
THIERS13 327 246 35862 32546 5.256 0.4046 0.0027 INVS15 217 691 18791 22451 4.164 0.1497 0.0012
Synthetic OPENABM-COVID19-2k-100 2000 100 1243551 198537 1.0 0.9927 0.0062 graphs OPENABM-COVID19-5k-20 5000 20 632523 99966 1.0 0.9997 0.0025
Baselines. We compare our approach with several baseline methods from the literature of timevarying graph embeddings, which learn time-stamped node representations: (1) DYANE (Sato et al., 2019), which learns temporal node embeddings with DEEPWALK, mapping a time-varying graph into a supra-adjacency representation; (2) DYNGEM (Goyal et al., 2017), a deep autoencoder architecture which dynamically reconstructs each graph snapshot initializing model weights with parameters learned in previous time frames; (3) DYNAMICTRIAD (Zhou et al., 2018), which captures structural information and temporal patterns of nodes, modeling the triadic closure process. Details about hyper-parameters used in each method can be found in the Supplementary Information.

heading: 4.2 DOWNSTREAM TASKS
text: Node Classification. The aim of this task is to classify nodes in epidemic states according to a SIR epidemic process with infection rate β and recovery rate µ. We simulated 30 realizations of the SIR process on top of each empirical graph with different combinations of parameters (β, µ). We used similar combinations of epidemic parameters and the same dynamical process to produce SIR states as described in Sato et al. (2019). Then we set a logistic regression to classify epidemic states S-I-R assigned to each active node i(k) during the unfolding of the spreading process. We combine the embedding vectors of HOSGNS using the Hadamard (element-wise) product wi ◦ tk. We compared with dynamic node embeddings learned from baselines. For fair comparison, all models produce time-stamped node representations with dimension d = 128 as input to the logistic regression.
Temporal Event Reconstruction. In this task, we aim to determine if an event (i, j, k) is in H = (V, E , T ), i.e., if there is an edge between nodes i and j at time k. We create a random timevarying graphH∗ = (V, E∗, T ) with same active nodes V(T ) and a number of |E| events that are not part of E . Embedding representations learned fromH are used as features to train a logistic regression to predict if a given event (i, j, k) is in E or in E∗. We combine the embedding vectors of HOSGNS as follows: for HOSGNS(stat), we use the Hadamard product wi ◦ cj ◦ tk; for HOSGNS(dyn) and HOSGNS(stat|dyn), we use wi◦cj ◦tk◦sk. For baseline methods, we aggregate vector embeddings to obtain link-level representations with binary operators (Average, Hadamard, Weighted-L1, WeightedL2 and Concat) as already used in previous works (Grover & Leskovec, 2016; Tsitsulin et al., 2018). For fair comparison, all models are required produce event representations with dimension d = 192.
Tasks were evaluated using train-test split. To avoid information leakage from training to test, we randomly split V and T in train and test sets (Vtr,Vts) and (Ttr, Tts), with proportion 70%− 30%. For node classification, only nodes in Vtr at times in Ttr were included in the train set, and only nodes in Vts at times in Tts were included in the test set. For temporal event reconstruction, only events with i, j ∈ Vtr and k ∈ Ttr were included in the train set, and only events with i, j ∈ Vts and k ∈ Tts were included in the test set.

heading: 4.3 RESULTS
text: In this section we first show downstream task performance results for the empirical datasets, leaving results for synthetic datasets to Supplementary Information. Synthetic datasets are used here to compare the performance of the different approaches in terms of training complexity, by measuring the number of trainable parameters and the training time with fixed number of training steps. All approaches were evaluated for both downstream tasks in terms of Macro-F1 scores in all datasets. 5 different runs of the embedding model are evaluated on 30 different train-test splits for both downstream tasks. We report the average score with standard error over all splits. In node classification, every SIR realization is assigned to a single embedding run to compute prediction scores. In event reconstruction, a different random realizationH∗ is assigned to each train-test subset. Results for the classification of nodes in epidemic states are shown in Table 2. We report here a subset of (β, µ) but other combinations are available on the Supplementary Information. DYNGEM and DYNAMICTRIAD have low scores, since they are not devised to learn from graph dynamics. HOSGNS(stat) is not able to capture the graph dynamics due to the static nature of P(stat). DYANE, HOSGNS(stat|dyn) and HOSGNS(dyn) show good performance, with these two HOSGNS variants outperforming DYANE in most of the combinations of datasets and SIR parameters.
Results for the temporal event reconstruction task are reported in Table 3. Temporal event reconstruction is not performed well by DYNGEM. DYNAMICTRIAD has better performance with Weighted-L1 and Weighted-L2 operators, while DYANE has better performance using Hadamard or WeightedL2. Since Hadamard product is explicitly used in Eq. (3.1) to optimize HOSGNS, all HOSGNS variants show best scores with this operator. HOSGNS(stat) outperforms all approaches, setting new state-of-the-art results in this task. The P(dyn) representation used as input to HOSGNS(dyn) does not focus on events but on dynamics, so the performance for event reconstruction is slightly below DYANE, while HOSGNS(stat|dyn) is comparable to DYANE. Results for HOSGNS models using other operators are available in the Supplementary Information. We observe an overall good performance of HOSGNS(stat|dyn) in both downstream tasks, being in almost all cases the second highest score, compared to the other two variants which excel in one task but have lower performance in the other.
Training Complexity. We report in Table 4 the number of trainable parameters and training time duration for each considered algorithm, when applied to an empirical graph (LYONSCHOOL) and to the synthetic ones. The proposed HOSGNS model requires a number of trainable parameters that is orders of magnitude smaller than other approaches, with a training time considerably shorter as the number of nodes increases, given a fixed number of training iterations. For this analysis, HOSGNS
sampling was implemented by picking positive and negative examples from a corpus of random walks sampled from a given graph. For HOSGNS(stat) random walks are sampled from the set of temporal snapshots {G(k)}k∈T with window size T = 1, and for HOSGNS(dyn) random walks are sampled from the supra-adjacency graph GH with window size T = 10. With these sampling strategies, positive examples are drawn from the same probability distributions as defined in Eq. (3.5) and Eq. (3.6).
We recall that HOSGNS, by learning disentangled representations of nodes and time intervals, uses a number of parameters in the order ofO(|V|+ |T |), while models that learn node-time representations (such as DYANE) need a number of parameters that is at least O(|V| · |T |). In the Supplementary Information we include plots with two dimensional projections of these embeddings, showing that the embedding matrices of HOSGNS approaches successfully capture both the structure and the dynamics of the time-varying graph.

heading: 5 CONCLUSIONS
text: In this paper, we introduce higher-order skip-gram with negative sampling (HOSGNS) for timevarying graph representation learning. We show that this method is able to disentangle the role of nodes and time, with a small fraction of the number of parameters needed by other methods. The embedding representations learned by HOSGNS outperform other methods in the literature and set new state-of-the-art results for predicting the outcome of dynamical processes and for temporal event reconstruction. We show that HOSGNS can be intuitively applied to time-varying graphs, but this methodology can be easily adapted to solve other representation learning problems that involve multi-modal data and multi-layered graph representations.

Review summaries:
In this paper, the authors propose learning node embeddings of time varying graphs. They extend the ideas from Skip Gram Negative Sampling (SGNS) to time varying graphs. They extend the relationship between SGNS and Matrix Factorization to a tensor setting. The key contribution seems to be learning a static embedding for each node and an embedding for a time step. These embeddings are combined to learn a time-aware node embedding. Experiments on multiple datasets show that the proposed method outperforms related benchmarks. 

In this paper, the authors studied the problem of time-varying graph embedding problems. The authors generalized skip-gram based graph embedding method to time-varying graphs. The authors show that the method can be used to factorize time-varying graphs as high-order tensors via negative sampling. The authors carried out experiments on several time-resolved proximity networks with comparison to several state-of-art baselines.

The paper proposes an implicit tensor factorization approach for learning time-varying node representations over dynamic networks.  The core method lifts the well-known skip gram based embedding approach from matrix to higher order tensors to support temporal dimensions. The authors claim that such tensor based treatment allows to disentangle the role of node and time. Negative sampling method ( similar to noise contrastive estimation) is extended the higher order tensor setting and incorporated in the cross entropy objective for training. In the experiments, the authors consider five variants of face-to-face proximity data that contains temporal interactions and focuses on tasks of node classification (predicting outcome of SIR epidemic process) and link prediction (in the form of event reconstruction). The proposed method has been compared against two discrete time graph representation learning model and a recently proposed tensor based method. The authors claim that the provided method shows comparable performance with requirement to train lesser number of parameters. Also, the authors provide qualitative analysis in terms of embedding visualizations and goodness of fit plots.