Continuous-time edge modelling using non-parametric point processes
Year: 2022

Abstract:
The mutually-exciting Hawkes process (ME-HP) is a natural choice to model reciprocity, which is an important attribute of continuous-time edge (dyadic) data. However, existing ways of implementing the ME-HP for such data are either inflexible, as the exogenous (background) rate functions are typically constant and the endogenous (excitation) rate functions are specified parametrically, or inefficient, as inference usually relies on Markov chain Monte Carlo methods with high computational costs. To address these limitations, we discuss various approaches to model design, and develop three variants of non-parametric point processes for continuous-time edge modelling (CTEM). The resulting models are highly adaptable as they generate intensity functions through sigmoidal Gaussian processes, and so provide greater modelling flexibility than parametric forms. The models are implemented via a fast variational inference method enabled by a novel edge modelling construction. The superior performance of the proposed CTEM models is demonstrated through extensive experimental evaluations on four real-world continuous-time edge data sets.

heading: 1 Introduction
text: Continuous-time edge (or dyadic) data, which comprise information on edges (between two individuals or nodes) and their times of occurrence, arise naturally in many applications, such as sending and receiving messages, posting and answering comments on social networks, and purchasing groceries in online stores. Reciprocity is an important attribute of edge data, in that many edges are induced as responses to previous oppositely-directed edges (e.g. a reply email is written in response to an initial message between two individuals). As a result, mutually-exciting Hawkes processes (MEHP) [13, 12], which regard the directed edges as events, are commonly used to analyse edge data [3, 8, 25, 29, 26, 20]. For example, [3] uses the ME-HP to model the reciprocity of group-wise edges, and [20] uses sparse exchangeable random measures [4, 28] to model the ME-HPs exogenous values.
However, these existing methods are either limited in modelling flexibility, or inefficient to implement. For the former, the exogenous (background rate) functions in the ME-HP are typically assumed to be static and the endogenous (excitation) functions are specified parametrically, which may be inadequate for modelling complex real-world systems. E.g., the exogenous intensities may exhibit daily or weekly periodicity, and the endogenous functions might not be strictly decreasing over time. For the latter, Markov chain Monte Carlo (MCMC) methods are usually adopted for model inference, which usually require high computational costs and large output storage. Without proper scalable inference methods, it is difficult to deploy these models to larger-scale datasets.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
In this work, to address the above limitations, we develop continuous-time edge models (CTEM) based on non-parametric point processes. The CTEM models also use point processes to model data specific to nodal pairs. CTEMs have two major differences compared to existing approaches: Firstly, they place sigmoidal Gaussian process priors on the individual exogenous functions and endogenous functions. As a result, the exogenous functions are able to evolve over time and account for gradual changes and variations in the background environment (such as, e.g., trends or periodicity). Similarly, non-parametric modelling of the endogenous functions enables more behaviour choices (e.g., the rate of decrease in the excitation functions may vary over time, or there might be increases in mutual-excitation effects). Secondly, we develop an efficient variational inference method for model inference, which is more suited for large-scale data modelling than MCMC methods.
By specifying K exogenous functions and K endogenous functions, we propose a family of three non-parametric CTEM models on N nodes, in which the directed edges from node i to node j are regarded as events for the (i, j)-th point process:
• Nodal pair-specific Poisson process (NP-PP): For each directed nodal pair (i, j), i, j ∈ {1, . . . , N}, the intensity function for their Poisson process is constructed from the linear combination of K exogenous functions and nodes i, j’s latent exogenous features.
• Nodal pair-specific mutually-exciting Hawkes process (NP-HP): We additionally include the combination of K endogenous functions and nodes i, j’s latent endogenous features as the endogenous functions for the (i, j)-th ME-HP.
• Mutually-exciting Hawkes process with edge clustering (NP-HP-C): Rather than using nodal pair-specific endogenous functions, we associate each edge with a latent label, and assume that those edges with the same label are allocated to the same endogenous function.
The last approach may be used to model the phenomenon that many edges are often associated with certain behavioural or other patterns. For instance, messages may be grouped into particular labelled events, or edges associated with different patterns may generate different triggering effects.
Inspired by inference techniques [5, 31] which adopt the Polya-Gamma trick to form conjugate models for sigmoidal Gaussian process-modulated Poisson processes and Hawkes processes, we develop mean-field variational algorithms to enable fast inference for the proposed models. In addition we also construct a fast variational inference solution to the Dirichlet-Hawkes process related models [7, 19, 24], which typically heavily rely on Sequential Monte Carlo [6] methods.
In summary, the merits of the proposed CTEM family include: (1) Modelling flexibility: The proposed models are highly flexible as they are able to describe both exogenous and endogenous functions using non-parametric sigmoidal Gaussian processes. The experimental results on four real-world data sets show clear periodic patterns in both exogenous and endogenous functions. (2) Inferential efficiency: To the best of our knowledge, this is the first work to develop fast variational inference for the continuous-time edge modelling problem, making it suitable for large-scale relational modelling settings; (3) Modelling choice exploration: The CTEM family permits different structural choices in modelling continuous-time edge data, including nodal pair-specific intensity functions and latent label determined endogenous functions.

