# Piecewise-Velocity Model for Learning Continuous-time Dynamic Node
Representations

Abdulkadir Çelikkanat, Nikolaos Nakis, Morten Mørup
Technical University of Denmark
Kongens Lyngby 2800, Denmark
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Networks have become indispensable and ubiquitous structures in many fields to
model the interactions among different entities, such as friendship in social
networks or protein interactions in biological graphs. A major challenge is to
understand the structure and dynamics of these systems. Although networks
evolve through time, most existing graph representation learning methods
target only static networks. Whereas approaches have been developed for the
modeling of dynamic networks, there is a lack of efficient continuous time
dynamic graph representation learning methods that can provide accurate
network characterization and visualization in low dimensions while explicitly
accounting for prominent network characteristics such as homophily and
transitivity. In this paper, we propose the Piecewise-Velocity Model (PiVeM)
for the representation of continuous-time dynamic networks. It learns dynamic
embeddings in which the temporal evolution of nodes is approximated by
piecewise linear interpolations based on a latent distance model with
piecewise constant node-specific velocities. The model allows for analytically
tractable expressions of the associated Poisson process likelihood with
scalable inference invariant to the number of events. We further impose a
scalable Kronecker structured Gaussian Process prior to the dynamics
accounting for community structure, temporal smoothness, and disentangled
(uncorrelated) latent embedding dimensions optimally learned to characterize
the network dynamics. We show that PiVeM can successfully represent network
structure and dynamics in ultra-low two-dimensional embedding spaces. We
further extensively evaluate the performance of the approach on various
networks of different types and sizes and find that it outperforms existing
relevant state-of-art methods in downstream tasks such as link prediction. In
summary, PiVeM enables easily interpretable dynamic network visualizations and
characterizations that can further improve our understanding of the intrinsic
dynamics of time-evolving networks.

##  1 Introduction

With technological advancements in data storage and production systems, we
have witnessed the massive growth of graph (or network) data in recent years,
with many prominent examples, including social, technological, and biological
networks from diverse disciplines [1]. They propose an exquisite way to store
and represent the interactions among data points and machine learning
techniques on graphs have thus gained considerable attention to extract
meaningful information from these complex systems and perform various
predictive tasks. In this regard, Graph Representation Learning (GRL)
techniques have become a cornerstone in the field through their exceptional
performance in many downstream tasks such as node classification and edge
prediction. Unlike the classical techniques relying on the extraction and
design of handcrafted feature vectors peculiar to given networks, GRL
approaches aim to design algorithms that can automatically learn features
optimally preserving various characteristics of networks in their induced
latent space.

Many networks evolve through time and are liable to modifications in structure
with newly arriving nodes or emerging connections, the GRL methods have
primarily addressed static networks, in other words, a snapshot of the
networks at a specific time. However, recent years have seen increasing
efforts toward modeling dynamic complex networks, see also [2] for a review.
Whereas most approaches have concentrated their attention on discrete-time
temporal networks, which have built upon a collection of time-stamped networks
(c.f. [3, 4, 5, 6, 7, 8, 9, 2, 10]) modeling of networks in continuous time
has also been studied (c.f. [11, 12, 13, 14]). These approaches have been
based on latent class [3, 11, 12, 13, 4] and latent feature modeling
approaches [5, 6, 7, 8, 9, 14, 2, 10], including advanced dynamic graph neural
network representations [15, 16].

Although these procedures have enabled the characterization of evolving
networks for downstream tasks such as link prediction and node classification,
existing dynamic latent feature models are either in discrete time or do not
explicitly account for network homophily and transitivity in terms of their
latent representations. Whereas latent class models typically provide
interpretable representations at the level of groups, latent feature models in
general rely on high-dimensional latent representations that are not easily
amenable to visualization and interpretation. A further complication of most
existing dynamic modeling approaches is their scaling typically growing in
complexity by the number of observed events and number of network dyads.

This work addresses the embedding problem of nodes in a continuous-time latent
space and seeks to model network interaction patterns using low-dimensional
representations accurately. We model the node interactions with Nonhomogeneous
Poisson Point Processes whose densities are defined based on the relative
distances among the node trajectories in the latent space. The node movements
are characterized by node-specific piecewise velocity vectors, such that each
node acquires a dynamic representation pursuing a continuous path in the
latent space throughout the timeline. The main contributions of the paper can
be summarized as follows:

  * •

We propose a novel scalable GRL method, the Piecewise-Velocity Model (PiVeM),
to flexibly learn continuous-time dynamic node representations. The temporal
evolutions of networks are represented by piecewise linear motions of the
nodes’ embeddings in the latent space.

  * •

We present a framework balancing the trade-off between the smoothness of node
trajectories in the latent space and model capacity accounting for the
temporal evolution.

  * •

We show that the PiVeM can embed nodes accurately in very low dimensional
spaces, i.e., $D=2$, such that it serves as a dynamic network visualization
tool facilitating human insights into networks’ complex, evolving structures.

  * •

The performance of the introduced approach is extensively evaluated in various
downstream tasks, such as network reconstruction and link prediction. We show
that it outperforms well-known baseline methods on a wide range of datasets.
Besides, we propose an efficient model optimization strategy enabling the
PiVeM to scale to large networks.

Source code and other materials. The datasets, the implementation, and all the
generated animations can be found at the address:
https://abdcelikkanat.github.io/projects/pivem/.

##  2 Related Work

The work on dynamic modeling of complex networks has spurred substantial
attention in recent years and covers approaches for the modeling of dynamic
structures at the level of groups (i.e., latent class models) and dynamic
representation learning approaches based on latent feature models, including
graph neural networks (GNNs). Whereas most attention has been given to
discrete-time dynamic networks, a substantial body of work has also covered
continuous-time modeling, as outlined below.

Dynamic Latent Class Models. Initial efforts for modeling continuously
evolving networks has combined latent class models defined by the stochastic
block models [17, 18] with Hawkes processes [19, 20]. In the work of [11], co-
dependent (through time) Hawkes processes were combined with the Infinite
Relational Model [21] (Hawkes IRM), yielding a non-parametric Bayesian
approach capable of expressing reciprocity between inferred groups of actors.
A drawback of such a model is the computational cost of the imposed Markov-
chain Monte-Carlo optimization, as well as, its limitation on modeling only
reciprocation effects. Scalability issues were addressed in [12] via the Block
Hawkes Model (BHM), which utilizes variational inference and simplifies the
Hawkes IRM model by associating only the inferred block structure pairs with a
univariate point process. Recently, the BHM model was extended to decoupling
interactions between different pairs of nodes belonging to the same block
pair, through the use of independent univariate Hawkes processes, defining the
Community Hawkes Independent Pairs model [13]. Whereas the above works have
been based on continuous time modeling of dynamic networks, the dynamic-IRM
(dIRM) of [3] focused on the modeling of discrete-time networks by inducing an
infinite Hidden Markov Model (IHMM) to account for transitions over time of
nodes between communities. In [4], a dynamic hierarchical block model was
proposed based on the modeling of change points admitting dynamic node
relocation within a Gibbs fragmentation tree. Despite the various advantages
of such models, networks are constrained to be regarded and analyzed at a
block level which in many cases is restrictive.

Dynamic Latent Feature Models. Prominent works around node-level
representations of continuous-time networks [22, 23] have originally
considered feature propagation within the discrete time network topology [5]
or extended the random-walk frameworks [6, 7] to the temporal case yielding
the CTDNE [24] model. CTDNE provides a single temporal-aware node embedding,
meaning that network and node evolution are unable to be visualized and
explored. A more flexible approach was designed in [15] (DyRep), where
temporal node embeddings are learned under a so-called latent mediation
process, combining an association process describing the dynamics of the
network with a communication process describing the dynamics on the network.
It uses deep recurrent architectures to parameterize the intensity function of
the point process, and thus the embedding space suffers from a lack of
explainability. HTNE [25] introduces a model utilizing a Hawkes process
relying on node embeddings. Unlike many approaches concentrating only on the
structural modifications occurring between nodes, MMDNE [26] explicitly
considers such pairwise micro, and network scale macro dynamics and uses a
temporal node representation learning algorithm relying on a temporal
attention point process. Graph neural networks (GNNs) can be extended to the
analysis of continuous networks via the Temporal Graph Network (TGN) [16]
where the classical encoder-decoder architecture is coupled with a memory
cell.

In the context of latent feature dynamic network models, Gaussian Processes
(GP) have been used to characterize the smoothness of the temporal dynamics.
This includes the discrete-time dynamic models considered in [8] in which
latent factors were endowed a GP prior based on radial basis kernels imposing
temporal smoothness within the latent representation. The approach was
extended in [9] to impose stochastic differential equations for the evolution
of latent factors. In [14], GPs were used for the modeling of continuous-time
dynamic networks based on Poisson and Hawkes processes, including exogenous as
well as endogenous features specified by a radial basis function prior.

Latent Distance Models (LDM) [27] have recently been shown to outperform
prominent GRL methods utilizing very-low dimensions in the static case [28,
29]. LDMs for temporal networks have been mostly studied in the discrete case
[2], considering mainly diffusion dynamics to make predictions, as firstly
studied in [30] and extended with popularity and activity effects [10]. While
all these models express homophily and transitivity in the dynamic case, they
fail to account for continuous dynamics.

Our work is inspired by these previous approaches for the modeling of dynamic
complex networks. Specifically, we make use of the latent distance model
formulation to account for homophily and transitivity, the Poisson Process for
the characterization of continuous-time dynamics, and a Gaussian Process prior
based on the radial-basis-function kernel to account for temporal smoothness
within the latent representation. Inspired by latent class models, we further
impose a structured low-rank representation of nodes based on soft-assigning
nodes to communities exhibiting similar temporal dynamics. Notably, we exploit
how LDMs as opposed to GNN approaches in general, can provide easily
interpretable yet accurate network representations in ultra-low dimensional
spaces ($D=2$), facilitating accurate dynamic network visualization and
interpretation.

##  3 Proposed Approach

Our main objective is to represent every node of a given network,
$\mathcal{G}=(\mathcal{V},\mathcal{E})$, into a low-dimensional metric space,
$(\mathsf{X},d_{\mathsf{X}})$, in which the pairwise node proximities will be
characterized by their distances in a continuous-time latent space (Objective
3.1). Since we address the continuous-time dynamic networks, the interactions
among nodes through time can vary, with new links appearing or disappearing at
any time. More precisely, we will presently consider undirected continuous-
time networks:

######  Definition 3.1.

A continuous-time dynamic undirected graph on a time interval
$\mathcal{I}_{T}\mathrel{\mathop{\mathchar 58\relax}}=[0,T]$ is an ordered
pair $\mathcal{G}=(\mathcal{V},\mathcal{E})$ where
$\mathcal{V}=\\{1,\ldots,N\\}$ is a set of nodes and
$\mathcal{E}\subseteq\\{\\{i,j,t\\}\in\mathcal{V}^{2}\times\mathcal{I}_{T}|1\leq
i<j\leq N\\}$ is a set of events or edges.

We will use the symbol, $N$, to denote the number of nodes in the vertex set
and $\mathcal{E}_{ij}[t_{l},t_{u}]\subseteq\mathcal{E}$ to indicate the set of
edges between nodes $i$ and $j$ occurring on the interval
$[t_{l},t_{u}]\subseteq\mathcal{I}_{T}$.

###  3.1 Nonhomogeneous Poisson Point Processes

The Poisson Point Processes (PPP)s are one of the natural choices widely used
to model the number of random events occurring in time or the locations in a
spatial space. PPPs are parameterized by a quantity known as the rate or the
intensity indicating the average density of the points in the underlying space
of the Poisson process. If the intensity depends on the time or location, the
point process is called Nonhomogeneous PPP (Defn. 3.2), and it is typically
adapted for applications in which the event points are not uniformly
distributed [31].

######  Definition 3.2.

[Nonhomogeneous PPP] A counting process $\\{M(t),t\geq 0\\}$ is called a
nonhomogeneous Poisson process with intensity function $\lambda(t)$, $t\geq 0$
if (i) $M(0)=0$, (ii) $M(t)$ has independent increments: i.e.,
$\left(M(t_{1})-M(t_{0})\right),\ldots,\left(M(t_{B})-M(t_{B-1})\right)$ are
independent random variables for each $0\leq t_{0}<\cdots<t_{B}$, and (iii)
$M(t_{u})-M(t_{l})$ is Poisson distributed with mean
$\int_{t_{l}}^{t_{u}}\lambda(t)dt$.