heading: 2 Background on mutually-exciting Hawkes processes for edge data
text: The mutually-exciting Hawkes process [12, 13] is a special case of the temporal point process (TPP). TPPs are realised through sets of events in a time interval [0, T ], and are parameterised by an intensity function λ(t), which states that the probability of an event occurring in any infinitesimal time interval [t, t+ dt) is λ(t)dt. TPPs which are independent of the history of the process (e.g. a Poisson point process) are not well-suited to describe edge data (which comprise directed events from one node to another), as the occurrence of edges may be the result of responses to previous oppositely-directed edges. For example, receiving a message from a person would likely increase the probability of sending a response message to that person over the following time period.
The ME-HP provides a natural way to model the reciprocity of edge data. Its intensity functions λij(t), which describes edge generation from node i to node j at time t, depend on the historical edges from node j to node i. We use the triplet (im, jm, tm) to denote an edge, meaning that among all edges, the m-th one is initiated from node im to node jm and occurred at time tm. Using fij(t)
and gij(t) as the exogenous function and endogenous function from node i to node j respectively, the intensity function λij(t) and the likelihood L(·) of the ME-HP over all M edges is
λij(t) = fij(t) + ∑
m′:im′=j∩jm′=i gij(t− tm′), (1)
L({λij(t)}(i,j):i̸=j |{(im, jm, tm)}Mm=1|) = ∏
(i,j):i ̸=j e− ∫ T0 λij(t)dt ∏ m′:im′=i∩jm′=j λij(tm′).  (2)
ME-HP with block models: Stochastic block models [22, 21, 14] are a classical way to describe static edge data, via partitioned nodal groups and group-wise edges. Several methods [3, 8, 25, 29, 26, 20] have incorporated stochastic block models within the ME-HP to describe edge data reciprocity. Among these, the HPGP-IRM [25] is the most similar model to ours, with the generative process
ν ∼ CRP(1 : N), λij(t) = γν(i)ν(j) + ∑
m′:im′=j∩jm′=i βij(Fi, Fj)gij(t− tm′) (3)
where ν is a Chinese Restaurant Process (CRP) [23] group partition on all N nodes, γν(i)ν(j) is the exogenous intensity from person i’s group ν(i) to person j’s group ν(j), βij(Fi, Fj) is the endogenous scalar measuring the feature similarity between person i and person j through their features Fi, Fj , and the occurrence times {tm}m follow a Hawkes process (HP) with intensity function λij(t). Similar to other approaches, the HPGP-IRM specifies parametric functions for both exogenous and endogenous functions, and its MCMC-based approach to inference is not practical in large-scale settings. Further, the absence of edge clustering capability means that intrinsic structure among persons, especially the sharing of endogenous functions, can not be explored.

heading: 3 Non-parametric point processes for continuous-time edge modelling
text: The family of Continuous-Time Edge Modelling (CTEM) methods is inspired by three motivations: constructing highly flexible models to adequately describe the complexities of real world edge data; enabling fast variational inference; and exploring different modelling choices for studying continuoustime edge data. The following sections introduce three non-parametric point processes, in which the intensity functions are generated from sigmoidal Gaussian processes.