In this paper, we consider continuous-time dynamic undirected networks such
that the events (or links/edges) among nodes can occur at any point in time.
As we will examine in the following sections, these interactions do not
necessarily exhibit any recurring characteristics; instead, they vary over
time in many real networks. In this regard, we assume that the number of
links, $M[t_{l},t_{u}]$, between a pair of nodes $(i,j)\in\mathcal{V}^{2}$
($i<j$ since the graph is undirected) follows a nonhomogeneous Poisson point
process (NHPP) with intensity function $\lambda_{ij}(t)$ on the time interval
$[t_{l},t_{u})$, and the log-likelihood function can be written by

$\displaystyle\mathcal{L}(\Omega)\mathrel{\mathop{\mathchar 58\relax}}=\log
p(\mathcal{G}|\Omega)=\sum_{\begin{subarray}{c}i<j\\\
i,j\in\mathcal{V}\end{subarray}}\left(\sum_{e_{ij\in\mathcal{E}_{ij}}}\log\lambda_{ij}(e_{ij})-\int_{0}^{T}\lambda_{ij}(t)dt\right)$
(1)

where $\mathcal{E}_{i,j}\subseteq\mathcal{E}[0,T]$ is the set of links of node
pair $(i,j)\in\mathcal{V}^{2}$ on the timeline
$\mathcal{I}_{T}\mathrel{\mathop{\mathchar 58\relax}}=[0,T]$ for a network
$\mathcal{G}=(\mathcal{V},\mathcal{E})$, and $\Omega=\\{\lambda_{ij}\\}_{1\leq
i<j\leq N}$ indicates the set of intensity functions.

###  3.2 Problem Formulation

Without loss of generality, it can be assumed that the timeline starts from
$0$ and is bounded by $T\in\mathbb{R}^{+}$. Since the interactions among nodes
can occur at any time point on $\mathcal{I}_{T}=[0,T]$, we would like to
identify an accurate continuous-time node representation
$\\{r(i,t)\\}_{(i,t)\in\mathcal{V}\times\mathcal{I}_{T}}$ defined using a low-
dimensional latent space $\mathbb{R}^{D}$ $(D\ll N)$ where
$\mathbf{r}\mathrel{\mathop{\mathchar
58\relax}}\mathcal{V}\times\mathcal{I}_{T}\rightarrow\mathbb{R}^{D}$ is a map
indicating the embedding or representation of node $i\in\mathcal{V}$ at time
point $t\in\mathcal{I}_{T}$. We define our objective more formally as follows:

######  Objective 3.1.

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be a continuous-time dynamic
network and $\bm{\lambda}^{*}\mathrel{\mathop{\mathchar
58\relax}}\mathcal{V}^{2}\times\mathcal{I}_{T}\longrightarrow\mathbb{R}$ be an
unknown intensity function of a nonhomogeneous Poisson point process. For a
given metric space $(\mathsf{X},d_{\mathsf{X}})$, our purpose is to learn a
function or representation $\mathbf{r}\mathrel{\mathop{\mathchar
58\relax}}\mathcal{V}\times\mathcal{I}_{T}\rightarrow\mathsf{X}$ satisfying

$\displaystyle\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}\psi^{+}\\!\\!\left(d_{\mathsf{X}}\big{(}\mathbf{r}(i,t),\mathbf{r}(j,t)\big{)}\right)dt\approx\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}\bm{\lambda}^{*}(i,j,t)dt$
(2)

for a continuous function $\psi^{+}\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}\rightarrow\mathbb{R}^{+}$ for all
$(i,j)\in\mathcal{V}^{2}$ pairs, and for every interval
$[t_{l},t_{u}]\subseteq\mathcal{I}_{T}$.

In this work, we consider the Euclidean metric on a $D$-dimensional real
vector space, $\mathsf{X}\mathrel{\mathop{\mathchar 58\relax}}=\mathbb{R}^{D}$
and the embedding of node $i\in\mathcal{V}$ at time $t\in\mathcal{I}_{T}$ will
be denoted by $\mathbf{r}_{i}(t)\in\mathbb{R}^{D}$.

###  3.3 PiVeM: Piecewise-Velocity Model For Learning Continuous-time
Embeddings

We learn continuous-time node representations by employing the canonical
exponential link-function defining the intensity function as

$\displaystyle\lambda_{ij}(t)\mathrel{\mathop{\mathchar
58\relax}}=\exp{\left(\beta_{i}+\beta_{j}-||\mathbf{r}_{i}(t)-\mathbf{r}_{j}(t)||^{2}\right)}$
(3)

where $\mathbf{r}_{i}(t)\in\mathbb{R}^{D}$ and $\beta_{i}\in\mathbb{R}$ denote
the embedding vector at time $t$ and the bias term of node $i\in\mathcal{V}$,
respectively. Importantly, for a pair of nodes, we would like to have
embeddings close enough to each other when they have high interactions during
a particular time interval and far away from each other if they have less or
no links. For given bias terms, it can be seen by the following Lemma 3.1,
that the definition of the intensity function provides a guarantee for our
goal given in Equation (2).

######  Lemma 3.1.

For given fixed bias terms $\\{\beta_{i}\\}_{i\in\mathcal{V}}$, the node
embeddings, $\\{\mathbf{r}_{i}(t)\\}_{i\in\mathcal{V}}$, learned by the
objective given in Equation 1 within a bounded set of radius $R_{t}$ during a
time interval $[t_{l},t_{u}]$ satisfy

$\displaystyle\log\\!\left(\frac{(t_{u}\\!-t_{l})}{-\log
p_{ij}^{0}}\right)\\!+\\!(\beta_{i}\\!+\\!\beta_{j})\leq\frac{1}{(t_{u}\\!-\\!t_{l})}\\!\int_{t_{l}}^{t_{u}}\\!\\!\\!\\!||\mathbf{r}_{i}(t)\\!-\\!\mathbf{r}_{j}(t)||^{2}dt\\!\leq\\!\log\\!\left(\\!\frac{(t_{u}-t_{l})}{-\\!\log(1-p_{ij}^{>})}\\!\\!\right)\\!\\!+\\!(\beta_{i}\\!+\\!\beta_{j})\\!+\\!R_{t}$

where $p_{ij}^{0}$ and $p_{ij}^{>}$ are the probabilities of having zero and
more than zero links for the nodes $i$ and $j$.

###### Proof.

Please see the appendix for the proof. ∎

When we take a look at the lower bound provided in Lemma 3.1 for a pair of
nodes having no links ($p_{ij}^{0}\rightarrow 1$) within a particular time
interval $[t_{l},t_{u}]$, the lower bound converges to infinity, so nodes are
positioned in distant locations. Similarly, the upper bound provides us an
intuition about the position of the nodes for the case in which they exhibit a
high number of interactions ($p_{ij}^{>}\rightarrow 1$). The log term squeezes
the distance, and the nodes are forced to be positioned close to each other.
As we will see in the following parts, we restrict the node representations
within a bounded region of radius $R_{t}$ by imposing a prior function in
order to restrain the movements in the latent space. We utilize the squared
Euclidean distance in Equation (3), which is not a metric, but we presently
impose it as a distance [29, 32] for computational convenience, see Lemma A.4
in Appendix.

Notably, constraining the approximation of the unknown intensity function by a
metric space imposes the homophily property (i.e., similar nodes in the graph
are placed close to each other in embedding space). It can also be seen that
the transitivity property holds up to some extent (i.e., if node $i$ is
similar to $j$ and $j$ similar to $k$, then $i$ should also be similar to $k$)
since we can bound the squared Euclidean distance [33, 29]. Note that the bias
terms $\\{\beta_{i}\\}_{i\in\mathcal{V}}$ are responsible for the node-
specific effects such as degree heterogeneity [28, 33], and they provide
additional flexibility to the model by acting as scaling factor for the
corresponding nodes so that, for instance, a hub node might have a high number
of interactions simultaneously without getting close to the others in the
latent space.

Since our primary purpose is to learn continuous-time node representations in
a latent space, we define the representation of node $i\in\mathcal{V}$ at time
$t\in\mathcal{I}_{T}$ based on a linear model by
$\mathbf{r}_{i}(t)\mathrel{\mathop{\mathchar
58\relax}}=\mathbf{x}_{i}^{(0)}+\mathbf{v}_{i}t$. Here, $\mathbf{x}_{i}^{(0)}$
can be considered as the initial position and $\mathbf{v}_{i}$ the velocity of
the corresponding node. However, the linear model provides a minimal capacity
for tracking the nodes and modeling their representations. Therefore, we
reinterpret the given timeline $\mathcal{I}_{T}\mathrel{\mathop{\mathchar
58\relax}}=[0,T]$ by dividing it into $B$ equally-sized bins,
$[t_{b-1},t_{b})$, $(1\leq b\leq B)$ such that
$[0,T]=[0,t_{1})\cup\cdots\cup[t_{B-1},t_{B}]$ where
$t_{0}\mathrel{\mathop{\mathchar 58\relax}}=0$ and
$t_{B}\mathrel{\mathop{\mathchar 58\relax}}=T$. By applying the linear model
for each subinterval, we obtain a piecewise linear approximation of general
intensity functions strengthening the models’ capacity. As a result, we can
write the position of node $i$ at time $t\in\mathcal{I}_{T}$ as follows:

$\displaystyle\mathbf{r}_{i}(t)\mathrel{\mathop{\mathchar
58\relax}}=\mathbf{x}^{(0)}_{i}+\Delta_{B}\mathbf{v}_{i}^{(1)}+\Delta_{B}\mathbf{v}_{i}^{(2)}+\cdots+(t\
\text{mod}(\Delta_{B}))\mathbf{v}_{i}^{\left(\lfloor
t/\Delta_{B}\rfloor+1\right)}$ (4)

where $\Delta_{B}$ indicates the bin widths, $T/B$, and $\text{mod}(\cdot)$ is
the modulo operation used to compute the remaining time. Note that the piece-
wise interpretation of the timeline allows us to track better the path of the
nodes in the embedding space, and it can be seen by Theorem 3.2 that we can
obtain more accurate trails by augmenting the number of bins.

######  Theorem 3.2.

Let $\mathbf{f}(t)\mathrel{\mathop{\mathchar
58\relax}}[0,T]\rightarrow\mathbb{R}^{D}$ be a continuous embedding of a node.
For any given $\epsilon>0$, there exists a continuous, piecewise-linear node
embedding, $\mathbf{r}(t)$, satisfying
$||\mathbf{f}(t)-\mathbf{r}(t)||_{2}<\epsilon$ for all $t\in[0,T]$ where
$\mathbf{r}(t)\mathrel{\mathop{\mathchar 58\relax}}=\mathbf{r}^{(b)}(t)$ for
all $(b-1)\Delta_{B}\leq t<b\Delta_{B}$,
$\mathbf{r}(t)\mathrel{\mathop{\mathchar 58\relax}}=\mathbf{r}^{(B)}(t)$ for
$t=T$ and $\Delta_{B}=T/B$ for some $B\in\mathbb{N}^{+}$.

###### Proof.

Please see the appendix for the proof. ∎