heading: 3.1 Nodal pair-specific non-parametric Poisson processes (NP-PP)
text: The first approach (NP-PP) is to use nodal pair-specific Poisson processes to model the continuoustime edge data. This involves K exogenous functions σ(fk(t)), where fk(t) is a random function generated from a Gaussian process and σ(x) = 1/(1 + exp(−x)) is the sigmoidal function to ensure the positiveness of the intensity values. Each node i is assumed to have K-length positive latent exogenous features πi ∈ [R+]1×K , in which the k-th entry πik represents node i’s affiliation with the k-th exogenous feature. The intensity function of the Poisson process for generating edges from node i to node j is thus a linear combination of {σ(fk(t))}k weighted by πi,πj , with the corresponding generative process described as follows:
1. Generate fk(t) ∼ GP(0, κ(f)k ) for group features k = 1, . . . ,K, where GP(0, κ) denotes a Gaussian process with mean 0 and covariance function κ(·, ·).
2. Generate {πik}Ni=1 ∼ Gam(aπ, bπ) for group features k = 1, . . . ,K, where Gam(a, b) denotes a Gamma distribution with mean a/b and variance a/b2.
3. Generate the time stamps of edges from node i to node j as: {(i, j, tm′)}m′ ∼ PoissonProcess(λij(t)), with intensity function λij(t) = ∑ k πikπjk · σ(fk(t)).
Under this setting, the intensity functions between any nodal pair are symmetric λij(t) = λji(t). Together, πik and πjk scale the k-th exogenous intensity onto the interval (0, πikπjk). Nodes with similar feature-affiliation values are more likely to generate edges. The L1-norm of the latent features of a node πi indicates its edge activity, since a larger value means it is more likely to induce edges.
This non-parametric setting provides greater modelling flexibility than those using static or piecewiseconstant exogenous functions. I.e., the NP-PP allows each exogenous function to change over time, while most existing approaches [3, 8, 25, 29, 26, 20] use a time-constant exogenous intensity value, even though real-world settings are known to vary in this manner. [8] also used Poisson processes to model the edges, however, their intensity functions are piecewise-constant with the changepoints and intensty values almost pre-fixed, which is substantially less flexible than our approach.
However, without incorporating endogenous functions, the NP-PP cannot describe the reciprocity property of continuous-time edges. A direct consequence might be overestimated exogenous intensities, as they need to explain all edges, and consequently lower model performance.

heading: 3.2 Nodal pair-specific non-parametric Hawkes process (NP-HP)
text: Our second approach is to use mutually-exciting Hawkes processes (ME-HP) to model the continuoustime edge data. In addition to the non-parametric exogenous functions in the first approach, the NP-HP includes nodal pair-specific endogenous functions to model the reciprocity of edge data. In particular, the endogenous function for the (i, j)-th nodal pair is structured as ∑ k ρikρjk · σ(gk(t)), where {σ(gk(t))}k are random functions from a sigmoidal Gaussian process and ρi ∈ [R+]1×K is node i’s K-length latent endogenous features. As a result, the occurrence of an edge at time tm is either instantiated by the exogenous effect (i.e. ∑ k πikπjk · σ(fk(tm))) or triggered by one of its
historical oppositely-directed edges (through the endogenous effect ∑
k ρikρjk · σ(gk(tm − tm′)), where m′ ∈ Hij(tm) and Hij(t) = {m′ : (tm′ < t) ∩ (im′ , jm′) = (j, i)} denotes the set of historical (i.e. occurring before time stamp t) oppositely-directed (i.e. from node j to node i) edges). Since the occurrence of an edge increases the intensity function value of its future oppositely-directed edges, the NP-HP can describe reciprocity well. It’s generative process is described as follows:
1. Generate fk(t) ∼ GP(0, κ(f)k ), gk(t) ∼ GP(0, κ (g) k ), k = 1, . . . ,K. 2. Generate πik ∼ Gam(aπ, bπ), ρik ∼ Gam(aρ, bρ), i = 1, . . . , N, k = 1, . . . ,K. 3. Simultaneously generate two sets of triplets {(i, j, tm′)}m′ and {(j, i, tm′)}m′ , which are
edges initiated from node i to node j, and from node j to node i respectively, from Hawkes processes with the following intensity functions:
λij(t) = ∑ k πikπjk · σ(fk(t)) + ∑ k ∑ m′:m′∈Hij(t) ρikρjk · σ(gk(t− tm′)), (4)
λji(t) = ∑ k πjkπik · σ(fk(t)) + ∑ k ∑ m′:m′∈Hji(t) ρjkρik · σ(gk(t− tm′)). (5)
Since the functions {σ(gk(t))}k are also drawn from sigmoidal Gaussian process, the NP-HP’s endogenous functions possess greater flexibility than previous approaches, in which the endogenous functions are typically specified as exponential functions. Incorporating the product of the nodes’ latent endogenous features ensures that these endogenous functions are nodal pair specific, as for the exogenous functions. However, NP-HP assumes that each nodal pair would use same endogenous functions for all its edges. This might be unable to describe the changes in the edges, such as emails with different topics between two persons might generate distinct endogenous effects.