Prior probability. In order to control the smoothness of the motion in the
latent space, we employ a Gaussian Process (GP) [34] prior over the initial
position $\mathbf{x}^{(0)}\in\mathbb{R}^{N\times D}$ and velocity vectors
$\mathbf{v}\in\mathbb{R}^{B\times N\times D}$. Hence, we suppose that
$\text{vect}(\mathbf{x^{(0)}})\oplus\text{vect}(\mathbf{v})\sim\mathcal{N}(\bm{0},\bm{\Sigma})$
where $\bm{\Sigma}\mathrel{\mathop{\mathchar
58\relax}}=\mathbf{\lambda}^{2}(\sigma_{\Sigma}^{2}\mathbf{I}+\mathbf{K})$ is
the covariance matrix with a scaling factor $\lambda\in\mathbb{R}$. We
utilize, $\sigma_{\Sigma\in\mathbb{R}}$, to denote the noise of the
covariance, and $\text{vect}(\mathbf{z})$ is the vectorization operator
stacking the columns to form a single vector. To reduce the number of
parameters of the prior and enable scalable inference, we define $\mathbf{K}$
as a Kronecker product of three matrices $\mathbf{K}\mathrel{\mathop{\mathchar
58\relax}}=\mathbf{B}\otimes\mathbf{C}\otimes\mathbf{D}$ respectively
accounting for temporal-, node-, and dimension specific covariance structures.
Specifically, we define $\mathbf{B}\mathrel{\mathop{\mathchar
58\relax}}=\big{[}c_{\mathbf{x}^{0}}\big{]}\oplus\big{[}\exp(-{(c_{b}-\tilde{c}_{\tilde{b}})^{2}}/{2\sigma_{\mathbf{B}}^{2}})\big{]}_{1\leq
b,\tilde{b}\leq B}$ is a ${(B+1)\times(B+1)}$ matrix intending to capture the
smoothness of velocities across time-bins where
$c_{b}=\frac{t_{b-1}+t_{b}}{2}$ is the center of the corresponding bin, and
the matrix is constructed by combining the radial basis function kernel (RBF)
with a scalar term $c_{\mathbf{x}^{0}}$ corresponding to the initial position
being decoupled from the structure of the velocities. The node specific
matrix, $\mathbf{C}\in\mathbb{R}^{N\times N}$, is constructed as a product of
a low-rank matrix $\mathbf{C}\mathrel{\mathop{\mathchar
58\relax}}=\mathbf{Q}\mathbf{Q}^{\top}$ where the row sums of
$\mathbf{Q}\in\mathbb{R}^{N\times k}$ equals to $1$ $(k\ll N)$, and it aims to
extract covariation patterns of the motion of the nodes. Finally, we simply
set the dimensionality matrix to the identity:
$\mathbf{D}\mathrel{\mathop{\mathchar
58\relax}}=\mathbf{I}\in\mathbb{R}^{D\times D}$ in order to have uncorrelated
dimensions.

To sum up, we can express our objective relying on the piecewise velocities
with the prior as follows:

$\displaystyle\hat{\Omega}=\operatorname*{arg\,max}_{\Omega}\sum_{\begin{subarray}{c}i<j\\\
i,j\in\mathcal{V}\end{subarray}}\left(\sum_{e_{ij\in\mathcal{E}_{ij}}}\log\lambda_{ij}(e_{ij})-\int_{0}^{T}\lambda_{ij}(t)dt\right)+\log\mathcal{N}\left(\left[\begin{array}[]{c}\mathbf{x}^{(0)}\\\
\mathbf{v}\end{array}\right];\mathbf{0},\bm{\Sigma}\right)$ (7)

where
$\Omega=\\{\bm{\beta},\mathbf{x}^{(0)},\mathbf{v},\sigma_{\Sigma},\sigma_{\mathbf{B}},c_{\mathbf{x}^{0}},\mathbf{Q}\\}$
is the set of hyper-parameters, and $\lambda_{ij}(t)$ is the intensity
function as defined in Equation (3) based on the node embeddings,
$\mathbf{r}_{i}(t)\in\mathbb{R}^{D}$.

Figure 1: Illustrative comparison of the ground-truth embeddings, the
adjacency matrices here for illustrative purposes constructed based on
aggregating the links appearing within the corresponding time intervals and
learned node representations.

We provide the general overview of the PiVeM method in Figure 1. The first row
shows how the ground truth node embeddings evolve through time, and the dashed
curves in the latent space show the paths they have followed. The middle row
represents the adjacency matrices of the network constructed by aggregating
the links occurring within the corresponding time intervals
$[t_{init},t_{last}]$ for illustrative purposes (notably, the model operates
in continuous time and accounts for the temporal position of each edge). Each
entry of the adjacency matrices is shaded with respect to the number of links
in the intervals, so darker regions represent a higher number of links.
Finally, the last row illustrates the learned representations and their motion
histories in the latent space.

###  3.4 Optimization

Our objective given in Equation (7) is not a convex function, so the learning
strategy that we follow is of great significance in order to escape from the
local minima and for the quality of the representations. We start by randomly
initializing the model’s hyper-parameters from $[-1,1]$ except for the
velocity tensor, which is set to $0$ at the beginning. We adapt the sequential
learning strategy in learning these parameters. In other words, we first
optimize the initial position and bias terms together,
$\\{\mathbf{x}^{(0)},\bm{\beta}\\}$, for a given number of epochs; then, we
include the velocity tensor, $\\{\mathbf{v}\\}$, in the optimization process
and repeat the training for the same number of epochs. Finally, we add the
prior parameters and learn all model hyper-parameters together. We have
employed Adam optimizer [35] with learning rate $0.1$.

Computational issues and complexity. Note that we need to evaluate the log-
intensity term in Equation (7) for each $(i,j)\in\mathcal{V}^{2}$ and event
time $e_{ij}\in\mathcal{E}_{ij}$. Therefore, the computational cost required
for the whole network is bounded by
$\mathcal{O}\left(|\mathcal{V}|^{2}|\mathcal{E}|\right)$. However, we can
alleviate the computational cost by pre-computing certain coefficients at the
beginning of the optimization process so that the complexity can be reduced to
$\mathcal{O}\left(|\mathcal{V}|^{2}B\right)$. We also have an explicit formula
for the computation of the integral term since we utilize the squared
Euclidean distance so that it can be computed in at most
$\mathcal{O}(|\mathcal{V}|^{2})$ operations. Instead of optimizing the whole
network at once, we apply a batching strategy over the set of nodes in order
to reduce the memory requirements. As a result, we sample $\mathcal{S}$ nodes
for each epoch. Hence, the overall complexity for the log-likelihood function
is $\mathcal{O}\left(\mathcal{S}^{2}B\mathcal{I}\right)$ where $\mathcal{I}$
is the number of epochs and $\mathcal{S}\ll|\mathcal{V}|$. Similarly, the
prior can be computed in at most $\mathcal{O}(B^{3}D^{3}K^{2}\mathcal{S})$
operations by using various algebraic properties such as Woodbury matrix
identity and Matrix Determinant lemma [36]. To sum up, the complexity of the
proposed approach is
$\mathcal{O}(B\mathcal{S}^{2}\mathcal{I}+B^{3}D^{3}K^{2}\mathcal{S}\mathcal{I})$
(Please see the appendix for the derivations and other details).

##  4 Experiments

In this section, we extensively evaluate the performance of the proposed
Piecewise-Velocity Model with respect to the well-known baselines in
challenging tasks over various datasets of sizes and types. We consider all
networks as undirected, and the event times of links are scaled to the
interval $[0,1]$ for the consistency of experiments. We use the finest
granularity level of the given input timestamps, such as seconds and
milliseconds. We provide a brief summary of the networks below, but more
details and various statistics are reported in Table 4 in the appendix. For
all the methods, we learn node embeddings in two-dimensional space $(D=2)$
since one of the objectives of this work is to produce dynamic node embeddings
facilitating human insights into a complex network.

Experimental Setup. We first split the networks into two sets, such that the
events occurring in the last $10\%$ of the timeline are taken out for the
prediction. Then, we randomly choose $10\%$ of the node pairs among all
possible dyads in the network for the graph completion task, and we ensure
that each node in the residual network contains at least one event keeping the
number of nodes fixed. If a pair of nodes only contains events in the
prediction set and if these nodes do not have any other links during the
training time, they are removed from the networks.

For conducting the experiments, we generate the labeled dataset of links as
follows: For the positive samples, we construct small intervals of length
$2\times 10^{-3}$ for each event time (i.e., $[e-10^{-3},e+10^{3}]$ where $e$
is an event time). We randomly sample an equal number of time points and
corresponding node pairs to form negative instances. If a sampled event time
is not located inside the interval of a positive sample, we follow the same
strategy to build an interval for it, and it is considered a negative
instance. Otherwise, we sample another time point and a dyad. Note that some
networks might contain a very high number of links, which leads to
computational problems for these networks. Therefore, we subsample $10^{4}$
positive and negative instances if they contain more than this.

Synthetic networks. We generate two artificial networks in order to evaluate
the behavior of the models in controlled experimental settings. (i)
Synthetic($\pi$) is sampled from the prior distribution stated in Subsection
3.2. The hyper-parameters, $\bm{\beta}$, $K$ and $B$ are set to $\bm{0}$, $20$
and $100$, respectively. (ii) Synthetic($\mu$) is constructed based on the
temporal block structures. The timeline is divided into $10$ sub-intervals,
and the nodes are randomly split into $20$ groups for each interval. The links
within each group are sampled from the Poisson distribution with the constant
intensity of $5$.

Real networks. The (iii) Hypertext network [37] was built on the radio badge
records showing the interactions of the conference attendees for $2.5$ days,
and each event time indicates $20$ seconds of active contact. Similarly, (iv)
the Contacts network [38] was generated concerning the interactions of the
individuals in an office environment. (v) Forum [39] is comprised of the
activity data of university students on an online social forum system. (vi)
College [40] indicates the private messages among the students on an online
social platform. Finally, (vii) Email [41] was constructed based on the
exchanged e-mail information among the members of European research
institutions.

Baselines. We compare the performance of our method with five baselines. We
include (i) Ldm with Poisson rate with node-specific biases [42, 33] since it
is a static method having the closest formulation to ours. We randomly
initialize the embeddings and bias terms, and train the model with the Adam
optimizer [35] for $500$ epochs and a learning rate of $0.1$. A very well-
known GRL method, (ii) Node2Vec [7] relies on the explicit generation of the
random walks by starting from each node in the network, then it learns node
embeddings by inspiring from the SkipGram [43] algorithm. It optimizes the
softmax function for the nodes lying within a fixed window region with respect
to a chosen center node over the produced node sequences. In our experiments,
we tune the model’s parameters ($p$, $q$) from $\\{0.25,0.5,1,2,4\\}$. Since
it has the ability to run over the weighted networks, we also constructed a
weighted graph based on the number of links through time and reported the best
score of both versions of the networks. (iii) CTDNE [24] is a dynamic node
embedding approach performing temporal random walks over the network. (iv)
HTNE [25] learns embeddings based on the Hawkes process modeling the
neighborhood formation sequence induced from the network structure. (v) MMDNE
[26] introduces a temporal attention point process to model the newly
established links and proposes a general dynamics equation relying on latent
node representations to capture the network scale evolutions.

The continuous-time baseline methods are unable to produce instantaneous node
representations and they produce embeddings only for a given time. Therefore,
we have utilized the last time of the training set to obtain the
representations. We have chosen the recommended values for the common hyper-
parameters of Node2Vec and CTDNE, so the number of walks, walk length, and
window size parameters have been set to $10$, $80$, and $10$, respectively. We
used the implementation provided by the StellarGraph Python package to produce
the embeddings for CTDNE. Similarly, we have adapted the suggested
hyperparameter settings for MMDNE and CTDNE with $100$ epochs.

For our method, we set the parameter $K=25$, and bins count $B=100$ to have
enough capacity to track node interactions. For the regularization term
($\lambda$) of the prior, we first mask $20\%$ of the dyads in the
optimization of Equation (7). Furthermore, we train the model by starting with
$\lambda=10^{6}$, and then we reduce it to one-tenth after each $100$ epoch.
The same procedure is repeated until $\lambda=10^{-6}$, and we choose the
$\lambda$ value minimizing the log-likelihood of the masked pairs. The final
embeddings are then obtained by performing this annealing strategy without any
mask until this $\lambda$ value. We repeat this procedure $5$ times, and we
consider the best-performing method in learning the embeddings. The
Coefficient of Variation (CV) of the experiments is always less than $0.5$,
and Figure 2(a) shows an illustrative example for tuning $\lambda$ over the
Synthetic($\pi$) dataset with $5$ random runs.

For the performance comparison of the methods, we provide the Area Under Curve
(AUC) scores for the Receiver Operating Characteristic (ROC) and Precision-
Recall (PR) curves [44]. We compute the intensity of a given instance for LDM
and PiVeM for the similarity measure of the node pair. Since Node2Vec and
CTDNE rely on the SkipGram architecture [43], we use cosine similarity for
them.

Table 1: The performance evaluation for the network reconstruction experiment
over various datasets.

| Synthetic($\pi$) | Synthetic($\mu$) | College | Contacts | Email | Forum | Hypertext
---|---|---|---|---|---|---|---
| ROC | PR | ROC | PR | ROC | PR | ROC | PR | ROC | PR | ROC | PR | ROC | PR
LDM | $.563$ | $.539$ |  $.669$ |  $.642$ | $.951$ |  $.944$ |  $.860$ |  $.835$ |  $.954$ |  $.948$ | $.909$ |  $.897$ |  $.818$ |  $.797$
Node2Vec | $.519$ | $.507$ | $.503$ | $.509$ | $.711$ | $.655$ | $.812$ | $.756$ | $.853$ | $.828$ | $.677$ | $.619$ | $.696$ | $.648$
CTDNE | $.613$ | $.580$ | $.539$ | $.544$ | $.661$ | $.622$ | $.787$ | $.760$ | $.854$ | $.840$ | $.657$ | $.622$ | $.725$ | $.725$
HTNE |  $.614$ |  $.591$ | $.599$ | $.571$ | $.721$ | $.683$ | $.846$ | $.823$ | $.871$ | $.867$ | $.723$ | $.691$ | $.775$ | $.787$
MMDNE | $.582$ | $.565$ | $.600$ | $.576$ | $.725$ | $.692$ | $.844$ | $.825$ | $.867$ | $.863$ | $.737$ | $.712$ | $.778$ | $.787$
PiVeM | $.762$ | $.713$ | $.905$ | $.869$ |  $.948$ | $.948$ | $.938$ | $.938$ | $.978$ | $.977$ |  $.907$ | $.902$ | $.830$ | $.823$

Table 2: The performance evaluation for the network completion experiment over
various datasets.

| Synthetic($\pi$) | Synthetic($\mu$) | College | Contacts | Email | Forum | Hypertext
---|---|---|---|---|---|---|---
| ROC | PR | ROC | PR | ROC | PR | ROC | PR | ROC | PR | ROC | PR | ROC | PR
LDM | $.535$ | $.529$ | $.646$ | $.631$ |  $.931$ |  $.926$ |  $.836$ | $.799$ |  $.948$ |  $.942$ |  $.863$ |  $.858$ |  $.761$ | $.738$
Node2Vec | $.519$ | $.511$ |  $.747$ |  $.677$ | $.685$ | $.637$ | $.787$ | $.744$ | $.818$ | $.777$ | $.635$ | $.592$ | $.596$ | $.588$
CTDNE |  $.608$ | $.573$ | $.531$ | $.539$ | $.601$ | $.556$ | $.752$ | $.703$ | $.831$ | $.812$ | $.568$ | $.539$ | $.554$ | $.537$
HTNE | $.605$ |  $.583$ | $.573$ | $.557$ | $.673$ | $.651$ | $.792$ | $.759$ | $.853$ | $.834$ | $.596$ | $.581$ | $.602$ | $.633$
MMDNE | $.587$ | $.570$ | $.592$ | $.571$ | $.677$ | $.662$ | $.819$ |  $.811$ | $.844$ | $.829$ | $.596$ | $.570$ | $.587$ | $.614$
PiVeM | $.750$ | $.696$ | $.874$ | $.851$ | $.935$ | $.934$ | $.873$ | $.864$ | $.951$ | $.953$ | $.879$ | $.875$ | $.770$ |  $.712$

Table 3: The performance evaluation for the link prediction experiment over
various datasets.

| Synthetic($\pi$) | Synthetic($\mu$) | College | Contacts | Email | Forum | Hypertext
---|---|---|---|---|---|---|---
| ROC | PR | ROC | PR | ROC | PR | ROC | PR | ROC | PR | ROC | PR | ROC | PR
LDM | $.562$ | $.539$ |  $.498$ | $.642$ | $.951$ | $.944$ | $.860$ | $.835$ |  $.954$ |  $.948$ | $.909$ | $.897$ | $.819$ | $.797$
Node2Vec | $.518$ | $.506$ |  $.498$ | $.502$ | $.705$ | $.676$ | $.783$ | $.716$ | $.825$ | $.807$ | $.635$ | $.605$ | $.748$ | $.739$
CTDNE |  $.680$ |  $.629$ | $.481$ | $.487$ | $.691$ | $.711$ | $.842$ | $.815$ | $.824$ | $.815$ | $.664$ | $.642$ | $.699$ | $.734$
HTNE | $.573$ | $.569$ | $.491$ | $.493$ | $.715$ | $.684$ | $.864$ | $.824$ | $.838$ | $.837$ | $.764$ | $.747$ | $.785$ | $.820$
MMDNE | $.591$ | $.575$ | $.506$ |  $.515$ | $.717$ | $.703$ |  $.874$ |  $.847$ | $.827$ | $.832$ | $.762$ | $.746$ |  $.795$ |  $.813$
PiVeM | $.716$ | $.689$ | $.474$ | $.485$ |  $.891$ |  $.887$ | $.876$ | $.884$ | $.964$ | $.964$ |  $.894$ |  $.895$ | $.756$ | $.767$

Network Reconstruction. Our goal is to see how accurately a model can capture
the interaction patterns among nodes and generate embeddings exhibiting their
temporal relationships in a latent space. In this regard, we train the models
on the residual network and generate sample sets as described previously. The
performance of the models is reported in Table 1. Comparing the performance of
PiVeM against the baselines, we observe favorable results across all networks,
highlighting the importance and ability of PiVeM to account for and detect
structure in a continuous time manner.

Network Completion. The network completion experiment is a relatively more
challenging task than the reconstruction. Since we hide $10\%$ of the network,
the dyads containing events are also viewed as non-link pairs, and the
temporal models should place these nodes in distant locations of the embedding
space. However, it might be possible to predict these events accurately if the
network links have temporal triangle patterns through certain time intervals.
In Table 2, we report the AUC-ROC and PR-AUC scores for the network completion
experiment. Once more, PiVeM outperforms the baselines (in most cases
significantly). We again discovered evidence supporting the necessity for
modeling and tracking temporal networks with time-evolving embedding
representations.

Future Prediction. Finally, we examine the performance of the models in the
future prediction task. Here, the models are asked to forecast the $10\%$
future of the timeline. For PiVeM, the similarity between nodes is obtained by
calculating the intensity function for the timeline of the training set (i.e.,
from $0$ to $0.9$), and we keep our previously described strategies for the
baselines since they generate the embeddings only for the last training time.
Table 3 presents the performances of the models. It is noteworthy that while
PiVeM outperforms the baselines significantly on the Synthetic($\pi$) network,
it does not show promising results on Synthetic($\mu$). Since the first
network is compatible with our model, it successfully learns the dominant link
pattern of the network. However, the second network conflicts with our model:
it forms a completely different structure for every $0.1$ second. For the real
datasets, we observe mostly on-par results, especially with LDM. Some real
networks contain link patterns that become "static" with respect to the future
prediction task.

(a) Annealing strategy

(b) Influence of dimension

(c) Influence of bin count

Figure 2: Influence of the model hyperparameters over the Synthetic($\pi$)
dataset.

We have previously described how we set the prior coefficient, $\lambda$, and
now we will examine the influence of the other hyperparameters over the
Synthetic($\pi$) dataset for network reconstruction.

Influence of dimension size ($D$). We report the AUC-ROC and AUC-PR scores in
Figure 2(b). When we increase the dimension size, we observe a constant
increase in performance. It is not a surprising result because we also
increase the model’s capacity depending on the dimension. However, the two-
dimensional space still provides comparable performances in the experiments,
facilitating human insights into networks’ complex, evolving structures.

Influence of bin count ($B$). Figure 2(c) demonstrates the effect of the
number of bins for the network reconstruction task. We generated the
Synthetic($\pi$) network from for $100$ bins, so it can be seen that the
performance stabilizes around $2^{6}$, which points out that PiVeM reaches
enough capability to model the interactions among nodes.

Latent Embedding Animation. It is of great significance from a human
perspective to embed nodes into a latent space and visualize them in order to
comprehend the intricate interactions among the entities and extract
noteworthy information from the data. Although many GRL methods show high
performance in the downstream tasks, in general, they require high dimensional
spaces, so a postprocessing step later has to be applied in order to visualize
the node representations in a small dimensional space. However, such processes
cause distortions in the embeddings, which can lead a practitioner to end up
with inaccurate arguments about the data.

As we have seen in the experimental evaluations, our proposed approach
successfully learns embeddings in the two-dimensional space, and it also
produces continuous-time representations. Therefore, it offers the ability to
animate how the network evolves through time and can play a crucial role in
grasping the underlying characteristics of the networks. As an illustrative
example, Figure 3 compares the ground truth representations of
Synthetic($\pi$) with the learned ones. The synthetic network consists of
small communities of $5$ nodes, and each color indicates these groups.
Although the problem does not have unique solutions, it can be seen that our
model successfully seizes the clustering patterns in the network. We refer the
reader to supplementary materials for the full animation.

Ground Truth

(a) $t=0.1$

(b) $t=0.4$

(c) $t=0.7$

(d) $t=0.9$

Learned Model

(e) $t=0.1$

(f) $t=0.4$

(g) $t=0.7$

(h) $t=0.9$

Figure 3: Comparisons of the ground truth and learned representations in two-
dimensional space.

##  5 Conclusion and Limitations

In this paper, we have proposed a novel continuous-time dynamic network
embedding approach, namely, Piecewise Velocity Model (PiVeM). Its performance
has been examined in various experiments, such as network reconstruction and
completion tasks over various networks with respect to the very well-known
baselines. We demonstrated that it could accurately embed the nodes into a
two-dimensional space. Therefore, it can be directly utilized to animate the
learned node embeddings, and it can be beneficial in extracting the networks’
underlying characteristics, foreseeing how they will evolve through time. We
theoretically showed that the model could scale up to large networks.

Although our model successfully learns continuous-time representations, it is
unable to capture temporal patterns in the network in terms of the GP
structure. Therefore, we are planning to employ different kernels instead of
RBF, such as periodic kernels in the prior. The optimization strategies of the
proposed method might be improved to escape from local minima. As a possible
future direction, the algorithm can also be adapted for other graph types,
such as directed and multi-layer networks.

## Acknowledgements

We would like to thank the reviewers for the constructive feedback and their
insightful comments. We would also like to thank Sune Lehmann, Louis
Boucherie, Lasse Mohr Mikkelsen, Simon Tommerup, August Semrau Andersen and
William Diedrichsen Marstrand for the valuable and fruitful discussions. We
gratefully acknowledge the Independent Research Fund Denmark for supporting
this work [grant number: 0136-00315B].