heading: 3.3 Non-parametric Hawkes process with edge-clustering effects (NP-HP-C)
text: The third approach introduces an edge clustering mechanism for the mutually-exciting Hawkes process with nodal pair-specific exogenous functions. Here each generated edge is assigned a latent group label. The NP-HP-C generative process for edges in the time period [0, T ] is then as follows:
1. Generate fk(t) ∼ GP(0, κ(f)k ), gk(t) ∼ GP(0, κ (g) k ), k = 1, . . . ,K. 2. Generate πik ∼ Gam(aπ, bπ), k = 1, . . . ,K, i = 1, . . . , N ; 3. Generate vk ∼ Gam(av, bv), k = 1, . . . ,K, where vk is the scaling variable for the k-th
endogenous function σ(gk(t)). 4. For each pair of nodes {(i, j) : 1 ≤ i < j ≤ N}, generate triplets
{(i, j, t(m ′) ij ), b (m′) ij , z (m′) ij }m′ and {(j, i, t (m′′) ji ), b (m′′) ji , z (m′′) ji }m′′ for the two directed edges
according to the following procedure, until tM ′ ≤ T < tM ′+1, tM ′′ ≤ T < tM ′′+1. Initialise the current time stamp at t∗ = 0 and set Mij = Mji = 0. Then:
(a) Starting from t∗, obtain two candidate event times t̄ij = min({t̄(l)ij }l|t̄ (l) ij > t∗), t̄ji =
min({t̄(l)ji }l|t̄ (l) ji > t∗), where the two time sequences {t̄ (l) ij }l, {t̄ (l) ji }l are generated from:
{t̄(l)ij }l ∼ PoissonProcess( K∑
k=1
πikπjkσ(fk(t)) + Mji∑ m′=1 v z (m′) ji · σ(g z (m′) ji (t− t(m ′) ji ))), (6)
{t̄(l)ji }l ∼ PoissonProcess( K∑
k=1
πikπjkσ(fk(t)) + Mij∑ m′=1 v z (m′) ij · σ(g z (m′) ij (t− t(m ′) ij ))). (7)
(Note that the second term in Eqs (6) and (7) is not involved when Mij = Mji = 0). If t̄ij > T ∩ t̄ji > T , break for this nodal pair and start generating edges for other pairs.
(b) If t̄ij < t̄ji, set Mij = Mij + 1, t∗ = t (Mij) ij = t̄ij and generate b (Mij) ij , z (Mij) ij
from (m′ = 1, . . . ,Mji)
P (b (Mij) ij , z (Mij) ij ) ∝  πikπjkσ(fk(t∗)), if b (Mij) ij = 0, z (Mij) ij = k;
v z (m′) ji σ(g z (m′) ji
(t∗ − t(m ′) ji )), if b (Mij) ij = m ′, z (Mij) ij = z (m′) ji ,
else, set Mji = Mji + 1, t∗ = t (Mji) ji = t̄ji and generate b (Mji) ji , z (Mji) ji from (m ′ = 1, . . . ,Mij)
P (b (Mji) ji , z (Mji) ji ) ∝  πikπjkσ(fk(t∗)), if b (Mji) ji = 0, z (Mji) ji = k;
v z (m′) ij σ(g z (m′) ij
(t∗ − t(m ′) ij )), if b (Mji) ji = m ′, z (Mji) ji = z (m′) ij .
where Mij is the number of edges from node i to node j, and where t (m) ij , b (m) ij ∈ {{0}∪Hij(t (m) ij )}, and z(m)ij ∈ {1, . . . ,K} denote the time, “source” factor (i.e., exogenous rate or endogenous rates) and latent label of the m-th edge from node i to node j respectively. The above generative process produces observations between nodal pairs in the form {{(i, j, t(m ′) ij ), b (m′) ij , z (m′) ij } Mij m′=1}i ̸=j . For notational convenience, we reformat these as {(im, jm, tm), bm, zm}Mm=1, with m denotes the edge number among all the edges, for the following discussion.
In this generative process, steps 1 & 2 generate the exogenous function random bases and the nodes’ latent features, in the same way as for the NP-HP. The exogenous function in the NP-HPC ∑
k πikπjk · σ(fk(t)) is also nodal-pair specific. Step 3 generates the scaling variables for K endogenous functions. Step 4.(a) generates two directed edges between each nodal pair. Starting from any time stamp t∗, the Poisson processes with the current intensity functions (Eqs (6) and (7)) are used to generate the candidate edge occurrence time stamps t̄ij , t̄ji for each edge direction. If both t̄ij , t̄ji are larger than T , which means that the two Poisson processes do not generate edges in (t∗, T ], we stop generating edges for this nodal pair. Otherwise, we use the earlier candidate time stamp (e.g. we take t̄ij if t̄ij < t̄ji) as the actually occurring edge time, and add its triggering effect to the oppositely-directed intensity function (from the actual occurrence time).
Step 4.(b) uses the idea of branching structure [17, 32] to describe the “source” factor of the edge triplet (im, jm, tm). That is, bm = 0 if the triplet (im, jm, tm) is instantiated by an exogenous intensity and bm = m′(m′ ∈ Himjm(tm)) if the triplet (im, jm, tm) is triggered by its historical oppositely-directed edge triplet (jm, im, tm′). At the same time, zm ∈ {1, . . . ,K} is the latent label of (im, jm, tm) which determines the pattern group to which (im, jm, tm)’s triggering effect belongs. The edges with the same label use the same endogenous function. In this way, the endogenous functions are shared among the edges of all nodal pairs. Those edges triggered by the same labelled edges (through the same labelled endogenous function) or instantiated by the same labelled exogenous component are associated with the same endogenous function. This process is visualised in Figure 1.
As previously mentioned, NP-HP-C’s clustering mechanism for the endogenous functions of edges can help modelling practical phenomena, such as urgent emails would usually result in fast response, regardless of the nodal pair. What is more, NP-HP-C uses a much smaller number (K) of endogenous functions than NP-HP. This small number might be able to help NP-HP-C avoid overfitting issues.