## References

  * Newman [2003] M. E. J. Newman.  The structure and function of complex networks.  _SIAM Review_ , 45(2):167–256, 2003.
  * Kim et al. [2018] Bomin Kim, Kevin H Lee, Lingzhou Xue, and Xiaoyue Niu.  A review of dynamic network models with latent variables.  _Statistics surveys_ , 12:105, 2018.
  * Ishiguro et al. [2010] Katsuhiko Ishiguro, Tomoharu Iwata, Naonori Ueda, and Joshua Tenenbaum.  Dynamic infinite relational model for time-varying relational data analysis.  _NeurIPS_ , 23, 2010.
  * Herlau et al. [2013] Tue Herlau, Morten Mørup, and Mikkel Schmidt.  Modeling temporal evolution and multiscale structure in networks.  In _ICML_ , pages 960–968, 2013.
  * Heaukulani and Ghahramani [2013] Creighton Heaukulani and Zoubin Ghahramani.  Dynamic probabilistic models for latent feature propagation in social networks.  In _ICML_ , pages 275–283, 2013.
  * Perozzi et al. [2014] Bryan Perozzi, Rami Al-Rfou, and Steven Skiena.  Deepwalk: Online learning of social representations.  In _KDD_ , page 701–710, 2014.
  * Grover and Leskovec [2016] Aditya Grover and Jure Leskovec.  Node2Vec: Scalable feature learning for networks.  In _KDD_ , pages 855–864, 2016.
  * Durante and Dunson [2014] Daniele Durante and David Dunson.  Bayesian Logistic Gaussian Process Models for Dynamic Networks.  In _AISTATS_ , volume 33, pages 194–201, 2014.
  * Durante and Dunson [2016] Daniele Durante and David B Dunson.  Locally adaptive dynamic networks.  _The Annals of Applied Statistics_ , 10(4):2203–2232, 2016.
  * Sewell and Chen [2015] Daniel K. Sewell and Yuguo Chen.  Latent space models for dynamic networks.  _JASA_ , 110(512):1646–1657, 2015.
  * Blundell et al. [2012] Charles Blundell, Jeff Beck, and Katherine A Heller.  Modelling reciprocating relationships with hawkes processes.  In _NeurIPS_ , volume 25, 2012.
  * Arastuie et al. [2020] Makan Arastuie, Subhadeep Paul, and Kevin Xu.  CHIP: A hawkes process model for continuous-time networks with scalable and consistent estimation.  In _NeurIPS_ , volume 33, pages 16983–16996, 2020.
  * Delattre et al. [2016] Sylvain Delattre, Nicolas Fournier, and Marc Hoffmann.  Hawkes processes on large networks.  _The Annals of Applied Probability_ , 26(1):216 – 261, 2016.
  * Fan et al. [2021] Xuhui Fan, Bin Li, Feng Zhou, and Scott SIsson.  Continuous-time edge modelling using non-parametric point processes.  _NeurIPS_ , 34:2319–2330, 2021.
  * Trivedi et al. [2019] Rakshit Trivedi, Mehrdad Farajtabar, Prasenjeet Biswal, and Hongyuan Zha.  Dyrep: Learning representations over dynamic graphs.  In _ICLR_ , 2019.
  * Rossi et al. [2020] Emanuele Rossi, Ben Chamberlain, Fabrizio Frasca, Davide Eynard, Federico Monti, and Michael M. Bronstein.  Temporal graph networks for deep learning on dynamic graphs.  _ICML 2020 Workshop_ , 2020.
  * Holland et al. [1983] Paul W Holland, Kathryn Blackmond Laskey, and Samuel Leinhardt.  Stochastic blockmodels: First steps.  _Social networks_ , 5(2):109–137, 1983.
  * Nowicki and Snijders [2001] Krzysztof Nowicki and Tom A. B Snijders.  Estimation and prediction for stochastic blockstructures.  _JASA_ , 96(455):1077–1087, 2001.
  * Hawkes [1971a] Alan G. Hawkes.  Spectra of some self-exciting and mutually exciting point processes.  _Biometrika_ , 58(1):83–90, 1971a.
  * Hawkes [1971b] Alan G. Hawkes.  Point spectra of some mutually exciting point processes.  _J. R. Stat. Soc_ , 33(3):438–443, 1971b.
  * Kemp et al. [2006] Charles Kemp, Joshua B. Tenenbaum, Thomas L. Griffiths, Takeshi Yamada, and Naonori Ueda.  Learning systems of concepts with an infinite relational model.  In _AAAI_ , page 381–388, 2006.
  * Xue et al. [2022] Guotong Xue, Ming Zhong, Jianxin Li, Jia Chen, Chengshuai Zhai, and Ruochen Kong.  Dynamic network embedding survey.  _Neurocomputing_ , 472:212–223, 2022.
  * Kazemi et al. [2020] Seyed Mehran Kazemi, Rishab Goel, Kshitij Jain, Ivan Kobyzev, Akshay Sethi, Peter Forsyth, and Pascal Poupart.  Representation learning for dynamic graphs: A survey.  _JMLR_ , 21(70):1–73, 2020.
  * Nguyen et al. [2018] Giang Hoang Nguyen, John Boaz Lee, Ryan A. Rossi, Nesreen K. Ahmed, Eunyee Koh, and Sungchul Kim.  Continuous-time dynamic network embeddings.  In _TheWebConf_ , page 969–976, 2018.
  * Zuo et al. [2018] Yuan Zuo, Guannan Liu, Hao Lin, Jia Guo, Xiaoqian Hu, and Junjie Wu.  Embedding temporal network via neighborhood formation.  In _KDD_ , page 2857–2866, 2018.
  * Lu et al. [2019] Yuanfu Lu, Xiao Wang, Chuan Shi, Philip S. Yu, and Yanfang Ye.  Temporal network embedding with micro- and macro-dynamics.  In _CIKM_ , page 469–478, 2019.
  * Hoff et al. [2002] Peter D Hoff, Adrian E Raftery, and Mark S Handcock.  Latent space approaches to social network analysis.  _JASA_ , 97(460):1090–1098, 2002.
  * Nakis et al. [2022a] Nikolaos Nakis, Abdulkadir Çelikkanat, Sune Lehmann Jørgensen, and Morten Mørup.  A hierarchical block distance model for ultra low-dimensional graph representations, 2022a.
  * Nakis et al. [2022b] Nikolaos Nakis, Abdulkadir Çelikkanat, and Morten Mørup.  HM-LDM: A hybrid-membership latent distance model, 2022b.
  * Sarkar and Moore [2005] Purnamrita Sarkar and Andrew Moore.  Dynamic social network analysis using latent space models.  In Y. Weiss, B. Schölkopf, and J. Platt, editors, _NeurIPS_ , volume 18, 2005.
  * Streit [2010] Roy L. Streit.  _The Poisson Point Process_ , pages 11–55.  Springer US, Boston, MA, 2010.
  * Nickel and Kiela [2017] Maximillian Nickel and Douwe Kiela.  Poincaré embeddings for learning hierarchical representations.  In _Neur_ , volume 30, 2017.
  * Krivitsky et al. [2009] Pavel N. Krivitsky, Mark S. Handcock, Adrian E. Raftery, and Peter D. Hoff.  Representing degree distributions, clustering, and homophily in social networks with latent cluster random effects models.  _Social Networks_ , 31(3):204 – 213, 2009.
  * Rasmussen and Williams [2005] Carl Edward Rasmussen and Christopher K. I. Williams.  _Gaussian Processes for Machine Learning_.  The MIT Press, 11 2005.
  * Kingma and Ba [2017] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2017.
  * Aggarwal [2020] C.C. Aggarwal.  _Linear Algebra and Optimization for Machine Learning: A Textbook_.  Springer International Publishing, 2020.
  * Isella et al. [2011] Lorenzo Isella, Juliette Stehlé, Alain Barrat, Ciro Cattuto, Jean-François Pinton, and Wouter Van den Broeck.  What’s in a Crowd? Analysis of Face-to-Face Behavioral Networks.  _Journal of Theoretical Biology_ , 271(1):166–180, 2011.
  * Génois and Barrat [2018] Mathieu Génois and Alain Barrat.  Can co-location be used as a proxy for face-to-face contacts?  _EPJ Data Science_ , 7(1):11, May 2018.
  * Opsahl [2010] Tore Opsahl.  Triadic closure in two-mode networks: Redefining the global and local clustering coefficients.  _Social Networks_ , 35, 06 2010.
  * Panzarasa et al. [2009] Pietro Panzarasa, Tore Opsahl, and Kathleen M. Carley.  Patterns and dynamics of users’ behavior and interaction: Network analysis of an online community.  _J. Assoc. Inf. Sci. Technol._ , 60(5):911–932, 2009.
  * Paranjape et al. [2017] Ashwin Paranjape, Austin R. Benson, and Jure Leskovec.  Motifs in temporal networks.  In _WSDM_ , page 601–610, 2017.
  * Hoff [2005] Peter D Hoff.  Bilinear mixed-effects models for dyadic data.  _JASA_ , 100(469):286–295, 2005.
  * Mikolov et al. [2013] Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg Corrado, and Jeffrey Dean.  Distributed representations of words and phrases and their compositionality.  In _NeurIPS_ , pages 3111–3119, 2013.
  * Pedregosa et al. [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  _JMLR_ , 12:2825–2830, 2011.

##  Appendix A Appendix

In the main paper, we could not clarify every aspect of the model due to the
restrictions on the number of pages. Hence, we will provide more detailed
explanations here about the experiments, the computational problems, the
proofs of the theoretical arguments, and possible extensions of the model
toward the bipartite, weighted, and directed networks.

###  A.1 Experiments

We consider all networks used in the experiments as undirected, and the event
times of links are scaled to the interval $[0,1]$ for the consistency of
experiments. We use the finest resolution level of the given input timestamps,
such as seconds and milliseconds. We provide a brief summary of the networks
below, and various statistics are reported in Table 4. The visualization of
the event distributions of the networks through time is depicted in Figure 4.

(a) Synthetic($\pi$)

(b) Synthetic($\mu$)

(c) College

(d) Contacts

(e) Email

(f) Forum

(g) Hypertext

Figure 4: Distribution of the links through time.

Synthetic datasets. We generate two artificial networks in order to evaluate
the behavior of the models in controlled experimental settings. (i)
Synthetic($\pi$) is sampled from the prior distribution stated in Subsection
3.2. The hyper-parameters, $\bm{\beta}$, $K$ and $B$ are set to $\bm{0}$, $20$
and $100$, respectively. (ii) Synthetic($\mu$) is constructed based on the
temporal block structures. The timeline is divided into $10$ intervals, and
the node set is split into $20$ groups. The links within each group are
sampled from the Poisson distribution with the constant intensity of $5$.

Real Networks. The (iii) Hypertext network [37] was built on the radio badge
records showing the interactions of the conference attendees for $2.5$ days,
and each event time indicates $20$ seconds of active contact. Similarly, (iv)
the Contacts network [38] was generated concerning the interactions of the
individuals in an office environment. (v) Forum [39] is comprised of the
activity data of university students on an online social forum system. The
(vi) CollegeMsg network [40] indicates the private messages among the students
on an online social platform. Finally, (vii) Eu-Email [41] was constructed
based on the exchanged e-mail information among the members of European
research institutions.

Table 4: Statistics of networks. $|\mathcal{V}|$: Number of nodes, $M$: Number of pairs having at least one link, $|\mathcal{E}|$: Total number of links, $\left|\mathcal{E}_{ij}\right|_{max}$: Max. number of links a pair of nodes has. | $\left|\mathcal{V}\right|$ | $M$ | $\left|\mathcal{E}\right|$ | $\left|\mathcal{E}_{ij}\right|_{max}$
---|---|---|---|---
Synthetic($\mu$) | 100 | 4,889 | 180,658 | 124
Synthetic($\pi$) | 100 | 3,009 | 22,477 | 32
College | 1,899 | 13,838 | 59,835 | 184
Contacts | 217 | 4,274 | 78,249 | 1,302
Hypertext | 113 | 2,196 | 20,818 | 1,281
Email | 986 | 16,064 | 332,334 | 4,992
Forum | 899 | 7,036 | 33,686 | 171

###  A.2 Computational Problems, Model Complexity and Optimization Strategy

Log-likelihood function. Note that we need to evaluate the log-intensity term
in Equation 7 for each $(i,j)\in\mathcal{V}^{2}$ $(i<j)$ and event time
$e_{ij}\in\mathcal{E}_{ij}$. Therefore, the computational cost required for
the whole network is bounded by
$\mathcal{O}\left(|\mathcal{V}|^{2}|\mathcal{E}|\right)$. However, we can
alleviate it by computing certain coefficients at the beginning of the
optimization process. If we define $\alpha_{ij}\mathrel{\mathop{\mathchar
58\relax}}=(e_{ij}-\Delta_{B}(b^{*}-1))$, then it can be seen that the sum
over the set of all events, $\mathcal{E}_{ij}^{b^{*}}$, lying inside
$b^{*}$’th bin (i.e., the events in $[\Delta_{B}(b^{*}-1),\ \
\Delta_{B}b^{*})$ can be rewritten by:

$\displaystyle\sum_{e_{ij}\in\mathcal{E}_{ij}^{b^{*}}}\\!\\!\\!\\!\log\lambda_{ij}(e_{ij})\\!$
$\displaystyle=\\!\\!\\!\\!\\!\sum_{e_{ij}\in\mathcal{E}_{ij}}\\!\\!\\!\\!\Big{(}\beta_{i}\\!+\\!\beta_{j}-||\mathbf{r}_{i}(e_{ij})-\mathbf{r}_{j}(e_{ij})||^{2}\Big{)}$
$\displaystyle=\\!\\!\\!\\!\\!\sum_{e_{ij}\in\mathcal{E}_{ij}^{b^{*}}}\\!\\!\\!\\!\Big{(}\beta_{i}\\!+\\!\beta_{j}\Big{)}\\!+\\!\\!\\!\sum_{e_{ij}\in\mathcal{E}_{ij}}\Big{|}\Big{|}\Delta\mathbf{x}^{(0)}_{ij}+\Delta_{B}\\!\\!\sum_{b=1}^{b^{*}-1}\Delta\mathbf{v}_{ij}^{(b)}+\Delta\mathbf{v}_{ij}^{(b^{*})}(e_{ij}\\!-\\!\Delta_{B}(b^{*}\\!-\\!1))\Big{|}\Big{|}^{2}$
$\displaystyle=\\!\\!\\!\\!\\!\sum_{e_{ij}\in\mathcal{E}_{ij}^{b^{*}}}\\!\\!\\!\\!\Big{(}\beta_{i}\\!+\\!\beta_{j}\Big{)}\\!+\\!\\!\\!\sum_{e_{ij}\in\mathcal{E}_{ij}}\Bigg{(}\alpha_{ij}^{2}\Big{|}\Big{|}\Delta\mathbf{v}_{ij}^{(b^{*})}\Big{|}\Big{|}^{2}+\Big{(}\Delta\mathbf{x}^{(0)}_{ij}+\Delta_{B}\\!\\!\sum_{b=1}^{b^{*}-1}\Delta\mathbf{v}_{ij}^{(b)}\Big{)}^{2}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad+2\alpha_{ij}\Big{\langle}\Delta\mathbf{x}^{(0)}_{ij}+\Delta_{B}\sum_{b=1}^{b^{*}-1}\Delta\mathbf{v}_{ij}^{(b)},\Delta\mathbf{v}_{ij}^{(b^{*})}\Big{\rangle}\Bigg{)}$
$\displaystyle=\left|\mathcal{E}_{ij}^{b^{*}}\right|(\beta_{i}\\!+\\!\beta_{j})\\!+\alpha_{2}\Big{|}\Big{|}\Delta\mathbf{v}_{ij}^{(b^{*})}\Big{|}\Big{|}^{2}+\sum_{e_{ij}\in\mathcal{E}_{ij}}\\!\\!\\!\\!\Big{(}\Delta\mathbf{x}^{(0)}_{ij}+\Delta_{B}\\!\\!\sum_{b=1}^{b^{*}-1}\Delta\mathbf{v}_{ij}^{(b)}\Big{)}^{2}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad+2\alpha_{1}\Big{\langle}\Delta\mathbf{x}^{(0)}_{ij}+\Delta_{B}\sum_{b=1}^{b^{*}-1}\Delta\mathbf{v}_{ij}^{(b)},\Delta\mathbf{v}_{ij}^{(b^{*})}\Big{\rangle}$

where $\alpha_{1}^{(b^{*})}\mathrel{\mathop{\mathchar
58\relax}}=\sum_{e_{ij}\in\mathcal{E}_{ij}}\alpha_{ij}$ and
$\alpha_{2}^{(b^{*})}\mathrel{\mathop{\mathchar
58\relax}}=\sum_{e_{ij}\in\mathcal{E}_{ij}}\alpha_{ij}^{2}$. We can follow the
same strategy for each bin, then the computational complexity can be reduced
to $\mathcal{O}\left(|\mathcal{V}|^{2}B\right)$

Since we use the squared Euclidean distance in the integral term of our
objective, we can derive the exact formula for the computation (please see
Lemma A.4 for the details). We need to evaluate it for all node pairs, so it
requires at most $\mathcal{O}\left(|\mathcal{V}|^{2}\right)$ operations.
Hence, the complexity of the log-likelihood function is
$\mathcal{O}\left(|\mathcal{V}|^{2}B\right)$. Instead of optimizing the whole
network at once, we are applying the batching strategy over the set of nodes
in order to reduce the memory requirements, so we sample $\mathcal{S}$ nodes
for each epoch. Hence, the overall complexity of the log-likelihood is
$\mathcal{O}\left(\mathcal{S}^{2}B\mathcal{I}\right)$ where $\mathcal{I}$ is
the number of epochs.

Computation of the prior function. The covariance matrix,
$\mathbf{\Sigma}\in\mathbb{R}^{BND\times BND}$, of the prior is defined by
$\mathbf{\Sigma}\mathrel{\mathop{\mathchar
58\relax}}=\lambda^{2}\left(\sigma_{\mathbf{\Sigma}}^{2}\mathbf{I}+\mathbf{K}\right)^{-1}$
with a scaling factor $\lambda\in\mathbb{R}$ and a noise variance
$\sigma_{\mathbf{\Sigma}}^{2}\in\mathbb{R}^{+}$. The multivariate normal
distribution is parametrized with a noise term
$\sigma_{\mathbf{\Sigma}}^{2}\mathbf{I}$ and a matrix
$\mathbf{K}\in\mathbb{R}^{BND\times BND}$ having a low-rank form. In other
words, $\mathbf{K}$ is written by
$\mathbf{B}\otimes\mathbf{C}\otimes\mathbf{D}$ where $\mathbf{B}$ is block
diagonal matrix combined with parameter $c_{\mathbf{x}^{0}}$ and the RBF
kernel
$\exp\left(-{(c_{b}-c_{b^{{}^{\prime}}})^{2}}/{\sigma_{\mathbf{B}}^{2}}\right)\in\mathbb{R}^{B\times
B}$ for $c_{b}\mathrel{\mathop{\mathchar 58\relax}}=(t_{b-1}-t_{b})/2$. The
matrix aiming for capturing the node interactions,
$\mathbf{C}\mathrel{\mathop{\mathchar
58\relax}}=\mathbf{Q}\mathbf{Q}^{\top}\in\mathbb{R}^{N\times N}$ is defined
with a low-rank matrix $\mathbf{Q}\in\mathbb{R}^{N\times k}$ whose rows equal
to $1$ $(k\ll N)$, and we set $\mathbf{D}\mathrel{\mathop{\mathchar
58\relax}}=\mathbf{I}\ \mathbf{I}^{\top}\in\mathbb{R}^{D\times D}$. By
considering the Cholesky decomposition [36] of
$\mathbf{B}\mathrel{\mathop{\mathchar 58\relax}}=\mathbf{L}\mathbf{L}^{\top}$
since $\mathbf{B}$ is symmetric positive semi-definite, we can factorize
$\mathbf{K}\mathrel{\mathop{\mathchar
58\relax}}=\mathbf{K}_{f}\mathbf{K}_{f}^{\top}$ where
$\mathbf{K}_{f}\mathrel{\mathop{\mathchar
58\relax}}=\mathbf{L}\otimes\mathbf{Q}\otimes\mathbf{I}$.

Note that the precision matrix, $\mathbf{\Sigma}^{-1}$, can be written by
using the Woodbury matrix identity [36] as follows:

$\displaystyle\mathbf{\Sigma}^{-1}=\lambda^{-2}\left(\sigma_{\mathbf{\Sigma}}^{2}\mathbf{I}+\mathbf{K}_{f}\mathbf{K}_{f}^{\top}\right)^{-1}=\lambda^{-2}\left({\sigma_{\mathbf{\Sigma}}^{2}}^{-1}\mathbf{I}-{\sigma_{\mathbf{\Sigma}}^{2}}^{-1}\mathbf{K}_{f}\mathbf{R}^{-1}\mathbf{K}_{f}^{\top}{\sigma_{\mathbf{\Sigma}}^{2}}^{-1}\right)$

where the capacitance matrix $\mathbf{R}\mathrel{\mathop{\mathchar
58\relax}}=\mathbf{I}_{BKD}+{\sigma_{\mathbf{\Sigma}}^{2}}^{-1}\mathbf{K}_{f}^{\top}\mathbf{K}_{f}$.

The log-determinant of $\lambda^{2}\mathbf{\Sigma}$ can be also simplified by
applying Matrix Determinant lemma [36]:

$\displaystyle\log({det(\mathbf{\Sigma})})$
$\displaystyle=(BND)\log\big{(}\lambda^{2}\big{)}+\log\big{(}det\big{(}\sigma_{\mathbf{\Sigma}}^{2}\mathbf{I}_{BND}+\mathbf{K}_{f}\mathbf{K}_{f}^{\top}\big{)}\big{)}$
$\displaystyle=(BND)\log\big{(}\lambda^{2}\big{)}+\log\big{(}det\big{(}\mathbf{I}_{BKD}+{\sigma_{\mathbf{\Sigma}}^{2}}^{-1}\mathbf{K}_{f}^{\top}\mathbf{K}_{f}\big{)}\big{)}+(BND)\log\big{(}\sigma_{\mathbf{\Sigma}}^{2}\big{)}$
$\displaystyle=(BND)\big{(}\log(\lambda^{2})+\log\big{(}\sigma_{\mathbf{\Sigma}}^{2}\big{)}\big{)}+\log(det(\mathbf{R}))$

Note that the most cumbersome points in the computation of the prior are the
calculations of the inverse and determinant of the terms and some matrix
multiplication operations. Since $R$ is a matrix of size $BKD\times BKD$, its
inverse and determinant can be found in at most $\mathcal{O}(B^{3}K^{3}D^{3})$
operations. We also need the term, $\mathbf{K}_{f}\mathbf{R}^{-1}\mathbf{R}$,
which can also be computed in $\mathcal{O}(B^{3}D^{3}K^{2}|\mathcal{V}|)$
steps, so the number of operations required for the prior can be bounded by
$\mathcal{O}(B^{3}D^{3}K^{2}|\mathcal{V}|)$. It is worth noticing that we
cannot directly apply the batching strategy for the computation of the inverse
of the capacitance matrix, $\mathbf{R}$. However, we can compute it once and
then we can utilize it for the calculation of the log-prior for different sets
of node samples, then we can recompute it when we decide to update the
parameters again. To sum up, the complexity of our proposed approach is
$\mathcal{O}(B\mathcal{I}\mathcal{S}^{2}+B^{3}D^{3}K^{2}\mathcal{S}\mathcal{I})$
where $\mathcal{S}$ is the batch size and $\mathcal{I}$ is the number of
epochs.

Optimization of the proposed approach. Our objective given in Equation (7) is
not a convex function, thus the learning strategy that we follow is of great
importance in order to escape from local minima of poor quality
representations. We start by randomly initializing the model’s hyper-
parameters from $[-1,1]$ except for the velocity tensor, which is set to $0$
at the beginning. We adapt a sequential learning strategy for the learning of
these parameters. In other words, we first optimize the initial position and
bias terms together, $\\{\mathbf{x}^{(0)},\bm{\beta}\\}$, for $33$ epochs;
then, we include the velocity tensor, $\\{\mathbf{v}\\}$, into the
optimization process and repeat the training for the same number of epochs.
Finally, we add the prior parameters and learn all model hyper-parameters
together. We have employed Adam optimizer [35] with a learning rate of $0.1$.

In our experiments, we set the parameter $K=25$, and bins count $B=100$ to
have enough capacity to track node interactions. In order to find an optimal
regularization term $\lambda$ value and to determine the influence of the
prior in the objective, we apply an annealing strategy for the model. We first
mask $20\%$ of the dyads during the optimization of Equation (7). Furthermore,
we train the model by starting with $\lambda=10^{6}$ and learn all parameters
using the sequential optimization strategy. We then gradually reduce $\lambda$
to one-tenth upon optimizing all model parameters for $100$ epochs. The same
procedure is repeated until $\lambda=10^{-6}$. We choose the $\lambda$ value
minimizing the log-likelihood of the masked pairs (i.e., based on the
predictive log-likelihood evaluated on these pairs).

The final node embeddings are then obtained by performing this annealing
strategy without any mask until the found ideal $\lambda$ value. We repeat
this procedure $5$ times with different initializations, and we consider the
best-performing method/seed value in learning the final embeddings. The
relative standard deviation of the experiments is always less than $0.5$ for
all the networks, and we display the negative log-likelihood of the masked
pairs for the annealing strategy with $5$ random runs in Figure 5. The blue
curves demonstrate the same annealing strategy but in the opposite order. In
other words, we start from a very restrictive model with low $\lambda$ value
and increase $\lambda$ to have a more flexible model.

(a) Synthetic($\pi$)

(b) Synthetic($\mu$)

(c) College

(d) Contacts

(e) Email

(f) Forum

(g) Hypertext

Figure 5: Negative log-likelihood of the masked pairs for the annealing
strategy applied for tuning $\lambda$ parameter with $5$ random runs.

The considered annealing strategy thereby quantifies the impact for different
strengths of imposing the GP prior. It corresponds to a highly constrained
model akin to static representations for small values of $\lambda$ in which
the GP prior has close to zero variance of the parameters to highly flexible
dynamic representations almost entirely driven by the likelihood function for
high values of $\lambda$. The annealing strategy thus highlights the GP
prior’s impact and the optimal regime imposing such prior.

###  A.3 Theoretical Results

######  Lemma A.1.

For given fixed bias terms $\\{\beta_{i}\\}_{i\in\mathcal{V}}$, the node
embeddings, $\\{\mathbf{r}_{i}(t)\\}_{i\in\mathcal{V}}$, learned by the
objective given in Equation 1 within a bounded set of radius $R_{t}$ during a
time interval $[t_{l},t_{u}]$ satisfy

$\displaystyle\log\\!\left(\frac{(t_{u}\\!-t_{l})}{-\log
p_{ij}^{0}}\right)\\!+\\!(\beta_{i}\\!+\\!\beta_{j})\leq\frac{1}{(t_{u}\\!-\\!t_{l})}\\!\int_{t_{l}}^{t_{u}}\\!\\!\\!\\!||\mathbf{r}_{i}(t)\\!-\\!\mathbf{r}_{j}(t)||^{2}dt\\!\leq\\!\log\\!\left(\\!\frac{(t_{u}-t_{l})}{-\\!\log(1-p_{ij}^{>})}\\!\\!\right)\\!\\!+\\!(\beta_{i}\\!+\\!\beta_{j})\\!+\\!R_{t}$

where $p_{ij}^{0}$ and $p_{ij}^{>}$ are the probabilities of having zero and
more than zero links for the nodes $i$ and $j$.

###### Proof.

Let $p_{ij}^{0}\mathrel{\mathop{\mathchar
58\relax}}=\mathbb{P}\\{X_{ij}=0\\}$, and it can be written that

$\displaystyle p_{ij}^{0}=\exp(-\lambda_{ij})$
$\displaystyle=\exp\left(-\int_{t_{l}}^{t_{u}}\exp\left(\beta_{i}+\beta_{j}-||\mathbf{r}_{i}(t)-\mathbf{r}_{j}(t)||^{2}\right)dt\right)$
$\displaystyle=\exp\left(-\exp(\beta_{i}+\beta_{j})\int_{t_{l}}^{t_{u}}\exp\left(-||\mathbf{r}_{i}(t)-\mathbf{r}_{j}(t)||^{2}\right)dt\right)$

Then, the lower bound can be derived by

$\displaystyle\frac{-\log(p_{ij}^{0})}{(t_{u}-t_{l})}\frac{1}{\exp(\beta_{i}+\beta_{j})}$
$\displaystyle=\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}\exp\left(-||\mathbf{r}_{i}(t)-\mathbf{r}_{j}(t)||^{2}\right)dt$
$\displaystyle\geq\exp\left(-\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}||\mathbf{r}_{i}(t)-\mathbf{r}_{j}(t)||^{2}dt\right)$

where the second line follows from the Jensen’s inequality, and it can be
written that:

$\displaystyle\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}||\mathbf{r}_{i}(t)-\mathbf{r}_{j}(t)||^{2}dt$
$\displaystyle\geq\log\left(\frac{(t_{u}-t_{l})}{-\log
p_{ij}^{0}}\right)+(\beta_{i}+\beta_{j}).$

Based on a similar strategy that we have applied for the lower bound, we can
write

$\displaystyle\frac{-\log(1-p_{ij}^{>})}{(t_{u}-t_{l})}\frac{1}{\exp(\beta_{i}+\beta_{j})}$
$\displaystyle=\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}\exp\left(-||\mathbf{r}_{i}(t)-\mathbf{r}_{j}(t)||^{2}\right)dt$
$\displaystyle\leq\exp\left(R_{t}\right)\exp\left(\\!\\!-\frac{1}{(t_{u}\\!-\\!t_{l})}\int_{t_{l}}^{t_{u}}\\!\\!\\!\\!||\mathbf{r}_{i}(t)\\!-\\!\mathbf{r}_{j}(t)||^{2}dt\right)$

where the inequality follows from Lemma A.2 in which $f(t)$ corresponds to
$||\mathbf{r}_{i}(t)-\mathbf{r}_{j}(t)||^{2}$ and
$R_{t}\mathrel{\mathop{\mathchar
58\relax}}=\sup_{t\in[t_{l},t_{u}]}||\mathbf{r}_{i}(t)-\mathbf{r}_{j}(t)||^{2}$.
Then, it can be concluded that