heading: 3.4 Mean-field variational inference
text: To implement faster posterior inference than a MCMC implementation, we use the strategy of [5, 31] to develop mean-field variational posterior approximation [2]. In particular, [5] integrated a PolyGamma trick into a variational approximation for the sigmoidal Gaussian process-modulated Poisson process, and [31] applied similar techniques in the simple Hawkes point process setting. For e.g. the NP-HP-C model, we aim to approximate the joint posterior distribution of sigmoidal Gaussian process modulated random functions {fk(t), gk(t)}k, the endogenous function scaling values {vk}k, the latent features {πik}i,k, and the branching variables {bm}m and latent labels {zm}m for all edges. Full details of the variational model and updates are provided in the Supplementary Material.
Variational distributions for bm, zm: Since bm, zm are discrete random variables, we can adopt categorical distributions as their variational distributions. It is noted that when calculating the likelihood for zm, which denotes the latent label of the m-th edge (im, jm, tm), we need to include all the edges triggered by the m-th edge (i.e. (jm′ , im′ , tm′) : bm′ = m ∩ (im, jm) = (jm′ , im′)). In order to reduce the computational cost of the Gaussian process regression, we use the sparse Gaussian process approximation [27] to calculate Eq[fk(t)],Eq[f2k (t)],Eq[gk(t)] and Eq[g2k(t)], making it scale to the number of edges [5].
Computational complexity: The CTEM models provide an efficient solution to large-scale continuous-time edge data modelling. Let NE denote the number of observed edges, NET =∑
m #{Fimjm(tm)} denote the number of potential triggering phenomena between all oppositelydirected edges, NIN denote the number of inducing points for sparse Gaussian process approximation and K denote the number of latent features. The computational cost of the CTEM models then scales with O((NET +NE)K +NEN2INK), which scales linearly to the number of edges and would be suitable for large-scale continuous-time edge datasets.

heading: 4 Related work
text: Of the various approaches which use the ME-HP to model continuous-time edge data, the Hawkes Infinite Relational Model (Hawkes-IRM) [3] was the first to apply the ME-HP on the edges between groups formed by the nodes partition. By using the stochastic block model [22, 21, 14], the Hawkes-IRM model adopted static exogenous rates. The Hawkes-Dual Latent Space Embedding model (Hawkes-DLS) [29] created latent space vectors for nodes and used the products over these vectors to construct the exogenous and endogenous functions in the ME-HP. [26] adopted the nestedCRP [11, 1] to construct a hierarchical structure in the exogenous function when using the ME-HP
for continuous-time edges. [20] applied the same ME-HP in the framework of sparse exchangeable random measures [4, 28], and obtained continuous-time sparse exchangeable random measures. [30] adopted the Hierarchical Gamma Process to construct static exogenous function rates.
We note that the marginal probability of zm, obtained by summing over the values of bm in their joint probability function P (zm, bm) in NP-HP-C, is similar to the one used in the Dirichlet-Hawkes process (DHP) and its variants [7, 19, 24]. The NP-HP-C model has three major differences with these DHP models: 1) These DHP models consider the particular latent feature in triggering the current event. In contrast, the NP-HP-C model simultaneously specifies the particular latent feature (through zm) and the particular historical oppositely-directed edges (through bm); 2) We are able to develop a fast variational inference algorithm for this model (Section 3.4), whereas inference for DHP models relies on Sequential Monte Carlo methods; 3) The NP-HP-C model describes continuous-time edge data, whereas DHP models are mainly used for streaming text data.
Instead of using point processes to model the continuous-time edges, the approaches of [9, 10] studied time-discretised networks and used Bernoulli emission distributions to model all the binary-valued edges at each observed timestamp. They adopted Gaussian processes/Nested Gaussian processes to generate continuous-time features for each node at the observed timestamps. The approach of [18] also considered continuous-time edges. It extended the Poisson IRM by using histogram or kernel approximations to model the time-varying exogenous rate function, whereas our NP-PP uses Sigmoidal Gaussian Processes approximation.

heading: 5 Experiments
text: We perform experiments on four continuous-time edge datasets (Table 1) collected from the Stanford Large Social Network Dataset [16]. In each dataset the occurrence time of an edge is recorded to the nearest second. We scale them onto the interval [0, T ] by dividing the recorded time by 24× 3600, so that T is the total number of observation days. The first two datasets (Email, College) comprise instant response-typed edges (replying to emails and sending messages), while the last two (Overflow, Ubuntu) comprise edges recording relatively slow replies (writing answers to questions or commenting on others’ answers).
Experimental setup: We use the Gaussian kernel functions, κ(f)k (ti, tj) = θ (f) k e −0.5(ti−tj)2/(δ(f)k ) 2 and κ(g)k (ti, tj) = θ (g) k e −0.5(ti−tj)2/(δ(g)k ) 2
, in the Gaussian processes for generating random functions {fk(t)}Kk=1 and {gk(t)}Kk=1 respectively, where θ (f) k , θ (g) k are scaling parameters and δ (f) k , δ (g) k are bandwidth parameters. These parameters are optimized by applying the ADAM algorithm [15] on the Evidence Lower Bound of the variational distributions. The parameters av, bv of the endogenous function scaling values vk and the parameters aπ, bπ of the latent feature πik are given a prior distribution of Gam(0.1, 0.1) and optimized through the variational inference. We usually set 400 iterations for the MF-VI algorithm. Without particular specifications, we set the number of latent features as K = 5.
Comparison methods: The CTEM models are compared to the following benchmark methods: (1) the Poisson-Infinite Relational Model (Poisson-IRM) [3]; (2) the Hawkes-Infinite Relational Model (Hawkes-IRM) [3]; (3) the Hawkes Dual Latent Space model (Hawkes-DLS) [29]; and (4) Hawkes process with overlapping communities (HawkesNetOC) [20]. We use the exponential function as the endogenous functions and implement Markov chain Monte Carlo (MCMC) methods for the Poisson-IRM and the Hawkes-IRM models to the best of our ability. Inference for the Hawkes-DLS model and Hawkes-NetOC is executed by using the authors’ provided code. The default settings of Hawkes-NetOC and and Hawkes-DLS are used in the comparison. We use the previously developed CTM model constructions, giving the non-parametric Poisson process (NP-PP),
non-parametric Hawkes process (NP-HP), and non-parametric Hawkes process with edge clustering (NP-HP-C).