$\displaystyle\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}||\mathbf{r}_{i}(t)-\mathbf{r}_{j}(t)||^{2}dt$
$\displaystyle\leq\log\left(\exp(\beta_{i}+\beta_{j})\exp\left(R_{t}\right)\frac{(t_{u}-t_{l})}{-\log(1-p_{ij}^{>})}\right)$
$\displaystyle\leq(\beta_{i}+\beta_{j})+R_{t}+\log\left(\frac{(t_{u}-t_{l})}{-\log(1-p_{ij}^{>})}\right)$

∎

######  Lemma A.2.

Let $f\mathrel{\mathop{\mathchar
58\relax}}[t_{l},t_{u}]\rightarrow\mathbb{R}^{+}$ be a real-valued non-
negative continuous function, then it satisfies

$\displaystyle\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}\exp(-f(t))dt\leq\exp(R_{t})\exp\left(-\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}f(t)dt\right)$

for some constant $R_{t}\mathrel{\mathop{\mathchar
58\relax}}=\sup_{t\in[t_{l},t_{u}]}f(t)\in\mathbb{R}$.

###### Proof.

Let $\hat{f}\mathrel{\mathop{\mathchar
58\relax}}=\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}f(t)dt$ for some
$t_{u},t_{l}\in\mathbb{R}$ where $t_{l}<t_{u}$, then it can be written that

$\displaystyle\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}\exp(-f(t)+\hat{f})dt$
$\displaystyle\leq\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}\exp(\hat{f})dt$
$\displaystyle=\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}\exp\left(\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}f(t^{{}^{\prime}})dt^{{}^{\prime}}\right)dt$
$\displaystyle\leq\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}\exp\left(\sup_{t^{*}\in[t_{l},t_{u}]}f(t^{*})\right)dt$
$\displaystyle=\exp\left(R_{t}\right)$

where $R_{t}\mathrel{\mathop{\mathchar
58\relax}}=\sup_{t\in[t_{l},t_{u}]}f(t)$, and we have the first inequality
since $f(t)$ is a non-negative function. In conclusion, we can obtain

$\displaystyle\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}\exp(-f(t))dt\leq\exp\left(R_{t}\right)\exp\left(-\frac{1}{(t_{u}-t_{l})}\int_{t_{l}}^{t_{u}}f(t)dt\right)$

∎

######  Theorem A.3.

Let $\mathbf{f}(t)\mathrel{\mathop{\mathchar
58\relax}}[0,T]\rightarrow\mathbb{R}^{D}$ be a continuous embedding of a node.
For any given $\epsilon>0$, there exists a continuous, piecewise-linear node
embedding, $\mathbf{r}(t)$, satisfying
$||\mathbf{f}(t)-\mathbf{r}(t)||_{2}<\epsilon$ for all $t\in[0,T]$ where
$\mathbf{r}(t)\mathrel{\mathop{\mathchar 58\relax}}=\mathbf{r}^{(b)}(t)$ for
all $(b-1)\Delta_{B}\leq t<b\Delta_{B}$,
$\mathbf{r}(t)\mathrel{\mathop{\mathchar 58\relax}}=\mathbf{r}^{(B)}(t)$ for
$t=T$ and $\Delta_{B}=T/B$ for some $B\in\mathbb{N}^{+}$.

###### Proof.

Let $\mathbf{f}(t)\mathrel{\mathop{\mathchar
58\relax}}[0,T]\rightarrow\mathbb{R}^{D}$ be a continuous embedding so it is
also uniformly continuous by the Heine–Cantor theorem since $[0,T]$ is a
compact set. Then, we can find some $B\in\mathbb{N}^{+}$ such that for every
$t,\tilde{t}\in[0,T]$ with
$|t-\tilde{t}|\leq\Delta_{B}\mathrel{\mathop{\mathchar 58\relax}}=T/B$ implies
$||\mathbf{f}(t)-\mathbf{f}(\tilde{t})||_{2}<\epsilon/2$ for any given
$\epsilon>0$.

Let us define
$\mathbf{r}^{(b)}(t)=\mathbf{r}^{(b-1)}\left((b-1)\Delta_{B}\right)+\mathbf{v}_{b}(t-(b-1)\Delta_{B})$
recursively for each $b\in\\{1,\ldots,B\\}$ where
$\mathbf{r}^{(0)}(0)\mathrel{\mathop{\mathchar
58\relax}}=\mathbf{x}_{0}=\mathbf{f}(0)$, and
$\mathbf{v}_{b}\mathrel{\mathop{\mathchar
58\relax}}=\frac{\mathbf{f}\left(b\Delta_{B}\right)-\mathbf{f}\left((b-1)\Delta_{B}\right)}{\Delta_{B}}$.
Then it can be seen that we have
$\mathbf{r}^{(b)}\left(b\Delta_{B}\right)=\mathbf{f}(b\Delta_{B})$ for all
$b\in\\{1,\ldots,B\\}$ because

$\displaystyle\mathbf{r}^{(b)}(b\Delta_{B})$
$\displaystyle=\mathbf{r}^{(b-1)}\left((b-1)\Delta_{B}\right)+\mathbf{v}_{b}\left(b\Delta_{B}-\Delta_{B}(b-1)\right)$
$\displaystyle=\mathbf{r}^{(b-1)}\left((b-1)\Delta_{B}\right)+\mathbf{v}_{b}\Delta_{B}$
$\displaystyle=\mathbf{r}^{(b-1)}\left((b-1)\Delta_{B}\right)+\left(\frac{\mathbf{f}\left(b\Delta_{B}\right)-\mathbf{f}\left((b-1)\Delta_{B}\right)}{\Delta_{B}}\right)\Delta_{B}$
$\displaystyle=\mathbf{r}^{(b-1)}\left((b-1)\Delta_{B}\right)+\left(\mathbf{f}\left(b\Delta_{B}\right)-\mathbf{f}\left((b-1)\Delta_{B}\right)\right)$
$\displaystyle=\mathbf{r}^{(b-2)}\left((b-2)\Delta_{B}\right)+\left(\mathbf{f}\left((b-1)\Delta_{B}\right)-\mathbf{f}((b-2)\Delta_{B})\right)+\left(\mathbf{f}\left(b\Delta_{B}\right)-\mathbf{f}\left((b-1)\Delta_{B}\right)\right)$
$\displaystyle=\mathbf{r}^{(b-2)}\left((b-2)\Delta_{B}\right)+\left(\mathbf{f}\left(b\Delta_{B}\right)-\mathbf{f}((b-2)\Delta_{B})\right)$
$\displaystyle=\cdots$
$\displaystyle=\mathbf{r}^{(0)}(0)+\left(\mathbf{f}\left(b\Delta_{B}\right)-\mathbf{f}(0)\right)$
$\displaystyle=\mathbf{f}\left(b\Delta_{B}\right)$

where the last line follows from the fact that
$\mathbf{r}^{(0)}(0)=\mathbf{x}_{0}=\mathbf{f}(0)$ by the definition.
Therefore, for any given point $t\in[0,T)$ for $b=\lfloor
t/\Delta_{b}\rfloor+1$, it can be seen that