heading: 5.1 Testing Log-likelihood and link-prediction performance
text: The performance of each model is evaluated through its model fitting performance to the Testing edges and its capability to predict future edges. The edges of each dataset are sorted according to their occurrence times, and the dataset split into a training set (the first 70% of edges) and a testing set (the remaining 30%). The comparisons between the CTEM and other models are based on two criteria: the log-likelihood for the edges in the testing time period, and the AUC values on randomly sampled sub-interval in the testing time period. We use Eq. (2) to calculate the log-likelihood, in which the intensity function λij(t) is either the exogenous function (for Poisson processes) or a combination of the exogenous and endogenous functions (for Hawkes processes).
For the AUC calculation, we first uniformly sample 100 time stamps ti in the time interval of the test dataset, and then evaluate the probability of an edge appearing between all nodal pairs in the time window [ti, ti + δ). We compute AUC values for each time window [ti, ti + δ). The final link-prediction performance score is the mean of these 100 AUC values on all time windows.
The results are displayed in Table 2. Due to the space limit, we have provided the standard error in the Supplementary Material. The CTEM models perform differently for the two types of datasets. For the instant-response datasets, the NP-HP and NP-HP-C produces better performance in both testing log-likelihood and AUC. This may indicate that our exogenous functions and endogenous functions are fitted better than the comparison methods. However, for the slow-response datasets, the CTEM models performs worse in the testing log-likelihoods, although it still produces better prediction performance. A possible reason is that these two datasets have relatively simple patterns (here low-response activities) which can be well captured by simple rate functions. For instance, the exogenous functions for these slow-response datasets did not show clear periodicities (see the right two panels in the top row of Figure 2). Simple rate functions like constant functions may well describe the exogenous effects. For the endogenous effects, they are dominated by one single endogenous function in both Overflow and Ubuntu (see the right two panels in the bottom row of Figure 2). One single exponential decay function may also capture these simple patterns.

heading: 5.2 Visualisations for exogenous functions and endogenous functions
text: The top row in Figure 2 presents the mean exogenous function values ∑ (i,j):i ̸=j ∑
k πikπjk ·σ(fk(t)) for the NP-PP, NP-HP and NP-HP-C models. With the vertical dotted lines separating the time lines into weeks, daily periodicity is quite clear for the instant-response datasets (i.e. Email, College), and approximate weekly periodicity for Overflow and Ubuntu. The weekend effect is also easily identifiable in Email, a gradual increasing trend in location and scale for College, and occasional weeks of peak high and low activity for the slow response datasets. For the latter, the large activity spikes could be explained by sudden interest in popular topics being discussed on the message boards.
In terms of performance difference within the CTEM models, the exogenous values of the NP-PP model (blue lines) are always higher than those of the NP-HP and NP-HP-C models, and the NP-HP (green lines) usually produces the smallest values. The former may occur as the NP-PP model is using the exogenous functions to explain all edges as it does not have endogenous functions. The latter may occur as the NP-HP models less edges than the NP-HP-C. For the slow-response datasets,
the average values are similar for the NP-PP and NP-HP-C models. This may be explained as the endogenous effects may be more important to the instant-response data than the slow-response data.
The bottom row in Figure 2 illustrates the most influential endogenous functions of the NP-HP-C model. Within our expectation that endogenous effects would become weaker along with the time, the values of all main endogenous functions have sharp decreases in the beginning period. Three of the datasets also exhibit (daily or weekly) periodicity, which indicates that for these data there has been an increase in triggered responses above and beyond the increase in events attributable to background (exogenous) rates. This is less obviously the case for the College dataset, where, despite some initial periodicity, many of the endogenous functions become constant after a certain time. We are unaware of any related discussions of periodic endogenous functions for these data in the literature.
We have provided visualisations for other variables in the Supplementary Material.

heading: 5.3 Intensity functions λij(t) and influence of the number of latent features
text: The left two panels of Figure 3 show the intensity functions λij(t) of the NP-HP-C for randomly selected nodal pairs on an instant-response (Email) and a slow-response (Overflow) dataset. The daily periodicity is again strongly apparent in λij(t) for Email, and more muted weekly periodicity for λij(t) in Overflow. We note that most of the observed edges (points on the x-axis) are located in regions where the intensity increases, indicating correct performance.
The right two panels of Figure 3 show the influence of different number of latent features K. Predictive performance appears less sensitive to K for the instant-response datasets, which perhaps have simpler structures, than the slow-response datasets, as the latter have improved performance with increasing K. Consistent with the complexity analysis, running times are linear in the number of edges. As the NP-HP-C with K = 5 requires ∼ 60 seconds for one iteration with 80 000 edges, we need ∼ 6 hours of for posterior inference on these edges (with 400 iterations).

heading: 6 Conclusion
text: In this work, we have explored three non-parametric point process models for continuous-time edge data. The sigmoidal Gaussian process priors for the exogenous and endogenous functions enable flexible modelling, and strong predictive performance. The developed mean-field variational
inference algorithm allows the models to work at scale. The experimental results on four real world data sets demonstrate superior performance compared to the state-of-the-art.

heading: Acknowledgements
text: Xuhui Fan and Scott A. Sisson are supported by the Australian Research Council through the Australian Centre of Excellence in Mathematical and Statistical Frontiers (ACEMS, CE140100049), and Scott A. Sisson through the Future Fellow Scheme (FT170100079). Bin Li is supported by the Program for Professor of Special Appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning. Feng Zhou is partially funded by China Postdoctoral Science Foundation.

Review summaries:
The paper proposes three  non-parametric point process variants for continuous-time edge modelling (CTEM) and the performance of the proposed framework evaluated on four real world dynamic networks. The three cases are  1) a Poisson process in which K time series are generated according to a (sigmoidal) GP and weighted by node specific feature weights utilizing these temporal evolutions. 2) An extension to the Hawkes process such that reciprocity is invoked by an additional sigmoidal GP weighted by other node specific feature weights defined through time to past reciprocative interactive events. 3) A further extension incorporating edge clustering rather than node specific weights for the reciprocating effects (i.e., endogeneous functions) reducing the parameter space and enabling specific endogenous effects according to the cluster used for the specific reciprocative event of each dyad.   VB is used for large scale inference and inducing points invoked to reduce complexity of the GP. The model is contrasted the Poisson and Hawkes IRM as well as Hawkes-DLS and HawkesNetOC and found on the four considered datasets to have superior performance as quantified using time-binned AUC but inferior results on two datasets for the (test I assume) log-likelihood. In particular, the edge clustering approach 3) forming the NP-HP-C is found to be the superior of the proposed approaches. 

The authors propose a family of models based on mutually-exciting Hawkes process that aim at modeling continuous time edges. The framework is based on a variational inference procedure for a fast and efficient computation. The models have been quantitatively and qualitatively evaluated with respect to baselines, showing that, in general, they outperform the baselines in terms of accuracy.

-This work develops a family of novel non-parametric point processes for continuous-time edge data. In details, the model captures the periodicity exhibited in its base rate using K sigmoidal Gaussian processes(SGPs) \sigmoid(f_k(t)).  -For each pair of two nodes, the nodal pair specific coefficients \p_{ik}, \pi_{jk} is used to rescale the pair-wise base intensity of (i,j) and (j,i).   -To  capture the reciprocity, the model employs a second set of functions \sigmoid(g_k(t)) drawn from SGPs, which has more flexibility in capturing exciting behaviors, compared to the exponential decaying kernel used in Hawkes processes. - A mean-field variational scheme is derived to perform model inference. -The model is demonstrated on four real data to illustrate its estimated periodically varying base intensity and the endogenous functions. They also show the superior performance of the model in terms of training log-likelihood and link prediction.