$\displaystyle||\mathbf{f}(t)-\mathbf{r}(t)||_{2}$
$\displaystyle=||\mathbf{f}(t)-\mathbf{r}^{(b)}(t)||_{2}$
$\displaystyle=\Big{|}\Big{|}\mathbf{f}(t)-\Big{(}\mathbf{r}^{(b-1)}((b-1)\Delta_{B})+\mathbf{v}_{b}(t-(b-1)\Delta_{B})\Big{)}\Big{|}\Big{|}_{2}$
$\displaystyle=\Bigg{|}\Bigg{|}\mathbf{f}(t)-\Bigg{(}\mathbf{r}^{(b-1)}\left((b-1)\Delta_{B}\right)+\left(\frac{\mathbf{f}(b\Delta_{B})-\mathbf{f}\left((b-1)\Delta_{B}\right)}{\Delta_{B}}\right)(t-(b-1)\Delta_{B})\Bigg{)}\Bigg{|}\Bigg{|}_{2}$
$\displaystyle=\Big{|}\Big{|}\left(\mathbf{f}(t)-\mathbf{r}^{(b-1)}\left((b-1)\Delta_{B}\right)\right)+\left(\mathbf{f}(b\Delta_{B})-\mathbf{f}\left((b-1)\Delta_{B}\right)\right)\left(\frac{t-(b-1)\Delta_{B}}{\Delta_{B}}\right)\Big{|}\Big{|}_{2}$
$\displaystyle\leq\Big{|}\Big{|}\mathbf{f}(t)-\mathbf{r}^{(b-1)}\left((b-1)\Delta_{B}\right)\Big{|}\Big{|}+\Big{|}\Big{|}\mathbf{f}(b\Delta_{B})-\mathbf{f}\left((b-1)\Delta_{B}\right)\Big{|}\Big{|}$
$\displaystyle<\frac{\epsilon}{2}+\frac{\epsilon}{2}$ $\displaystyle=\epsilon$

where the inequality in the fifth line holds since we have
$\left|\frac{t-(b-1)\Delta_{B}}{\Delta_{B}}\right|\leq 1$ ∎

######  Lemma A.4 (Integral Computation).

The integral of the intensity function, $\lambda_{ij}(t)$, from $t_{l}$ to
$t_{u}$ is equal to

$\displaystyle\int_{t_{l}}^{t_{u}}\exp\left(\beta_{ij}-\mathinner{\\!\left\lVert\Delta\mathbf{x}_{ij}+\Delta\mathbf{v}_{ij}t\right\rVert}^{2}\right)=\frac{\sqrt{\pi}\exp{\left(\beta_{ij}+r_{ij}^{2}-\mathinner{\\!\left\lVert\Delta\mathbf{x}_{ij}\right\rVert}^{2}\right)}}{2\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}}\mathrm{erf}\Big{(}\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}t+r_{ij}\Big{)}\Bigg{|}_{t=t_{l}}^{t=t_{u}}$

where $\beta_{ij}\mathrel{\mathop{\mathchar 58\relax}}=\beta_{i}+\beta_{j}$,
$\Delta\mathbf{x}_{ij}\mathrel{\mathop{\mathchar
58\relax}}=\mathbf{x}_{i}^{(0)}-\mathbf{x}_{j}^{(0)}$,
$\Delta\mathbf{v}_{ij}\mathrel{\mathop{\mathchar
58\relax}}=\mathbf{v}_{i}^{(1)}-\mathbf{v}_{j}^{(1)}$ and
$r\mathrel{\mathop{\mathchar
58\relax}}=\frac{\langle\Delta\mathbf{v}_{ij},\Delta\mathbf{x}_{ij}\rangle}{\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}}$.

###### Proof.

$\displaystyle\int_{t_{l}}^{t_{u}}\exp\left(-\mathinner{\\!\left\lVert\Delta\mathbf{x}_{ij}+\Delta\mathbf{v}_{ij}t\right\rVert}^{2}\right)$
$\displaystyle=\int_{t_{l}}^{t_{u}}\exp\left(-\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}^{2}t^{2}-2\left\langle\Delta\mathbf{x}_{ij},\Delta\mathbf{v}_{ij}\right\rangle
t-\mathinner{\\!\left\lVert\Delta\mathbf{x}_{ij}\right\rVert}^{2}\right)\mathrm{d}t$
$\displaystyle=\int_{t_{l}}^{t_{u}}\exp\left(-\left(\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}t+r_{ij}\right)^{2}+r_{ij}^{2}-\mathinner{\\!\left\lVert\Delta\mathbf{x}_{ij}\right\rVert}^{2}\right)\mathrm{d}t$
(8)

where $r_{ij}\mathrel{\mathop{\mathchar
58\relax}}=\frac{\langle\Delta\mathbf{v}_{ij},\Delta\mathbf{x}_{ij}\rangle}{\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}}$.
The substitution
$u=\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}t+r_{ij}$ yields
$\mathrm{d}u=\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}\mathrm{d}t$.
Furthermore, we have

$\displaystyle\int_{t_{l}}^{t_{u}}\\!\\!\\!\\!\\!\exp\left(-\left(\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}t\\!+\\!r_{ij}\right)^{2}\right)\mathrm{d}t$
$\displaystyle=\frac{1}{\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}}\int_{\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}t_{l}+r_{ij}}^{\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}t_{u}+r_{ij}}\exp\left(-u^{2}\right)\\!\mathrm{d}u$
$\displaystyle=\frac{1}{\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}}\frac{\sqrt{\pi}}{2}\left(\frac{2}{\sqrt{\pi}}\int_{\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}t_{l}+r_{ij}}^{\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}t_{u}+r_{ij}}\exp\left(-u^{2}\right)\mathrm{d}u\right)$
$\displaystyle=\frac{\sqrt{\pi}}{2\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}}\mathrm{erf}\Big{(}\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}t+r_{ij}\Big{)}\Bigg{|}_{t=t_{l}}^{t=t_{u}}$
(9)

By using Equations 8 and 9, it can be obtained that

$\displaystyle\int_{t_{l}}^{t_{u}}\exp\left(-\mathinner{\\!\left\lVert\Delta\mathbf{x}_{ij}+\Delta\mathbf{v}_{ij}t\right\rVert}^{2}\right)$
$\displaystyle=\exp{\left(r_{ij}^{2}-\mathinner{\\!\left\lVert\Delta\mathbf{x}_{ij}\right\rVert}^{2}\right)}\int_{t_{l}}^{t_{u}}\exp\left(-\left(\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}t+r_{ij}\right)^{2}\right)\mathrm{d}t$
$\displaystyle=\frac{\sqrt{\pi}\exp{\left(r_{ij}^{2}-\mathinner{\\!\left\lVert\Delta\mathbf{x}_{ij}\right\rVert}^{2}\right)}}{2\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}}\Bigg{[}\mathrm{erf}\Big{(}\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}t+r_{ij}\Big{)}\Bigg{|}_{t=t_{l}}^{t=t_{u}}$

Therefore, we can conclude that

$\displaystyle\int_{t_{l}}^{t_{u}}\exp\left(\beta_{ij}-\mathinner{\\!\left\lVert\Delta\mathbf{x}_{ij}+\Delta\mathbf{v}_{ij}t\right\rVert}^{2}\right)=\frac{\sqrt{\pi}\exp{\left(\beta_{ij}+r_{ij}^{2}-\mathinner{\\!\left\lVert\Delta\mathbf{x}_{ij}\right\rVert}^{2}\right)}}{2\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}}\mathrm{erf}\Big{(}\mathinner{\\!\left\lVert\Delta\mathbf{v}_{ij}\right\rVert}t+r_{ij}\Big{)}\Bigg{|}_{t=t_{l}}^{t=t_{u}}$

∎

###  A.4 Extension to Weighted, Directed, and Bipartite Networks

In this section, we will discuss how the proposed approach, PiVeM, can be
extended for weighted, directed, and bipartite networks.

Weighted networks. Our approach can be simply adapted for positive integer-
weighted networks by replacing each weighted link in the network with multiple
unit events corresponding to the integer weight at the specific time point of
the integer-weighted link. Then, PiVeM can be run as is for the reinterpreted
version of the network without making any modifications to the structure of
the approach.

Directed and Bipartite networks. Let
$\mathcal{G}=((\mathcal{V},\mathcal{U}),\mathcal{E})$ be a bipartite network
with the parts $\mathcal{V}=\\{v_{1},\ldots,v_{N_{1}}\\}$ and
$\mathcal{U}=\\{u_{1},\ldots,u_{N_{2}}\\}$. We can rewrite the objective given
in Equation 7 by considering only the pairs
$(i,j)\in\mathcal{V}\times\mathcal{U}$ belonging to different parts:

$\displaystyle\hat{\Omega}=\operatorname*{arg\,max}_{\Omega}\sum_{i\in\mathcal{V}}\sum_{j\in\mathcal{U}}\left(\sum_{e_{ij\in\mathcal{E}_{ij}}}\log\lambda_{ij}(e_{ij})-\int_{0}^{T}\lambda_{ij}(t)dt\right)+\log\mathcal{N}\left(\left[\begin{array}[]{c}\mathbf{x}^{(0)}\\\
\mathbf{v}\end{array}\right];\mathbf{0},\bm{\Sigma}\right)$ (12)

where the intensity function, $\lambda_{ij}(e_{ij})$ is defined as follows:

$\displaystyle\lambda_{ij}(t)\mathrel{\mathop{\mathchar
58\relax}}=\exp{\left(\beta_{i}+\eta_{j}-||\mathbf{r}^{*}_{i}(t)-\mathbf{r}^{**}_{j}(t)||^{2}\right),}$
(13)

where $\beta_{i}$ and $\eta_{j}$ indicate the bias/random effect terms for the
nodes belonging respectively to $\mathcal{V}$ an $\mathcal{U}$. Similarly, we
can introduce distinct initial position $\mathbf{x}^{*}_{i}$,
$\mathbf{x}^{**}_{j}$ and velocity tensors $\mathbf{v}^{*}_{i}$,
$\mathbf{v}^{**}_{j}$ to define the node representations,
$\mathbf{r}^{*}_{i}(t)\in\mathbb{R}^{D}$ and
$\mathbf{r}^{**}_{j}(t)\in\mathbb{R}^{D}$ at time $t$. Note that we can write
$\mathbf{x}_{i}^{(0)}=\mathbf{x}^{*}_{i}\oplus\mathbf{x}^{**}_{i}$, and
$\mathbf{v}_{i}=\mathbf{v}^{*}_{i}\oplus\mathbf{v}^{**}_{i}$ where $\oplus$
indicates the tensor concatenation operation.

For the directed case, we specify the model similar to the bipartite case but
define the likelihood function as

$\displaystyle\hat{\Omega}=\operatorname*{arg\,max}_{\Omega}\sum_{i\neq
j}\left(\sum_{e_{ij\in\mathcal{E}_{ij}}}\log\lambda_{ij}(e_{ij})-\int_{0}^{T}\lambda_{ij}(t)dt\right)+\log\mathcal{N}\left(\left[\begin{array}[]{c}\mathbf{x}^{(0)}\\\
\mathbf{v}\end{array}\right];\mathbf{0},\bm{\Sigma}\right).$ (16)

###  A.5 Table of Symbols

The detailed list of the symbols used throughout the manuscript and their
corresponding definitions can be found in Table 5.

Table 5: Table of symbols Symbol | Description
---|---
$\mathcal{G}$ | Graph
$\mathcal{V}$ | Vertex set
$\mathcal{E}$ | Edge set
$\mathcal{E}_{ij}$ | Edge set of node pair $(i,j)$
$N$ | Number of nodes
$D$ | Dimension size
$\mathcal{I}_{T}$ | Time interval
$T$ | Time length
$B$ | Number of bins
$\beta_{i}$ | Bias term of node $i$
$\mathbf{x}$ | Initial position matrix
$\mathbf{v}^{(b)}$ | Velocity matrix for bin $b$
$\mathbf{r}_{i}(t)$ | Position of node $i$ at time $t$
$\lambda_{ij}(t)$ | Intensity of node pair $(i,j)$ at time $t$
$e_{ij}$ | An event time of node pair $(i,j)$
$\Sigma$ | Covariance matrix
$\lambda$ | Scaling factor of the covariance
$\sigma_{\mathbf{\Sigma}}$ | Noise variance
$\sigma_{\mathbf{B}}$ | Lengthscale variable of RBF kernel
$\otimes$ | Kronecker product
$\mathbf{I}$ | Identity matrix
$\mathbf{B}$ | Bin interaction matrix
$\mathbf{C}$ | Node interaction matrix
$\mathbf{D}$ | Dimension interaction matrix
$\mathbf{R}$ | Capacitance matrix
$K$ | Latent dimension of $\mathbf{C}$