# Deep graph kernel point processes

Zheng Dong1, Matthew Repasky1, Xiuyuan Cheng2, and Yao Xie 1 Email:
<EMAIL_ADDRESS>
(1 H. Milton Stewart School of Industrial & Systems Engineering, Georgia
Institute of Technology
2 Department of Mathematics, Duke University )

###### Abstract

Point process models are widely used for continuous asynchronous event data,
where each data point includes time and additional information called “marks”,
which can be locations, nodes, or event types. This paper presents a novel
point process model for discrete event data over graphs, where the event
interaction occurs within a latent graph structure. Our model builds upon
Hawkes’s classic influence kernel-based formulation in the original self-
exciting point processes work to capture the influence of historical events on
future events’ occurrence. The key idea is to represent the influence kernel
by Graph Neural Networks (GNN) to capture the underlying graph structure while
harvesting the strong representation power of GNNs. Compared with prior works
focusing on directly modeling the conditional intensity function using neural
networks, our kernel presentation herds the repeated event influence patterns
more effectively by combining statistical and deep models, achieving better
model estimation/learning efficiency and superior predictive performance. Our
work significantly extends the existing deep spatio-temporal kernel for point
process data, which is inapplicable to our setting due to the fundamental
difference in the nature of the observation space being Euclidean rather than
a graph. We present comprehensive experiments on synthetic and real-world data
to show the superior performance of the proposed approach against the state-
of-the-art in predicting future events and uncovering the relational structure
among data.

Keywords: Point processes over graphs, Deep kernel, Graph neural networks

##  1 Introduction

Asynchronous discrete event data, where each data point includes time and
additional information called “marks”, are ubiquitous in modern applications
such as crime [59], health care [42], earthquake events [31, 57], traffic
modeling [56], and so on. In contemporary applications, the collection of
discrete events often reveals an underlying latent graph structure, leading to
the widespread adoption of models incorporating graph structures for various
purposes. Point processes over latent graphs are a popular model for such
data, where the graph nodes can be introduced to capture event marks, such as
locations or event types.

Classic one- and multi-dimensional temporal, self-exciting point process
models introduced by Hawkes [22] leverage an event influence kernel function
to capture the impact of historical events on future events’ occurrence. The
influence kernel takes an exponentially decaying form, often not expressive
enough to capture complex influence mechanisms. Recently, there have been many
successes in deep point process models that represent the influence kernel
using neural networks for temporal-only kernels [60], spatio-temporal kernels
[32], and non-stationary kernels [12, 58]. Such works achieve competitive
performance through efficient modeling of the influence kernel compared to
point processes that model the conditional intensity function using neural
networks. The key lies in that modeling event influence through a kernel
captures the repeated influence patterns more effectively by combining
statistical and deep models, thus achieving better model estimation/learning
efficiency and superior predictive performance.

Despite much success in prior work on deep kernel modeling of point processes,
there has been limited work in exploiting underlying graph structures of
multi-dimensional point processes when characterizing the influence kernel.
Meanwhile, recent developments in Graph Neural Networks (GNNs) offer potential
solutions to modeling deep graph kernels with point process data by harvesting
their representation power. However, how to properly adopt GNNs into point
processes while preserving the statistical model interpretability remains an
open question.

This paper presents a novel point process model with a deep graph influence
kernel for discrete event data over graphs. The graph kernel is represented by
a properly encoded GNN architecture to capture the underlying dynamics of
event influence. Moreover, the kernel architecture proposes a unified
framework for integrating various GNN structures with point processes via
localized graph filter basis functions. We also present a computationally
efficient and flexible learning scheme by considering two types of loss
functions: the commonly used maximum likelihood estimation (MLE) and a new
least-square estimation (LSE) scheme. We can allow general types of non-
stationary, multi-hop influence (which can be negative), and a log-barrier
penalty in the loss function ensures the non-negative constraint for the
conditional intensity function.

The modeling of the deep graph influence kernel also enables the
characterization of dynamic graph influence through asynchronous discrete
point process observations. The proposed kernel can capture the influence
among discrete events originating from various nodes at different times with
the encoded GNN architectures. The temporal relationships of these events can
further represent the time-evolving interactions among graph nodes. Meanwhile,
leveraging the power of Graph Attention Networks (GATs) [41], our kernel can
recover the underlying graph structure through the learned attention weights
between nodes.

Our work significantly extends the existing deep spatio-temporal kernel for
point process data, which cannot be applied in our setting since the
observation space is fundamentally different: the spatio-temporal point
processes are for events that occurred in physical Euclidean space, and point
processes over graphs represent vicinity using nodes and edges. Our
contributions can be summarized as follows:

  1. 1.

Our proposed method explicitly models the influence kernel in point processes
via GNNs instead of typical intensity-based models. This permits greater
expressivity of inter-event-category contributions, including non-stationary,
multi-hop exciting, and inhibiting effects. Furthermore, the deep graph kernel
can capture dynamic graph influence, yielding clear information about the
relational structure in graph point processes.

  2. 2.

The localized graph filter basis allows the deep kernel to go beyond simple
distance-based influence for graphs representing events in space, providing a
model structure for non-spatial graphs such as traffic or social networks.
Meanwhile, a larger class of GNN models can be incorporated within our
framework, enabling broader versatility in real-world applications.

  3. 3.

Comprehensive experiments demonstrate the benefits of including latent graph
structures in the deep kernel modeling over the state-of-the-art in simulated
and real data settings. Our method applies to a wide array of point process
data settings, including events generated by infrastructural, climatic, and
social phenomena.

###  1.1 Related Works

, There are various deep learning point process models based on modeling the
conditional intensity function using neural networks (rather than the
influence kernel), such as recurrent neural networks (RNNs) [14, 28]. Due to
advances in attention models for sequential data modeling [40], RNN approaches
have been surpassed by self-attention-based approaches [50, 62]. These RNN and
self-attention methods provide expressive models for the conditional
intensity; however, they often lack statistical interpretability coming from
the original “influence kernel” approach, and they do not consider graph
structures.

Modeling point processes with graph information has been explored in early
attempts. Classical multivariate Hawkes processes [36] assume a parametric
form of the conditional intensity. A recent work [16] develops the novel Group
Network Hawkes Process that can account for heterogeneous nodal
characteristics using a parametric network model. Our approach differs from
parametric point processes by assuming a deep graph-based influence kernel.
Many modern approaches [34, 48, 51, 62] adopt non-graph neural networks
incorporated with the graph structure, rather than directly incorporating GNN
structures in the models. For example, A-NHP [48] and THP-S [62] adopt
attention-based mechanisms with additional constraints on the learned
attention weights posted by the graph structure, which differs from the Graph
Attention Network (GAT) [41]. These approaches mainly consider the single-hop,
adjacency-based influence over the latent graph. The Geometric Hawkes Process
(GHP) [38] integrates graph convolution RNNs with the Hawkes process to
estimate the parameters of the conditional intensity function with a
parametric (exponentially decaying) kernel, achieving enhanced model
expressiveness. In our work, incorporating GNN architectures in the influence
kernel without any parametric constraints permits the flexible and
interpretable recovery of complex event dependency, such as multi-hop
influence mechanisms over the graph topology.

Our work is related to GNNs, which have seen wide applications in areas
including temporal phenomena on graphs [27, 44]. Graph convolutions
incorporate spatial or spectral convolutions based upon the adjacency matrix
or graph Laplacian [5]. An early attempt applying spatial convolutions is the
Diffusion-Convolutional Neural Network (DCNN) [1]. Other spatial-based
approaches include attention models [4, 15, 41]. Prototypical spectral
convolutions are utilized in Chebnet [11] and graph convolutional networks
(GCN) [26], with recent extensions including auto-regressive moving average
(ARMA) spectral convolutions [3] and the Simple Spectral Graph Convolution
(SSGC) [55]. Modern approaches incorporate both local and global features,
such as L3Net [9] and General, Powerful, Scalable (GPS) Graph Transformer
[35].

Prior work has been done aiming to infer graph influence or structures given
observed data, such as spectral and graph signal processing methods via
optimization of the graph Laplacian under certain assumptions [45]. Deep
learning methods and GNNs have also been employed in graph structure learning
[8, 18, 61]. The components of the adjacency matrix can be treated as
learnable parameters, subject to downstream tasks [17, 23]. The utilization of
point process data in graph structure/influence learning has been explored
[47, 52, 53]. The interaction between two nodes is represented by a learnable
scalar, which may not be expressive enough to capture the complexity of the
interactions, especially in modern applications. Our model can learn
complicated and dynamic nodal interactions by adaptively learning the graph
influence kernel from point process observations, facilitated by a carefully
crafted framework that enables the incorporation of various GNN architectures.

##  2 Background

##### Self-exciting point process.

A self-exciting point process [36] models the occurrence of time-stamped
discrete events that depend on the observed history. Consider a simple
temporal point process that only models event times. Let
$\mathcal{H}=\\{t_{1},\dots,t_{n}\\}$ be an observed event sequence, where
$t_{i}\in[0,T]\subset\mathbb{R}$ is the time of the $i$-th event. We denote
the history before a given time $t$ as $\mathcal{H}_{t}=\\{t_{i}|t_{i}<t\\}$.
The conditional intensity of events is defined as $\lambda(t)=\lim_{\Delta
t\downarrow 0}\mathbb{E}\left[\mathbb{N}([t,t+\Delta
t])|\mathcal{H}_{t}\right]/\Delta t,$ where the counting measure $\mathbb{N}$
is defined as the number of events occurring in $[t,t+\Delta t]$. For
notational simplicity, we omit the dependency of history $\mathcal{H}_{t}$ in
$\lambda(t)$. The well-known Hawkes process [22] models the self-excitation
effect from history in an additive manner. The conditional intensity function
is defined as

$\lambda(t)=\mu+\sum_{t^{\prime}\in\mathcal{H}_{t}}k(t^{\prime},t),$

where $\mu$ is the background intensity and $k$ is the so-called influence
kernel measuring the effects of historical events.

In a marked point process, each event is associated with an additional
attribute called mark denoted by $v\in V$. The mark represents specific
characteristics of the event and can be either continuous or categorical, such
as event location or event type. Let
$\mathcal{H}=\\{(t_{i},v_{i})\\}_{i=1}^{n}$ and
$\mathcal{H}_{t}=\\{(t_{i},v_{i})|t_{i}<t\\}$ be the observed event sequence
and history before time $t$, respectively. The conditional intensity with
influence kernel $k$ can be written as:

$\lambda(t,v)=\mu+\sum_{(t^{\prime},v^{\prime})\in\mathcal{H}_{t}}k(t^{\prime},t,v^{\prime},v).$
(1)

The influence kernel is crucial when learning the conditional intensity
$\lambda(t,v)$ from event sequences. A standard and simplified way [36] is to
represent the kernel $k(t^{\prime},t,v^{\prime},v)$ to be a product of spatial
interaction with temporal kernel in the form of
$a_{v,v^{\prime}}f(t-t^{\prime})$, where the coefficient $a_{v,v^{\prime}}$
captures the influence of node $v^{\prime}$ on $v$ through a graph kernel, and
$f$ is a stationary temporal kernel. We consider a general graph kernel
representation using GNNs that goes beyond the parametric form in the classic
point process literature, thus enabling better characterization of event
dynamics (Appendix C).

##### Graph convolution.

Graph convolutions in GNNs [44] extend the convolution strategy to graphs. It
falls into two categories: spectral- and spatial-based models. Spectral graph
convolutions introduce graph filters $g_{\theta}$ based on the full eigen-
decomposition of the graph Laplacian. The graph signal $X$ is convoluted by
$X\ast_{G}g_{\theta}=Ug_{\theta}U^{T}X$, where $U$ is the matrix of the
eigenvectors of the graph Laplacian ordered by eigenvalues. For instance, in
Spectral Convolutional GNNs [5], the graph filter $g_{\theta}=\Theta_{i,j}$
contains a set of learnable parameters that characterize the relations between
node pairs. On the other hand, spatial-based graph convolution is performed by
information propagation along edges. The weight matrix in each layer is
constructed based on the node’s spatial relations (i.e., adjacency matrix).
Either the localized filter or the weight matrix plays a pivotal role in
capturing the nodal dependencies.

##  3 Point processes on graphs

###  3.1 Problem definition

We aim to develop a point process model for the occurrence of multiple types
of events within a latent graph structure. Let $G=(V,E)$ denote the underlying
graph, where each node $v\in V$ represents one event type. An undirected edge
connecting nodes $u$ and $v$ indicates the existence of potential interaction
between type-$u$ and type-$v$ events. Note that the edges merely suggest the
support of possible inter-event-category interactions without dictating the
directions.

Consider a set of event sequences
$\mathcal{S}=\\{\mathcal{H}^{1},\mathcal{H}^{2},\dots,\mathcal{H}^{|\mathcal{S}|}\\}$,
where each $\mathcal{H}^{s}=\\{(t_{i}^{s},v_{i}^{s})\\}_{i=1}^{n_{s}}$ is a
collection of events $(t_{i}^{s},v_{i}^{s})$ occurring on node $v_{i}^{s}$ at
time $t_{i}^{s}$. Our proposed graph point process is expected to (i) jointly
predict the times and types of forthcoming events based on the observed
historical data and (ii) provide an interpretable understanding of the event
generation process by revealing the interdependences among multiple types of
events and uncovering the latent graph structure with no prior information.
Toward this end, we adopt the statistical formulation of conditional intensity
in (1) and introduce an influence kernel built on localized graph filters in
GNNs, aiming to explicitly characterize the complicated contributing
relationship between any binary event pair (e.g., excitation, inhibition, or
other dynamic influences).

###  3.2 Deep temporal graph kernel

Modeling the multi-dimensional influence kernel $k$ for intricate event
dependency is crucial yet challenging. To go beyond simple parametric forms of
the kernel while maintaining the model efficiency, we represent the multi-
dimensional kernel by taking advantage of the kernel singular value
decomposition (SVD) [29, 30]. Specifically, the influence kernel
$k(t^{\prime},t,v^{\prime},v)$ in (2) is decomposed into basis kernel
functions as follows:

$k(t^{\prime},t,v^{\prime},v)=\sum_{d=1}^{D}\sigma_{d}g_{d}(t^{\prime},t-t^{\prime})h_{d}(v^{\prime},v),$
(2)

where $\\{g_{d},h_{d}\\}_{d=1}^{D}$ are sets of basis kernels in terms of
event time and type, respectively. The scalar $\sigma_{d}$ is the
corresponding weight (or “singular value”) at each rank $d$. Instead of
directly learning the multi-dimensional event dependency, we simplify the task
by “separately” modeling specific modes of event dependency over time or graph
using different basis kernels. It is worth noting that the weighted
combination of basis kernels covers a broad range of non-stationary influence
kernels used in point processes, and our kernel $k$ is not decoupled over time
and graph space. While functional SVD is usually infinite-dimensional, in
practice, we can truncate the decomposition as long as the singular values
$\sigma_{k}$ decay sufficiently fast, only considering a finite rank
representation (see an example in Appendix G).

The temporal basis kernels are carefully designed to capture the heterogeneous
temporal dependencies between past and future events. First, the
parametrization of temporal kernels $\\{g_{d}\\}_{d=1}^{D}$ using
displacements $t-t^{\prime}$ instead of $t$ provides us a low-rank way to
approximate general kernels [12]. To proceed, we approximate
$\\{g_{d}\\}_{d=1}^{D}$ using shared basis functions:

$g_{d}(t^{\prime},t-t^{\prime})=\sum_{l=1}^{L}\beta_{dl}\psi_{l}(t^{\prime})\varphi_{l}(t-t^{\prime}),\quad\forall
d=1,\dots,D.$

Here $\\{\psi_{l},\varphi_{l}:[0,T]\to\mathbb{R}\\}_{l=1}^{L}$ are two sets of
one-dimensional basis functions characterizing the temporal impact of an event
occurring at $t^{\prime}$ and the pattern of that impact spread over
$t-t^{\prime}$. The scalar $\beta_{d,l}$ is the corresponding weight. The
basis functions $\\{\psi_{l},\varphi_{l}\\}_{l=1}^{L}$ are represented by
independent fully connected neural networks, leveraging their universal
approximation power.

###  3.3 Graph basis kernel with localized graph filters

We develop a novel framework for the graph basis kernels by leveraging the
localized graph filters in graph convolution to extract informative inter-
event-category patterns from graph-structured data. Specifically, the basis
kernels $\\{h_{d}\\}_{d=1}^{D}$ are represented as follows:

$h_{d}(v^{\prime},v)=\sum_{r=1}^{R}\gamma_{dr}B_{r}(v^{\prime},v),\quad\forall
d=1,\dots,D,$

where $\\{B_{r}(v^{\prime},v):V\times V\to\mathbb{R}\\}_{r=1}^{R}$ are $R$
bases of localized graph filters, and $\gamma_{dr}$ is the corresponding
weight for each $B_{r}$. The bases can be constructed from a spatial or a
spectral approach, corresponding to two categories of commonly seen graph
convolutions.

Formally, the influence kernel $k$ can be represented as:

$k(t^{\prime},t,v^{\prime},v)=\sum_{r=1}^{R}\sum_{l=1}^{L}\alpha_{rl}\psi_{l}(t^{\prime})\varphi_{l}(t-t^{\prime})B_{r}(v^{\prime},v),$
(3)

where
$\alpha_{rl}=\sum_{d=1}^{D}\sum_{r=1}^{R}\sum_{l=1}^{L}\sigma_{d}\beta_{dl}\gamma_{dr}$.
To showcase the model flexibility of our graph basis kernel to incorporate
various GNN structures, we implement two examples of L3Net [9] and GAT [41] in
our numerical experiments (Section 4). L3Net provides a unified framework for
both spatial- and spectral-based graph convolutions, and GAT can predict the
presence or absence of an edge in a graph. Note that our framework is
compatible with general GNN layer types, such as those in Chebnet [11] and GPS
Graph Transformer [15]. An example of L3net incorporation is given in Appendix
B. Technical details of other GNN incorporations are presented in Appendix C.

The benefits of our design for the influence kernel $k$ are as follows. (i)
The kernel enables the adoption of various spectral and spatial filter bases.
The combination of $R$ bases allows us to represent complex local and global
patterns of inter-node influences with great model expressiveness. (ii) Our
framework substantially reduces the number of model parameters to
$\mathcal{O}(RC|V|)$ for modeling graph-structured point process data with
$|V|$ event types, while classic multivariate point processes and other neural
point processes typically require more than $\mathcal{O}(|V|^{2})$ parameters.
Here $C$ represents the average local patch size [9]. In practice, we have
$C,R\ll|V|$ when dealing with sparse graphs and considering only up to $o$-hop
influence (commonly 2 or 3), which significantly improves the scalability of
our model when applied to large graphs. Details of the complexity analysis can
be found in Appendix C.

Choice of the network hyperparameters. In practice, we can treat $\alpha_{rl}$
as one learnable parameter and the model parameter $D$ is absorbed. We also
provide two strategies for determining the kernel rank $L$ and $R$ in Appendix
G.

###  3.4 Model estimation

Previous studies of point processes primarily use two types of loss function
for model estimation, including the (negative) log-likelihood function (NLL)
[13, 36] and the least square loss function (LS) [2, 6]. We present the
incorporation of these two losses in our framework, which is not specific to a
particular choice of loss and, in fact, can be quite general. We present
derivations of two loss functions in Appendix D for completeness.

##### Negative Log-Likelihood (NLL).

The model parameters can be estimated by minimizing the negative log-
likelihood of observing event sequences $\mathcal{S}$ on $[0,T]\times V$ [36],
which is expressed as:

$\ell_{\text{NLL}}(\theta):=\frac{1}{|\mathcal{S}|}\sum_{s=1}^{|\mathcal{S}|}\left(\sum_{v\in
V}\int_{0}^{T}\lambda(t,v)dt-\sum_{i=1}^{n_{s}}\log\lambda\left(t_{i}^{s},v_{i}^{s}\right)\right).$

Note that the model parameter $\theta$ is incorporated into the intensity
function $\lambda$. This approach is equivalent to the maximum likelihood
estimation approach, and the log-likelihood indicates how well the model fits
the occurrences of events.

##### Least Square (LS) loss.

Another approach to the point process estimation is based on least square loss
[21] of observed sequences $\mathcal{S}$:

$\ell_{\text{LS}}(\theta):=\frac{1}{|\mathcal{S}|}\sum_{s=1}^{|\mathcal{S}|}\left(\sum_{v\in
V}\int_{0}^{T}\lambda^{2}(t,v)dt-\sum_{i=1}^{n_{s}}2\lambda(t_{i}^{s},v_{i}^{s})\right).$

The least-square loss function can be derived from the empirical risk
minimization principle [19]. Intuitively, we expect the integral of intensity
over infinitesimal time intervals containing event times to be approximately
one and approximately zero during non-event times.

Since negative values of the influence kernel are allowed for learning
inhibiting effects, an additional constraint is required for ensuring the non-
negativity of the conditional intensity function. To achieve this, we leverage
the log-barrier method for optimization in point processes [12], which
maintains the model interpretability while being computationally efficient.
The final objective function to be minimized is formulated as

$\displaystyle\min_{\theta}\mathcal{L}_{\text{NLL}}(\theta):=\ell_{\text{NLL}}(\theta)+\frac{1}{w}p(\theta,b)$
and $\displaystyle\min_{\theta}\mathcal{L}_{\text{LS}}(\theta)\ \
:=\ell_{\text{LS}}(\theta)+\frac{1}{w}p(\theta,b)$

for two loss functions, respectively. Here
$p(\theta,b)=-\frac{1}{|\mathcal{S}||\mathcal{U}_{\text{bar},t}\times
V|}\sum_{s=1}^{|\mathcal{S}|}\sum_{t\in\mathcal{U}_{\text{bar},t}}\sum_{v\in
V}\log(\lambda(t,v)-b)$ is the log-barrier function with regularization
parameter $w>0$. The $\mathcal{U}_{\text{bar},t}$ is a dense enough grid over
$[0,T]$ that the intensity function is penalized on, and $b$ is a lower bound
of the intensity value, which will be tuned to 0, to ensure the initialization
is within the feasible domain of the logarithm function. Both loss functions
can be efficiently computed numerically, as illustrated in Appendix E.

##  4 Experiment

In this section, we compare our method with a deep graph kernel, referred to
as GraDK, with seven state-of-the-art point process methods on both synthetic
and real data sets.

##### Baselines.

We include two groups of baselines with distinctive model characteristics.
Models in the first group treat the associated node information of each event
as one-dimensional event marks without considering the graph structure,
including (i) Recurrent Marked Temporal Point Process (RMTPP) [14] that uses a
recurrent neural network to encode dependence through time; (ii) Fully Neural
Network model (FullyNN) [33] that models the cumulative distribution via a
neural network; and (iii) Deep Non-Stationary Kernel [12] with a low-rank
neural marked temporal kernel (DNSK-mtpp). The second group includes two
models that encode the latent graph structure information when modeling event
times and marks, including (iv) Structured Transformer Hawkes Process (THP-S)
[62] and (v) Graph Self-Attentive Hawkes Process (SAHP-G) [50] with a given
graph structure, which both use self-attention mechanisms to represent the
conditional intensity. Comparisons with more methods are presented in Appendix
A.

##### Experimental setup.

We demonstrate the adaptability of the proposed method to various advanced GNN
architectures and loss functions using three different architectures: (i)
GraDK with L3Net and NLL (GraDK+L3net+NLL), (ii) GraDK with L3Net and LS
(GraDK+L3net+LS), and (iii) GraDK with GAT and NLL (GraDK+GAT+NLL). Details
about the experimental setup and model architectures can be found in Appendix
H.

Table 1: Synthetic data results.

| 3-node graph with negative influence | 16-node graph with 2-hop influence | 50-node graph | 50-node graph with non-separable kernel | 225-node graph
---|---|---|---|---|---
Model | Testing $\ell$ | Time MAE | Type KLD | Testing $\ell$ | Time MAE | Type KLD | Testing $\ell$ | Time MAE | Type KLD | Testing $\ell$ | Time MAE | Type KLD | Testing $\ell$ | Time MAE | Type KLD
RMTPP | $-3.473_{(0.087)}$ | $0.528$ | $0.093$ | $-7.239_{(0.193)}$ | $0.301$ | $0.142$ | $-27.915_{(1.251)}$ | $37.666$ | $0.103$ | $-9.560_{(0.513)}$ | $41.606$ | $0.371$ | $-16.294_{(1.047)}$ | $29.329$ | $0.252$
FullyNN | $-2.086_{(0.009)}$ | $0.291$ | $0.006$ | $-3.347_{(0.018)}$ | $0.198$ | $0.018$ | $-1.736_{(0.019)}$ | $13.295$ | $0.058$ | $-3.098_{(0.176)}$ | $25.080$ | $0.054$ | $-3.241_{(0.026)}$ | $9.672$ | $0.174$
DNSK-mtpp | $-2.127_{(0.003)}$ | $0.149$ | $0.012$ | $-3.005_{(0.002)}$ | $0.085$ | $0.002$ | $-1.165_{(0.003)}$ | $1.074$ | $0.076$ | $-1.595_{(0.003)}$ | $1.043$ | $0.005$ | $-2.511_{(0.001)}$ | $3.958$ | $0.086$
THP-S | $-2.089_{(0.008)}$ | $0.413$ | $0.006$ | $-3.079_{(0.004)}$ | $0.108$ | $0.011$ | $-1.091_{(0.005)}$ | $3.940$ | $0.019$ | $-1.587_{(0.004)}$ | $3.110$ | $0.003$ | $-2.550_{(0.008)}$ | $4.109$ | $0.087$
SAHP-G | $-2.113_{(0.005)}$ | $0.172$ | $0.003$ | $-3.036_{(0.008)}$ | $0.155$ | $0.005$ | $-1.099_{(0.004)}$ | $1.119$ | $0.014$ | $-1.591_{(0.005)}$ | $3.482$ | $0.009$ | $-2.506_{(0.005)}$ | $3.578$ | $0.094$
GraDK+L3net+NLL | $-2.062_{(0.003)}$ | $0.048$ | <0.001 | $-2.995_{(0.004)}$ | $0.028$ | $0.002$ | $-1.065_{(0.005)}$ | $1.023$ | 0.004 | $-1.583_{(0.002)}$ | 0.244 | <0.001 | $-2.487_{(0.005)}$ | $1.779$ | $0.015$
GraDK+L3net+LS | $\textbf{--2.059}_{(0.002)}$ | 0.021 | <0.001 | $\textbf{--2.993}_{(0.003)}$ | 0.001 | 0.001 | $\textbf{--1.056}_{(0.003)}$ | 0.957 | $0.005$ | $\textbf{--1.582}_{(0.002)}$ | $0.253$ | <0.001 | $\textbf{--2.485}_{(0.003)}$ | 0.146 | 0.012
GraDK+GAT+NLL | $-2.073_{(0.004)}$ | $0.068$ | <0.001 | $-2.997_{(0.005)}$ | $0.148$ | <0.001 | $-1.690_{(0.001)}$ | $0.690$ | $0.007$ | $-1.584_{(0.003)}$ | $0.354$ | $0.003$ | $-2.493_{(0.006)}$ | $2.458$ | $0.092$

  * •

*Numbers in parentheses are standard errors for three independent runs.

###  4.1 Synthetic data

We evaluate the efficacy of our model on five synthetic data sets generated by
point processes with the following kernels and latent graph structures: (i) a
non-stationary temporal kernel on a 3-node graph with negative influence; (ii)
a non-stationary temporal kernel on a 16-node graph with 2-hop graph
influence; (iii) an exponentially decaying temporal kernel on a 50-node graph
with central nodes; (iv) a non-separable diffusion kernel on a 50-node graph
with central nodes; and (v) an exponentially decaying temporal kernel on a
225-node graph. Details regarding synthetic data are presented in Appendix H.

####  4.1.1 Known graph topology

We demonstrate the exceptional performance of GraDK in modeling point process
data by leveraging observed latent graph structures using synthetic data sets.
These experiments aim to simulate scenarios where prior knowledge of the
latent space is available.

##### Kernel and intensity recovery.

Figure 1 contains the recovered graph kernel by each method for the synthetic
data generated by the kernel on a 16-node ring graph with 2-hop influence. Our
method and DNSK directly learn the kernel, and the graph kernel in SAHP-G is
constructed as in the original paper [50] by computing the empirical mean of
the learned attention weights between nodes. It is worth noting that our model
learns an accurate representation, recovering the self-exciting and multi-hop
influential structures in the kernel over the latent graph. At the same time,
SAHP-G only captures the mutual dependencies within one-hop neighbors,
restricted by their model formulation. Our methods also recover the event
intensity successfully (the bottom row of Figure 1). The conditional
intensities of SAHP-G and DNSK, however, either fail to capture the magnitude
of this interdependence or do not accurately decay node dependence along the
ring-structure connections.

Figure 1: Graph kernel, inter-event dependence, and conditional intensity
recovery for the 16-node synthetic data set with 2-hop graph influence. The
first column reflects the ground truth, while the subsequent columns reflect
the results obtained by GraDK (with three architectures), SAHP-G, and DNSK,
respectively. Figure 2: Dynamic influence among nodes learned by GraDK on
synthetic data generated by the 50-node graph with non-separable diffusion
kernel. Both true and learned kernels are evaluated at three-time lags of
$0.3,0.8$, and $1.5$.

##### Event dependency.

Our model exhibits exceptional performance in capturing sequential event
dependencies. The second row of Figure 1 visualizes the learned inter-event
dependency given a testing sequence generated on the 16-node graph with 2-hop
influence. The dependency between a prior and a future event is characterized
by the influence kernel (2) in GraDK, DNSK, and the true model. For SAHP-G,
the event dependency is indicated by the scaled self-attention weight
(Equation 10 [50]). While SAHP-G is capable of discovering long-term
dependencies, the decaying influence of recent events is not represented. DNSK
captures the decaying influence but fails on the long-term effects of certain
event types. Our method learns both of these features resembling those of the
true model. Similar results on event dependencies on other synthetic data sets
are shown in the second row of Figure H5-H8, indicating the superior
performance of our model against SAHP-G and DNSK.

##### Predictive ability.

The superior predictive performance of GraDK is further demonstrated through a
comprehensive evaluation. Apart from assessing the fitted log-likelihood
($\ell$) of the testing data, for each data set, we generate $100$ event
sequences using each learned model (one of three independent runs) and provide
two metrics [24]: (i) the mean absolute error of predicted event frequency
(Time MAE) compared to that in the testing data, and (ii) the Kullback–Leibler
Divergence of predicted event types (Type KLD), which compares the empirical
distributions of event types (nodes) in the testing data and generated
sequences. These metrics reflect the model’s predictive capacity for future
events, as opposed to individual event prediction accuracy, which tends to be
noisy when applied to large graphs. The quantitative results in Table 1
demonstrate that the GraDK method excels in fitting sequential data on a
latent graph. It achieves the highest log-likelihood across all datasets and
significantly outperforms all baseline methods in predicting future events,
which holds immense importance within the domain of point process modeling.

##### Dynamic graph influence.

One salient feature of the proposed deep graph kernel is the capability of
recovering dynamic interactions among graph nodes, i.e., the evolution of
event influence over the graph. In Figure 2, we visualized the true non-
separable diffusion kernel used in the synthetic experiments (first row) and
the learned kernel (second row) over the graph with the time lag
$t-t^{\prime}$ fixed at three values (the diffusion kernel is stationary;
thus, we can fix $t^{\prime}$ at an arbitrary value without impacting the
results). The interactions among graph nodes vary over time, and by using our
deep influence kernel, we can successfully recover the dynamic nodal influence
out of discrete point process observations. The quantitative results in Table
1 also showcase the superior model performance of GraDK to capture the dynamic
and non-separable event influence.

##### Comparison of NLL and LS.

The empirical results indicate that both loss functions can be well-suited for
learning graph point processes. In particular, models with LS show
consistently better performance by a slight margin. In terms of the model
complexity, we show in Appendix E that both loss functions enjoy efficient
computation of complexity $\mathcal{O}(n)$, where $n$ is the total number of
events, and their computation times are similar.

####  4.1.2 Unknown graph topology

We now justify the capability of the proposed framework to recover the
structure of the event dependency when we have no prior information about the
latent graph structure. This is achieved by incorporating GAT into the graph
kernel and uncovering the graph support through the learned attention weights.
In particular, we exploit a fully connected graph support in the experiments
of GraDK+GAT+NLL, assuming possible interactions between any pair of graph
nodes. The estimated localized graph filters can indicate the knowledge about
the node interactions that the model has learned from the data. As we can see
in Figure 1, the recovered graph kernels and event dependencies by
GraDK+GAT+NLL closely resemble the ground truth. The model’s capability to
capture the underlying patterns of node interactions or dependencies is
crucial in various real-world applications in which one does not have access
to the latent graph structure.

###  4.2 Real data

We test the performance of our proposed approach on real-world point process
data. Since the applications of graph point processes involve discrete events
over networks with asynchronous time information, most of the traditional
benchmark graph data sets are not applicable. In the following, we collect
four publicly available data sets for numerical experiments: (i) traffic
congestion data in Atlanta, (ii) Sepsis data by PhysioNet Challenge, (iii)
wildfire data in California, and (iv) theft data in Valencia, Spain. Details
of the real data sets can be found in Appendix H.

####  4.2.1 Ablation study

For real-world applications involving discrete event data observed in
geographic space, spatio-temporal point processes (STPPs) with Euclidean-
distance-based influence kernels can be used. Nevertheless, through an
ablation study where we compare GraDK with three STPP baselines, we showcase
that our proposed graph point process model is more flexible and powerful in
capturing the more complicated event influence across discretized locations
(represented by nodes) that are distant but strongly influence each other.
Such influences are captured by direct edges between nodes or through multiple
hops, which can be harder to capture using spatio-temporal kernels.

Table 2: Ablation study: Comparison with Spatio-Temporal Point Processes
(STPPs).

|  | Wildfire (25 nodes) | Theft (52 nodes)
---|---|---|---
Model | #parameters | Testing $\ell$ | Time MAE | Type KLD | Testing $\ell$ | Time MAE | Type KLD
MHP | $625$ | $-3.846_{(0.003)}$ | $1.103$ | $0.072$ | $-3.229_{(0.005)}$ | $1.912$ | $0.483$
ETAS | $2$ | $-3.702_{(0.002)}$ | $1.134$ | $0.385$ | $-3.049_{(0.001)}$ | $3.750$ | $0.685$
DNSK-stpp | $4742$ | $\textbf{--3.647}_{(0.005)}$ | $0.861$ | $0.214$ | $-3.004_{(0.002)}$ | $1.342$ | $0.600$
GraDK | $411$ | $-3.650_{(0.002)}$ | 0.580 | 0.018 | $\textbf{--2.998}_{(0.002)}$ | 0.127 | 0.292

  * •

*Numbers in parentheses are standard errors for three independent runs.

Table 3: Real data results.

| Traffic congestion (5 nodes) | Sepsis (13 nodes) | Wildfire (25 nodes) | Theft (52 nodes)
---|---|---|---|---
Model | Testing $\ell$ | Time MAE | Type KLD | Testing $\ell$ | Time MAE | Type KLD | Testing $\ell$ | Time MAE | Type KLD | Testing $\ell$ | Time MAE | Type KLD
RMTPP | $-5.197_{(0.662)}$ | $2.348$ | $0.021$ | $-7.395_{(0.613)}$ | $3.159$ | $1.124$ | $-6.155_{(1.589)}$ | $1.180$ | $0.178$ | $-11.496_{(1.474)}$ | $5.871$ | $0.124$
FullyNN | $-3.292_{(0.108)}$ | $0.511$ | $0.012$ | $-4.666_{(0.133)}$ | $0.302$ | $0.345$ | $-4.717_{(0.119)}$ | $0.817$ | $0.026$ | $-3.468_{(0.068)}$ | $6.457$ | $1.169$
DNSK-mtpp | $-2.401_{(0.011)}$ | $0.934$ | $0.010$ | $-4.088_{(0.012)}$ | $0.712$ | $0.057$ | $-3.706_{(0.008)}$ | $0.711$ | $0.083$ | $-3.347_{(0.012)}$ | $0.507$ | $0.177$
THP-S | $-2.254_{(0.007)}$ | $0.378$ | $0.003$ | $-3.543_{(0.008)}$ | $0.113$ | $0.067$ | $-4.523_{(0.018)}$ | $1.183$ | $0.134$ | $-2.982_{(<0.001)}$ | $0.739$ | $0.189$
SAHP-G | $-2.453_{(0.013)}$ | $0.729$ | $0.021$ | $\textbf{--3.134}_{(0.106)}$ | $0.858$ | $1.061$ | $-3.919_{(0.040)}$ | $0.551$ | $0.032$ | $\textbf{--2.970}_{(0.032)}$ | $0.464$ | $0.096$
GraDK+L3net+NLL | $-2.178_{(0.005)}$ | $0.314$ | 0.001 | $-3.303_{(0.008)}$ | $0.082$ | $0.018$ | $\textbf{--3.625}_{(0.002)}$ | 0.207 | 0.006 | $-2.980_{(0.003)}$ | $0.640$ | $0.079$
GraDK+L3net+LS | $\textbf{--2.159}_{(0.004)}$ | 0.247 | 0.001 | $-3.225_{(0.005)}$ | 0.036 | 0.011 | $-3.628_{(0.002)}$ | $0.347$ | $0.013$ | $-2.982_{(0.004)}$ | 0.391 | 0.067
GraDK+GAT+NLL | $-2.281_{(0.011)}$ | $0.356$ | $0.015$ | $-3.556_{(0.009)}$ | $0.110$ | $0.025$ | $-3.629_{(0.007)}$ | $0.898$ | $0.085$ | $-2.995_{(0.006)}$ | $0.942$ | $0.173$

  * •

*Numbers in parentheses are standard errors for three independent runs.

We compare GraDK against three STPP baselines including (i) Multivariate
Hawkes process (MHP) [22]; (ii) Epidemic Type Aftershock Sequence (ETAS) model
[31]; and (iii) Deep Non-Stationary Kernel with a low-rank spatio-temporal
kernel (DNSK-stpp). Each baseline represents one type of approach for modeling
the spatial effect of historical events. MHP discretizes the entire space by
several geographic units and models the dependency among units using a spatial
kernel matrix. On the contrary, ETAS and DNSK-stpp adopt parametric and
neural-network-based spatial kernels and model the event dependencies over
continuous Euclidean space, respectively. We estimate each of the four models
with a fixed exponentially decaying temporal kernel on wildfire and theft data
sets which are originally observed within Euclidean geographic space.

Table 2 presents the quantitative results of three metrics used in synthetic
data sets for each model, which demonstrate the superior performance of GraDK
in fitting the data and predicting future events. It is worth noting that the
performance gain of our model does not rely on the increasing of parameters,
indicating the model benefits of the proposed graph kernel framework.

####  4.2.2 Comparison with baselines

Figure 3: Learned graph kernels for the theft data set evaluated at
$t^{\prime}=0$ and $t-t^{\prime}=1$; our proposed method can capture complex
inter-node dependence compared with prior work DNSK using a spatio-temporal
kernel. Figure 4: Dynamic graph influence learned by GraDK on the Sepsis data.
The learned kernel is evaluated at four time lags $2,4,6$, and $8$. The node
radii are proportional to the background intensity on each node, and the edge
widths are proportional to the influence magnitude. Blue and red edges
represent the excitation and inhibition effects, respectively. The history
time $t^{\prime}$ is set to be $0$ for all four panels.

Results in Table 3 underscore the efficacy of the GraDK approach in acquiring
knowledge about graph point processes across a diverse array of real-world
domains. These settings cover diverse event dependency dynamics, as the
influence mechanisms include infrastructure (roadways for traffic patterns),
healthcare (Granger causality among Sepsis-associated derangements), nature
(weather and climate for wildfire patterns), and social dynamics (criminal
behavior for theft patterns). Despite the inherent complexity of these
scenarios, our method excels in providing a robust framework capable of
capturing the intricate dependencies and facilitating accurate predictions,
demonstrated by the low Time MAE and Type KLD from our method in each setting,
which is better than or comparable to the best baselines in each of the four
real data sets. Note that we adopt GraDK+GAT+NLL to learn the latent graph
structure for each real data set and evaluate the models using the recovered
graph supports, which are presented in Appendix H.

Figure 3 visualizes the learned graph kernels by (a) GraDK+L3net+NLL, (b)
GraDK+GAT+NLL, (c) SAHP-G, and (d) DNSK on the theft data. The third panel
reveals that SAHP-G learns a noisy graph kernel, resulting in a conditional
intensity that depends slightly on inter-event influence, thus learning a
homogeneous Poisson process for each node with a relatively high likelihood.
DNSK (last panel) fails to present meaningful or discernible patterns of self
or inter-node dependencies. Lastly, GraDK+GAT captures intricate self-
influence and inter-node dependencies, and GraDK+L3net recovers multi-hop
event influence with the aid of the flexible graph kernel, indicating the
complex inhomogeneous dynamics among data.

####  4.2.3 Dynamic graph influence on Sepsis data

We use the ICU time-series data for sepsis prediction to demonstrate the model
capability of learning dynamic graph influence. Figure 4 visualizes four
snapshots of the influence learned by GraDK among different medical indices.
The snapshots are arranged in an order with the increase of the time lag
$t-t^{\prime}$, showcasing the temporal evolution of the graph influence. As
we can see, the excitation effects on Sepsis onsets by the observations of
certain medical indices (e.g., RI, OCD, DCO, OD(vs), Inf(vs)) are stronger in
the near future (e.g., at $t-t^{\prime}=2$) and would decay over time. Some
medical indices (e.g., Shock, and OD(Lab)) would have more interactions when
the Sepsis does not occur. Similar decaying patterns exist among the
inhibition effects between certain index pairs (e.g., OCD to EI, OCD to RI,
Inf(Lab) to RI). Full names of the medical indices can be found in Table H5.

##  5 Conclusion

We develop a novel deep kernel for graph point processes that can incorporate
various GNN architectures through localized graph filters, which efficiently
learn intricate and non-stationary event dynamics over latent graphs. The
kernel modeling enhances interpretability, as one can parse the learned kernel
to understand the dynamic event type interdependence. Our approach outperforms
existing methods in terms of dependency recovery and event prediction in both
simulation and real-data experiments.

## Acknowledgement

The work is supported by NSF DMS-2134037. Z.D., M.R. and Y.X. are partially
supported by an NSF CAREER CCF-1650913, NSF DMS-2134037, CMMI-2015787,
CMMI-2112533, DMS-1938106, and DMS-1830210. X.C. is partially supported by NSF
CAREER DMS-2237842, Simons Foundation 814643, and the Alfred P. Sloan
Foundation.

## References

  * Atwood and Towsley, [2016] Atwood, J. and Towsley, D. (2016).  Diffusion-convolutional neural networks.  Advances in Neural Information Processing Systems, 29.
  * Bacry et al., [2020] Bacry, E., Bompaire, M., Gaïffas, S., and Muzy, J.-F. (2020).  Sparse and low-rank multivariate hawkes processes.  Journal of Machine Learning Research, 21(50):1–32.
  * Bianchi et al., [2021] Bianchi, F. M., Grattarola, D., Livi, L., and Alippi, C. (2021).  Graph neural networks with convolutional arma filters.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 44(7):3496–3507.
  * Brody et al., [2021] Brody, S., Alon, U., and Yahav, E. (2021).  How attentive are graph attention networks?  In International Conference on Learning Representations.
  * Bruna et al., [2014] Bruna, J., Zaremba, W., Szlam, A., and LeCun, Y. (2014).  Spectral networks and deep locally connected networks on graphs.  In International Conference on Learning Representations, ICLR 2014.
  * Cai et al., [2022] Cai, B., Zhang, J., and Guan, Y. (2022).  Latent network structure learning from high-dimensional multivariate point processes.  Journal of the American Statistical Association, pages 1–14.
  * [7] Chen, R. T., Amos, B., and Nickel, M. (2020a).  Neural spatio-temporal point processes.  In International Conference on Learning Representations.
  * [8] Chen, Y., Wu, L., and Zaki, M. (2020b).  Iterative deep graph learning for graph neural networks: Better and robust node embeddings.  Advances in neural information processing systems, 33:19314–19326.
  * Cheng et al., [2021] Cheng, X., Miao, Z., and Qiu, Q. (2021).  Graph convolution with low-rank learnable local filters.  In International Conference on Learning Representations.
  * Daley and Vere-Jones, [2007] Daley, D. J. and Vere-Jones, D. (2007).  An introduction to the theory of point processes: volume II: general theory and structure.  Springer Science & Business Media.
  * Defferrard et al., [2016] Defferrard, M., Bresson, X., and Vandergheynst, P. (2016).  Convolutional neural networks on graphs with fast localized spectral filtering.  Advances in Neural Information Processing Systems, 29.
  * [12] Dong, Z., Cheng, X., and Xie, Y. (2023a).  Spatio-temporal point processes with deep non-stationary kernels.  In International Conference on Learning Representations.
  * [13] Dong, Z., Zhu, S., Xie, Y., Mateu, J., and Rodríguez-Cortés, F. J. (2023b).  Non-stationary spatio-temporal point process modeling for high-resolution covid-19 data.  Journal of the Royal Statistical Society Series C: Applied Statistics, 72(2):368–386.
  * Du et al., [2016] Du, N., Dai, H., Trivedi, R., Upadhyay, U., Gomez-Rodriguez, M., and Song, L. (2016).  Recurrent marked temporal point processes: Embedding event history to vector.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data mining, pages 1555–1564.
  * Dwivedi and Bresson, [2021] Dwivedi, V. P. and Bresson, X. (2021).  A generalization of transformer networks to graphs.  In AAAI’21 Workshop on Deep Learning on Graphs.
  * Fang et al., [2023] Fang, G., Xu, G., Xu, H., Zhu, X., and Guan, Y. (2023).  Group network hawkes process.  Journal of the American Statistical Association, (just-accepted):1–78.
  * Fatemi et al., [2021] Fatemi, B., El Asri, L., and Kazemi, S. M. (2021).  Slaps: Self-supervision improves structure learning for graph neural networks.  Advances in Neural Information Processing Systems, 34:22667–22681.
  * Gao et al., [2020] Gao, X., Hu, W., and Guo, Z. (2020).  Exploring structure-adaptive graph learning for robust semi-supervised classification.  In 2020 ieee international conference on multimedia and expo (icme), pages 1–6. IEEE.
  * Geer, [2000] Geer, S. A. (2000).  Empirical processes in M-estimation, volume 6.  Cambridge University Press.
  * Goldberger et al., [2000] Goldberger, A. L., Amaral, L. A., Glass, L., Hausdorff, J. M., Ivanov, P. C., Mark, R. G., Mietus, J. E., Moody, G. B., Peng, C.-K., and Stanley, H. E. (2000).  Physiobank, physiotoolkit, and physionet: components of a new research resource for complex physiologic signals.  circulation, 101(23):e215–e220.
  * Hansen et al., [2015] Hansen, N. R., Reynaud-Bouret, P., and Rivoirard, V. (2015).  Lasso and probabilistic inequalities for multivariate point processes.
  * Hawkes, [1971] Hawkes, A. G. (1971).  Spectra of some self-exciting and mutually exciting point processes.  Biometrika, 58(1):83–90.
  * Jin et al., [2020] Jin, W., Ma, Y., Liu, X., Tang, X., Wang, S., and Tang, J. (2020).  Graph structure learning for robust graph neural networks.  In Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining, pages 66–74.
  * Juditsky et al., [2020] Juditsky, A., Nemirovski, A., Xie, L., and Xie, Y. (2020).  Convex parameter recovery for interacting marked processes.  IEEE Journal on Selected Areas in Information Theory, 1(3):799–813.
  * Kingma and Ba, [2014] Kingma, D. P. and Ba, J. (2014).  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * Kipf and Welling, [2016] Kipf, T. N. and Welling, M. (2016).  Semi-supervised classification with graph convolutional networks.  In International Conference on Learning Representations.
  * Longa et al., [2023] Longa, A., Lachi, V., Santin, G., Bianchini, M., Lepri, B., Lio, P., Scarselli, F., and Passerini, A. (2023).  Graph neural networks for temporal graphs: State of the art, open challenges, and opportunities.  Transactions on Machine Learning Research.
  * Mei and Eisner, [2017] Mei, H. and Eisner, J. M. (2017).  The neural hawkes process: A neurally self-modulating multivariate point process.  Advances in Neural Information Processing Systems, 30.
  * Mercer, [1909] Mercer, J. (1909).  Xvi. functions of positive and negative type, and their connection the theory of integral equations.  Philosophical Transactions of the Royal Society of London. Series A, containing papers of a mathematical or physical character, 209(441-458):415–446.
  * Mollenhauer et al., [2020] Mollenhauer, M., Schuster, I., Klus, S., and Schütte, C. (2020).  Singular value decomposition of operators on reproducing kernel hilbert spaces.  In Advances in Dynamics, Optimization and Computation: A volume dedicated to Michael Dellnitz on the occasion of his 60th birthday, pages 109–131. Springer.
  * Ogata, [1988] Ogata, Y. (1988).  Statistical models for earthquake occurrences and residual analysis for point processes.  Journal of the American Statistical Association, 83(401):9–27.
  * Okawa et al., [2021] Okawa, M., Iwata, T., Tanaka, Y., Toda, H., Kurashima, T., and Kashima, H. (2021).  Dynamic hawkes processes for discovering time-evolving communities’ states behind diffusion processes.  In Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pages 1276–1286.
  * Omi et al., [2019] Omi, T., Aihara, K., et al. (2019).  Fully neural network based model for general temporal point processes.  Advances in Neural Information Processing Systems, 32.
  * Pan et al., [2023] Pan, Z., Wang, Z., and Zhe, S. (2023).  Graph-informed neural point process with monotonic nets.
  * Rampášek et al., [2022] Rampášek, L., Galkin, M., Dwivedi, V. P., Luu, A. T., Wolf, G., and Beaini, D. (2022).  Recipe for a general, powerful, scalable graph transformer.  Advances in Neural Information Processing Systems, 35:14501–14515.
  * Reinhart, [2018] Reinhart, A. (2018).  A review of self-exciting spatio-temporal point processes and their applications.  Statistical Science, 33(3):299–318.
  * Reyna et al., [2019] Reyna, M. A., Josef, C., Seyedi, S., Jeter, R., Shashikumar, S. P., Westover, M. B., Sharma, A., Nemati, S., and Clifford, G. D. (2019).  Early prediction of sepsis from clinical data: the physionet/computing in cardiology challenge 2019.  In 2019 Computing in Cardiology (CinC), pages Page–1. IEEE.
  * Shang and Sun, [2019] Shang, J. and Sun, M. (2019).  Geometric hawkes processes with graph convolutional recurrent neural networks.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 4878–4885.
  * Shchur et al., [2019] Shchur, O., Biloš, M., and Günnemann, S. (2019).  Intensity-free learning of temporal point processes.  In International Conference on Learning Representations.
  * Vaswani et al., [2017] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., Kaiser, Ł., and Polosukhin, I. (2017).  Attention is all you need.  Advances in Neural Information Processing Systems, 30.
  * Velickovic et al., [2017] Velickovic, P., Cucurull, G., Casanova, A., Romero, A., Lio, P., Bengio, Y., et al. (2017).  Graph attention networks.  stat, 1050(20):10–48550.
  * Wei et al., [2023] Wei, S., Xie, Y., Josef, C. S., and Kamaleswaran, R. (2023).  Granger causal chain discovery for sepsis-associated derangements via continuous-time hawkes processes.  In Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pages 2536–2546.
  * [43] Wu, W., Liu, H., Zhang, X., Liu, Y., and Zha, H. (2020a).  Modeling event propagation via graph biased temporal point process.  IEEE Transactions on Neural Networks and Learning Systems.
  * [44] Wu, Z., Pan, S., Chen, F., Long, G., Zhang, C., and Philip, S. Y. (2020b).  A comprehensive survey on graph neural networks.  IEEE Transactions on Neural Networks and Learning Systems, 32(1):4–24.
  * Xia et al., [2021] Xia, F., Sun, K., Yu, S., Aziz, A., Wan, L., Pan, S., and Liu, H. (2021).  Graph learning: A survey.  IEEE Transactions on Artificial Intelligence, 2(2):109–127.
  * Xia et al., [2022] Xia, W., Li, Y., and Li, S. (2022).  Graph neural point process for temporal interaction prediction.  IEEE Transactions on Knowledge and Data Engineering.
  * Xu et al., [2016] Xu, H., Farajtabar, M., and Zha, H. (2016).  Learning granger causality for hawkes processes.  In International conference on machine learning, pages 1717–1726. PMLR.
  * Yang et al., [2022] Yang, C., Mei, H., and Eisner, J. (2022).  Transformer embeddings of irregularly spaced events and their participants.  In International Conference on Learning Representations.
  * Yuan et al., [2023] Yuan, Y., Ding, J., Shao, C., Jin, D., and Li, Y. (2023).  Spatio-temporal diffusion point processes.  In ACM SIGKDD Conference on Knowledge Discovery and Data Mining.
  * Zhang et al., [2020] Zhang, Q., Lipani, A., Kirnap, O., and Yilmaz, E. (2020).  Self-attentive hawkes process.  In International Conference on Machine Learning, pages 11183–11193. PMLR.
  * Zhang et al., [2021] Zhang, Q., Lipani, A., and Yilmaz, E. (2021).  Learning neural point processes with latent graphs.  In Proceedings of the Web Conference 2021, pages 1495–1505.
  * [52] Zhou, K., Zha, H., and Song, L. (2013a).  Learning social infectivity in sparse low-rank networks using multi-dimensional hawkes processes.  In Artificial Intelligence and Statistics, pages 641–649. PMLR.
  * [53] Zhou, K., Zha, H., and Song, L. (2013b).  Learning triggering kernels for multi-dimensional hawkes processes.  In International conference on machine learning, pages 1301–1309. PMLR.
  * Zhou et al., [2022] Zhou, Z., Yang, X., Rossi, R., Zhao, H., and Yu, R. (2022).  Neural point process for learning spatiotemporal event dynamics.  In Learning for Dynamics and Control Conference, pages 777–789. PMLR.
  * Zhu and Koniusz, [2021] Zhu, H. and Koniusz, P. (2021).  Simple spectral graph convolution.  In International Conference on Learning Representations.
  * [56] Zhu, S., Ding, R., Zhang, M., Van Hentenryck, P., and Xie, Y. (2021a).  Spatio-temporal point processes with attention for traffic congestion event modeling.  IEEE Transactions on Intelligent Transportation Systems, 23(7):7298–7309.
  * [57] Zhu, S., Li, S., Peng, Z., and Xie, Y. (2021b).  Imitation learning of neural spatio-temporal point processes.  IEEE Transactions on Knowledge and Data Engineering, 34(11):5391–5402.
  * Zhu et al., [2022] Zhu, S., Wang, H., Dong, Z., Cheng, X., and Xie, Y. (2022).  Neural spectral marked point processes.  In International Conference on Learning Representations.
  * Zhu and Xie, [2022] Zhu, S. and Xie, Y. (2022).  Spatiotemporal-textual point processes for crime linkage detection.  The Annals of Applied Statistics, 16(2):1151–1170.
  * [60] Zhu, S., Zhang, M., Ding, R., and Xie, Y. (2021c).  Deep fourier kernel for self-attentive point processes.  In International Conference on Artificial Intelligence and Statistics, pages 856–864. PMLR.
  * [61] Zhu, Y., Xu, W., Zhang, J., Du, Y., Zhang, J., Liu, Q., Yang, C., and Wu, S. (2021d).  A survey on graph structure learning: Progress and opportunities.  arXiv preprint arXiv:2103.03036.
  * Zuo et al., [2020] Zuo, S., Jiang, H., Li, Z., Zhao, T., and Zha, H. (2020).  Transformer hawkes process.  In International Conference on Machine Learning, pages 11692–11702. PMLR.

## Appendix

##  Appendix A Comparison with existing approaches

To further elaborate on the contributions of our paper, we provide a detailed
comparison between our method and existing approaches that study the point
processes in the context of graphs, neural networks, and graph neural
networks. Table A1 compares our GraDK with existing methods in terms of four
aspects that play crucial roles in establishing an effective and high-
performing model for marked temporal point processes over graphs, which are
also the main contributions of our paper that distinguish the proposed method
from those in previous studies.

Table A1: Comparison between our method GraDK with other point processes with
graphs, neural networks, and graph neural networks.

Model | No parametric constraint | Modeling influence kernel | Joint modeling event time and mark on graph | Using GNN
---|---|---|---|---
1\. GNHP [16] | ✗ | ✓ | ✓ | ✗
2\. GINPP [34] | ✓ | ✗ | ✗ | ✗
3\. RMTPP [14] | ✗ | ✗ | ✗ | ✗
4\. NH [28] | ✓ | ✗ | ✗ | ✗
5\. FullyNN [33] | ✓ | ✗ | ✗ | ✗
6\. SAHP [50] | ✓ | ✗ | ✗ | ✗
7\. DNSK [12] | ✓ | ✓ | ✗ | ✗
8\. LogNormMix [39] | ✓ | ✗ | ✗ | ✗
9\. DeepSTPP [54] | ✗ | ✗ | ✗ | ✗
10\. NSTPP [7] | ✓ | ✗ | ✗ | ✗
11\. DSTPP [49] | ✓ | ✗ | ✗ | ✗
12\. THP-S [62] | ✓ | ✗ | ✓ | ✗
13\. SAHP-G [51] | ✓ | ✗ | ✓ | ✗
14\. GHP [38] | ✗ | ✓ | ✗ | ✓
15\. GNPP [46] | ✓ | ✗ | ✗ | ✓
16\. GBTPP [43] | ✗ | ✗ | ✗ | ✓
17\. GraDK (Ours) | ✓ | ✓ | ✓ | ✓

We have compared our model with Methods 3, 5, 7, 12, and 13 in the above table
by experiments and showed our advantageous performance. The comparisons
between Methods 3, 5, and 7 demonstrate the benefits of considering latent
graph information when modeling graph point processes. The comparisons between
Methods 12 and 13 (they are deep neural point processes that incorporate known
graph structures through the adjacency matrix) showcase the advantages of
incorporating GNNs for flexibly capturing the complicated multi-hop, non-
Euclidean event influence. We do not compare with (i) Methods 1 and 14, for
they aim to develop theoretical results in the context of parametric point
processes with graph/network structures; (ii) Method 2, for it is an extension
of Method 5; (iii) Methods 4 and 8 for their similarity with Methods 3 and 5,
respectively; (iv) Method 6 for it is another version of Method 13 without
considering the graph structure; (v) Methods 9, 10, and 11 for they focus on
spatio-temporal point processes, and we have demonstrated the advantages of
our graph point process model against STPPs using ablation study in Section
4.2.1; (vi) Methods 15 and 16 for they focus on the task of temporal link
prediction rather than modeling discrete marked events.

Being the first paper that incorporates GNNs in the influence kernel without
any parametric constraints, four aspects of our paper’s contributions include
(also highlighted in Section 1):

  1. 1.

We have no parametric model constraints, which becomes increasingly important
in modern applications for the model to capture complex event dependencies.
Prior works (e.g., Methods 1, 3, 14) typically assume a specific parametric
form for the model, although they can be data- and computationally efficient,
which may restrict the model expressiveness for modern applications with
complex and large datasets.

  2. 2.

We model the influence kernel instead of the conditional intensity, which
effectively leverages the repeated influence patterns by combining statistical
and deep models with excellent model interpretability. As shown in Table A1,
most previous studies fall into the category of modeling the conditional
intensity function instead of the influence kernel. While such models are
expressive, they ignore the benefits of leveraging the repeated influence
patterns, and the model interpretability may be limited.

  3. 3.

We jointly model event time and mark over a latent graph structure, which most
existing approaches fail to achieve. Some approaches (e.g., Methods 5, 8)
enjoy great model flexibility for modeling event times, but they only focus on
temporal point processes and can hardly be extended to mark space. Another
body of research (e.g., Methods 2, 3, 8) decouples the time and mark space
when modeling event distributions for model simplicity and efficiency, which
sometimes limits the model flexibility and leads to sub-optimal solutions.
More studies model marked events without any graph structure considered. One
example is the family of existing spatio-temporal point process approaches
(e.g., Methods 9, 10, 11) that are inapplicable in the setting of graph point
processes.

  4. 4.

We adopt GNNs in the influence kernel, which faces non-trivial challenges and
requires a carefully designed model architecture. As shown in Table A1, there
has been limited work in exploiting underlying graph structures by leveraging
GNNs. Previous studies (e.g., Methods 2, 12, 13) combine non-graph neural
network architectures (e.g., fully-connected neural networks, RNNs, etc.)
along with certain graph information (e.g., adjacency matrices), rather than
directly incorporating GNN structures in the models. Another attempt (Method
14) integrates graph convolutional RNNs with the influence kernel. However,
they focus on using GNNs to estimate the parameters of the conditional
intensity function with a parametric (exponentially decaying) kernel. Two
previous studies (Methods 15, 16) on modeling point processes resort to
message-passing GNNs. However, they do not use influence kernels and focus on
a different problem of graph node interaction prediction rather than modeling
the dependencies and occurrences of discrete marked events.

It is worth noting that extending the influence kernel onto the graph
structure with GNN architectures incorporated is crucial and non-trivial for
modeling point processes over graphs. Because (i) the lack of a
straightforward definition of distance between event marks makes the distance-
based spatial kernel non-sensible on a graph structure, and (ii) simply
replacing the distance-based kernel with some scalar coefficients to capture
the interaction between each pair of nodes will significantly limit the model
expressiveness, especially in modern applications. In this paper, we design
the influence kernel motivated by kernel SVD and Mercer’s theorem [29] and
innovatively introduce the concept of localized graph filter basis into the
kernel to fully take advantage of the universal representation power of GNNs.

##  Appendix B Example: Graph filter bases in L3Net

Figure B1 illustrates the mechanism of the graph filter bases in L3Net to
capture the event dependencies during the modeling of sequential events on an
8-node graph. The neighborhood orders are highlighted as the superscripts of
each graph filter basis.

(a) The latent graph structure

(b) Events dependencies over graph nodes

Figure B1: An example of modeling events on an 8-node graph using graph filter
bases in L3Net: (a) The latent graph structure. Blue and red nodes represent
the $1$st and $2$nd order neighbors of $v_{0}$, denoted by $N_{v_{0}}^{(1)}$
and $N_{v_{0}}^{(2)}$, respectively. (b) Three graph filter bases $B^{(0)}$,
$B^{(1)}$, and $B^{(2)}$ capture the dependencies between events. Hollow
circles represent events observed on each node. Colored lines indicate the
potential influence of the earliest type-$v_{0}$ event on future events
captured by different bases.

##  Appendix C Incorporation of localized graph filters in graph kernel

We provide the details of incorporating different localized graph filters in
graph convolutions into our proposed graph basis kernel, including the model
formulation and complexity (number of trainable parameters). A graph
convolution layer maps the input node feature $X$ to the output feature $Y$.
Particularly, for the graph convolutions in the graph neural networks
discussed in Section 3.3, we have:

  1. (1)

Chebnet [11]: The filtering operation in Chebnet is defined as

$\vspace{-0.05in}Y=\sum_{r=0}^{R-1}\theta_{r}T_{r}(\tilde{L})X,$

where $T_{r}$ is the Chebyshev polynomial of order $r$ evaluated at the scaled
and normalized graph Laplacian $\tilde{L}=2L/\lambda_{\text{max}}-I$, where
$L=I-D^{-1/2}AD^{-1/2}$, $D$ is the degree matrix, $A$ is the adjacency
matrix, and $\lambda_{\text{max}}$ is the largest eigenvalue of $L$.

We let $B_{r}=T_{r-1}(\tilde{L})\in\mathbb{R}^{|V|\times|V|}$, which can be
computed by the recurrence relation of Chebyshev polynomial
$B_{r}=2LB_{r-1}-B_{r-2}$. Since $\\{B_{r}\\}_{r=1}^{R}$ are calculated by the
graph Laplacian, the trainable parameters are polynomial coefficients (in our
case, the weights that combine basis kernels) of $\mathcal{O}(R)$.

  2. (2)

L3Net [9]: In L3Net, the graph convolution layer for single channel graph
signal is specified as (omitting the bias, non-linear mapping)

$\vspace{-0.1in}Y=\sum_{r=1}^{R}a_{r}B_{r}^{(o_{r})}(v^{\prime},v)X.$

Each $B_{r}^{(o_{r})}(v^{\prime},v)\neq 0$ only if $v\in
N_{v^{\prime}}^{(o_{r})}$.

In our model, we choose an integer $o_{r}$ for each $B_{r}$ and formulate them
the same way in L3Net. Thus, each $B_{r}(v^{\prime},v)$ models the influence
from $v^{\prime}$ to its $o_{r}$-th order neighbors on the graph. All bases
$\\{B_{r}\\}_{r=1}^{R}$ are learnable. Assuming the average number of
neighbors of graph nodes is $C$, the total number of trainable parameters is
of $\mathcal{O}(RC|V|)$.

  3. (3)

GAT [41]: Considering GAT with $R$ attention heads, the graph convolution
operator in one layer can be written as

$\vspace{-0.1in}Y=\sum_{r=1}^{R}\mathcal{A}^{(r)}X\Theta_{r},\quad\mathcal{A}_{v^{\prime},v}^{(r)}=\frac{e^{c_{v^{\prime}v}^{(r)}}}{\sum_{v\in
N_{v^{\prime}}^{(1)}}e^{c_{v^{\prime}v}^{(r)}}},\quad
c_{v^{\prime}v}^{(r)}=\sigma((a^{(r)})^{\top}[W^{(r)}X_{v^{\prime}}||W^{(r)}X_{v}]).$

Here $\\{a^{(r)},W^{(r)}\\}_{r=1}^{R}$ are trainable model parameters, and
$\Theta_{r}=W^{(r)}C^{(r)}$, where $C^{(r)}$ is a fixed matrix for
concatenating outputs from $R$ attention heads. Mask attention is applied to
inject the graph structure, that is, $\mathcal{A}_{v^{\prime},v}^{(r)}\neq 0$
only when $v\in N_{v^{\prime}}^{(1)}$.

Therefore, each $B_{r}(v^{\prime},v)$ in our graph basis kernel can be a
learnable localized graph filter with positive entries and column-wise
summation normalized to one to mimic the affinity matrix $\mathcal{A}^{(r)}$.
When global attention is allowed (i.e., every node can attend on every other
node), the number of trainable parameters is of $\mathcal{O}(R|V|^{2})$. When
mask attention is applied (which is commonly used in GAT), the total number of
trainable parameters is $\mathcal{O}(RC|V|)$.

  4. (4)

GPS Graph Transformer [35]: At each layer of GPS Graph Transformer, the graph
node features are updated by aggregating the output of a message-passing graph
neural network (MPNN) layer and a global attention layer. Assuming a sum
aggregator in the MPNN layer, the filter operation can be expressed as

$\vspace{-0.1in}Y=\sum_{r=1}^{R}\left(\mathcal{A}^{(r)}XW_{1}^{(r)}+XW_{2}^{(r)}\right)C^{(r)}.$

Here $\mathcal{A}^{(r)}$ is the affinity matrix in the attention layer, and
$W_{1}^{(r)},W_{2}^{(r)}$ are weight matrices in the attention layer and MPNN
layer, respectively. The fixed matrices $\\{C^{(r)}\\}$ concatenate the $R$
attention heads and MPNN layers.

We can integrate such a GPS Graph Transformer structure by introducing $2R$
localized graph filter bases. For each $r$, $B_{r}$ takes the form of a
learnable affinity matrix, and $B_{2r}=A$, where $A$ is the graph adjacency
matrix. Here $\\{B_{r}\\}_{r=1}^{R}$ are learnable and
$\\{B_{r}\\}_{r=R+1}^{2R}$ are fixed. The total number of trainable parameters
with mask attention adopted is of $\mathcal{O}(RC|V|)$.

##  Appendix D Derivation of two loss functions

Considering one event sequence with $n$ events, we present the detailed
derivation of loss functions NLL and LS in the following:

##### Negative log-likelihood (NLL).

The model likelihood of point processes can be derived from the conditional
intensity (1). For the $i+1$-th event at $(t,v)$, by definition, we can re-
write $\lambda(t,v)$ as follows:

$\displaystyle\lambda(t,v)$
$\displaystyle=\mathbb{E}\left(\mathbb{N}_{v}([t,t+\Delta t]\times
v))|\mathcal{H}_{t}\right)/dt$
$\displaystyle=\mathbb{P}\\{t_{i+1}\in[t,t+\Delta
t],v_{i+1}=v|\mathcal{H}_{t_{i+1}}\cup\\{t_{i+1}\geq t\\}\\}/dt$
$\displaystyle=\frac{\mathbb{P}\\{t_{i+1}\in[t,t+\Delta
t],v_{i+1}=v,t_{i+1}\geq
t,v_{i+1}=v|\mathcal{H}_{t_{i+1}}\\}/dt}{\mathbb{P}\\{t_{i+1}\geq
t|\mathcal{H}_{t_{i+1}}\\}}.$ $\displaystyle=\frac{f(t,v)}{1-F(t)}.$

Here $\mathbb{N}_{v}$ is the counting measure on node $v$,
$F(t)=\mathbb{P}(t_{i+1}<t|\mathcal{H}_{t_{i+1}})$ is the conditional
cumulative probability and $f(t,v)$ is the corresponding conditional
probability density for the next event happening at $(t,v)$ given observed
history. We multiply the differential $dt$ on both sides of the equation and
sum over all the graph nodes $V$:

$\displaystyle dt\cdot\sum_{v\in V}\lambda(t,v)=\frac{dt\cdot\sum_{v\in
V}f(t,v)}{1-F(t)}=-d\log{(1-F(t))}.$

Integrating over $t$ on $[t_{i},t)$ with $F(t_{i})=0$ leads to the fact that

$F(t)=1-\exp\left(-\int_{t_{i}}^{t}\sum_{v\in V}\lambda(t,v)dt\right),$

then we have

$f(t,v)=\lambda(t,v)\cdot\exp\left(-\int_{t_{i}}^{t}\sum_{v\in
V}\lambda(t,v)dt\right).$

Using the chain rule, we have the model log-likelihood to be:

$\displaystyle\ell_{\text{NLL}}$ $\displaystyle=-\log
f((t_{1},v_{1}),\dots,(t_{n},v_{n}))=-\log\prod_{i=1}^{n}f(t_{i},v_{i})$
$\displaystyle=\sum_{v\in
V}\int_{0}^{T}\lambda(t,v)dt-\sum_{i=1}^{n}\lambda(t_{i},v_{i}).$

##### Least square (LS) loss.

We expect that the integral of the conditional intensity over infinitesimal
time intervals containing event times should be approximately one, while
during non-event times, it should be approximately zero. Specifically,
involving each graph node’s counting process $\mathbb{N}_{v}$, we have

$\displaystyle\ell_{\text{LS}}$ $\displaystyle=\sum_{v\in
V}\int_{0}^{T}\left(\lambda(t,v)-1\right)^{2}d\mathbb{N}_{v}(t)+\int_{0}^{T}\left(\lambda(t,v)-0\right)^{2}(dt-d\mathbb{N}_{v}(t))$
$\displaystyle=\sum_{v\in
V}\int_{0}^{T}\lambda^{2}(t,v)d\mathbb{N}_{v}(t)-2\sum_{v\in
V}\int_{0}^{T}\lambda(t,v)d\mathbb{N}_{v}(t)+n+\sum_{v\in
V}\int_{0}^{T}\lambda^{2}(t,v)dt$ $\displaystyle\quad-\sum_{v\in
V}\int_{0}^{T}\lambda^{2}(t,v)d\mathbb{N}_{v}(t)$ $\displaystyle=\sum_{v\in
V}\int_{0}^{T}\lambda^{2}(t,v)dt-\sum_{i=1}^{n}2\lambda(t_{i},v_{i})+n.$

Omitting the constant, we can have the $\ell_{\text{LS}}$.

The log-likelihood of $|\mathcal{S}|$ observed event sequences can be
conveniently obtained by summing over the loss over all the sequences with
counting measure $\mathbb{N}_{v}$ replaced by $\mathbb{N}_{s}$ for $s$-th
sequence.

##  Appendix E Efficient model computation

In our optimization problem, the log-barrier term $p(\theta,b)$ guarantees the
validity of the learned models by penalizing the value of the conditional
intensity function on a dense enough grid over space, denoted as
$\mathcal{U}_{\text{bar},t}\times V$ where
$\mathcal{U}_{\text{bar},t}\subset[0,T]$. The optimization problems with two
loss functions are expressed as

$\displaystyle\min_{\theta}\mathcal{L}_{\text{NLL}}(\theta):=$ (E1)
$\displaystyle\frac{1}{|\mathcal{S}|}\sum_{s=1}^{|\mathcal{S}|}\left(\sum_{v\in
V}\int_{0}^{T}\lambda(t,v)dt-\sum_{i=1}^{n_{s}}\log\lambda\left(t_{i}^{s},v_{i}^{s}\right)\right)-\frac{1}{w|\mathcal{S}||\mathcal{U}_{\text{bar},t}\times
V|}\sum_{s=1}^{|\mathcal{S}|}\sum_{t\in\mathcal{U}_{\text{bar},t}}\sum_{v\in
V}\log(\lambda(t,v)-b),$

and

$\displaystyle\min_{\theta}\mathcal{L}_{\text{LS}}(\theta):=$ (E2)
$\displaystyle\frac{1}{|\mathcal{S}|}\sum_{s=1}^{|\mathcal{S}|}\left(\sum_{v\in
V}\int_{0}^{T}\lambda^{2}(t,v)dt-\sum_{i=1}^{n_{s}}2\lambda(t_{i}^{s},v_{i}^{s})\right)-\frac{1}{w|\mathcal{S}||\mathcal{U}_{\text{bar},t}\times
V|}\sum_{s=1}^{|\mathcal{S}|}\sum_{t\in\mathcal{U}_{\text{bar},t}}\sum_{v\in
V}\log(\lambda(t,v)-b),$

respectively. The computational complexity associated with calculating the
objective function has consistently posed a limitation for neural point
processes. The evaluations of neural networks are computationally demanding,
and traditional numerical computation requires model evaluations between every
pair of events. Consider one event sequence $\\{(t_{i},v_{i})\\}_{i=1}^{n}$ in
$\mathcal{S}$ with a total number of $n$ events; traditional methods have a
complexity of $\mathcal{O}(n^{2})$ for neural network evaluation. In the
following, we showcase our efficient model computation of complexity
$\mathcal{O}(n)$ for two different loss functions using a domain
discretization strategy.

##### Negative log-likelihood.

We identify three distinct components in the optimization objective (E1) as
log-summation, integral, and log-barrier, respectively. We introduce a uniform
grid $\mathcal{U}_{t}$ over time horizon $[0,T]$, and the computation for each
term can be carried out as follows:

  * •

Log-summation. We plug in the influence kernel (3) to the log-summation term
and have

$\sum_{i=1}^{n}\log\lambda(t_{i},v_{i})=\sum_{i=1}^{n}\log{\left(\mu+\sum_{t_{j}<t_{i}}\sum_{r=1}^{R}\sum_{l=1}^{L}\alpha_{lr}\psi_{l}(t_{j})\varphi_{l}(t_{i}-t_{j})B_{r}(v_{j},v_{i})\right)}.$

Each $\psi_{l}$ is only evaluated at event times $\\{t_{i}\\}_{i=1}^{n}$. To
prevent redundant evaluations of the function $\varphi_{l}$ for every pair of
event times $(t_{i},t_{j})$, we restrict the evaluation of $\varphi_{l}$ on
the grid $\mathcal{U}_{t}$. By adopting linear interpolation between two
neighboring grid points of $t_{i}-t_{j}$, we can determine the value of
$\varphi_{l}(t_{i}-t_{j})$. In practice, the influence of past events is
limited to a finite range, which can be governed by a hyperparameter
$\tau_{\text{max}}$. Consequently, we can let $\varphi_{l}(\cdot)=0$ when
$t_{i}-t_{j}>\tau_{\text{max}}$ without any neural network evaluation. The
computation of $B_{r}(v_{j},v_{i})$ is accomplished using matrix indexing, a
process that exhibits constant computational complexity when compared to the
evaluation of neural networks.

  * •

Integral. The efficient computation of the integral benefits from the
formulation of our conditional intensity function. We have

$\displaystyle\sum_{v\in V}\int_{0}^{T}\lambda(t,v)dt$
$\displaystyle=\mu|V|T+\sum_{i=1}^{n}\sum_{v\in
V}\int_{0}^{T}I(t_{i}<t)k(t_{i},t,v_{i},v)dt$
$\displaystyle=\mu|V|T+\sum_{i=1}^{n}\sum_{r=1}^{R}\sum_{v=1}^{V}B_{r}(v_{i},v)\sum_{l=1}^{L}\alpha_{rl}\psi_{l}(t_{i})\int_{0}^{T-t_{i}}\varphi_{l}(t)dt.$

Similarly, the evaluations of $B_{r}(v_{i},v)$ are the extractions of
corresponding entries in the matrices, and $\psi_{l}$ is only evaluated at
event times $\\{t_{i}\\}_{i=1}^{n}$. We leverage the evaluations of
$\varphi_{l}$ on the grid $\mathcal{U}_{t}$ to facilitate the computation of
the integral of $\varphi_{l}$. Let
$F_{l}(t):=\int_{0}^{t}\varphi_{l}(\tau)d\tau$. The value of $F_{l}$ on $j$-th
grid point in $\mathcal{U}_{t}$ equals the cumulative sum of $\varphi_{l}$
from the first grid point up to the $j$-th point, multiplied by the grid
width. Then, $F_{l}(T-t_{i})$ can be computed by the linear interpolation of
values of $F_{l}$ at two adjacent grid points of $T-t_{i}$.

  * •

Log-barrier. The computation can be carried out similarly to the computation
of the log-summation term by replacing $(t_{i},v_{i})$ with $(t,v)$ where
$t\in\mathcal{U}_{\text{bar},t},v\in V$.

##### Least square loss.

Similarly, we term the three components in the objective function (E2) as
square integral, summation, and log-barrier, respectively. The terms of
summation and log-barrier are calculated in the same way as the log-summation
and log-barrier in (E1), respectively, except for no logarithm in the
summation. Since expanding the integral after squaring the conditional
intensity function is complicated, we leverage the evaluations of the
intensity function on the dense enough grid for log-barrier penalty
$\mathcal{U}_{\text{bar},t}\times V$ and use numerical integration for
computing the square integral.

Table E2: Comparison of the computation time for two loss functions on each
synthetic data set.

Model |  3-node graph with negative influence |  16-node graph with 2-hop influence |  50-node graph |  50-node graph with non-separable kernel |  225-node graph
---|---|---|---|---|---
GraDK+NLL | $0.931$ | $3.993$ | $25.449$ | $30.348$ | $4.948$
GraDK+LS | $0.985$ | $3.927$ | $26.182$ | $29.912$ | $5.619$

  * •

*Unit: seconds per epoch.

We provide the analysis of the computational complexity of $\mathcal{O}(n)$.
For objective (E1), the evaluation of $\\{\psi_{l}\\}_{l=1}^{L}$ has a
complexity of $\mathcal{O}(Ln)$, while the evaluation of
$\\{\varphi_{l}\\}_{l=1}^{L}$ requires $\mathcal{O}(L|\mathcal{U}_{t}|)$
complexity. On the other hand, $\\{B_{r}\\}_{r=1}^{R}$ are computed with
$\mathcal{O}(1)$ complexity. Therefore, the total complexity of our model
computation is $\mathcal{O}(n+|\mathcal{U}_{t}|)$. Moreover, the grid
selection in practice is flexible, striking a balance between learning
accuracy and efficiency. For instance, the number of grid points in
$\mathcal{U}_{t}$ and $\mathcal{U}_{\text{bar},t}$ can be both chosen around
$\mathcal{O}(n)$, leading to an overall computational complexity of
$\mathcal{O}(n)$. Similar analysis can be carried out for objective (E2). It
is worth noting that all the evaluations of neural networks and localized
filter bases can be implemented beforehand, and both optimization objectives
are efficiently calculated using basic numerical operations. Table E2 shows
the wall clock times for model training with each loss function on all the
synthetic data sets, indicating that the computation time for NLL and LS are
similar.

##  Appendix F Algorithm

The weight $w$ and the lower bound $b$ need to be adjusted accordingly during
optimization to learn the valid model. For example, non-sensible solutions
(negative intensity) appear frequently at the early stage of the training
process, thus the lower bound $b$ should be close enough to the minimum of all
the evaluations on the grid to effectively steer the intensity functions
towards non-negative values, and the weight $w$ can be set as a small value to
magnify the penalty influence. When the solutions successfully reach the
neighborhood of the ground truth and start to converge (thus no negative
intensity would appear with a proper learning rate), the $b$ should be away
from the intensity evaluations (e.g., upper-bounded by $0$), and $w$ should be
large enough to remove the effect of the log-barrier penalty. We provide our
training algorithm using stochastic gradient descent in the following.

Input: Training set $X$, batch size $M$, epoch number $E$, learning rate
$\gamma$, constant $a>1$ to update $w$ in (E1) or (E2).

Initialization: model parameter $\theta_{0}$, first epoch $e=0$, $s=s_{0}$.
Conditional intensity lower bound $b_{0}$.

while $e<E$ do

Set $b_{temp}=0$.

for each batch with size $M$ do

1\. Compute $\ell(\theta)$,
$\\{\lambda(t_{c_{t}},s_{c_{s}})\\}_{c_{t}=1,\dots,|\mathcal{U}_{\text{bar},t}|,{c_{s}}=1,\dots,|V|}$.

2\. Compute $\mathcal{L}(\theta)=-\ell(\theta)+\frac{1}{w}p(\theta,b_{e})$.

3\. Update
$\theta_{e+1}\leftarrow\theta_{e}-\gamma\frac{\partial\mathcal{L}}{\partial\theta_{e}}$.

4\. Set
$b_{temp}=\min\left\\{\min\\{\\{\lambda(t_{c_{t}},s_{c_{s}})\\}_{c_{t}=1,\dots,|\mathcal{U}_{\text{bar},t}|,c_{s}=1,\dots,|\mathcal{U}_{\text{bar},s}|}-\epsilon,b_{temp}\right\\}$,
where $\epsilon$ is a small value to guarantee logarithm feasibility.

end for$e\leftarrow e+1,w\leftarrow w\cdot a,b_{e}=b_{temp}$

end while

Algorithm 1 Model learning

##  Appendix G Choice of kernel rank

The kernel rank can be chosen using two approaches. The first approach treats
the kernel rank as a model hyperparameter, and its selection is carried out
through cross-validation. This process aims to optimize the log-likelihood on
validation datasets. Each choice of $(L,R)\in[1,2]\times[1,2,3]$ has been
assessed for GraDK+L3net+NLL on synthetic and real-world datasets. The order
of the $r$-th filter is set to be $r$ (Tables G3 and G4). The log-likelihood
values in the synthetic data show only minor differences given different
hyperparameter choices once the parameters are chosen such that they
sufficiently capture temporal or multi-hop graph influence (typically
requiring $R>1$). In the real-world data, $L=1$ with $R=3$ yields the highest
log-likelihood on the validation dataset, which are the parameters used in
Section 4.2.

The rank of the kernel can also be viewed as a level of model complexity to be
learned from the data. Generally, under certain regularity assumptions of the
kernel in kernel SVD, the singular values will decay towards 0, resulting in a
low-rank kernel approximation. The singular values are captured by the
coefficients $\alpha_{rl}$. The kernel rank is chosen by preserving
significant coefficients while discarding higher-order ones; we can treat each
$\alpha_{rl}$ as one learnable parameter without the need to choose $D$. The
effectiveness of this approach is showcased by the results of the wildfire
data. We set $L=2$ and $R=3$ at first, and then fit the model, and visualize
the learned kernel coefficients in Figure G2, which suggests that $L=1$ and
$R=3$ would be an optimal choice for the kernel rank. This choice is used in
subsequent model fitting.

Table G3: Testing log-likelihood for GraDK with different kernel ranks.

| 16-node graph | 50-node graph | Wildfire | Theft
---|---|---|---|---
Kernel rank | $R=1$ | $R=2$ | $R=3$ | $R=1$ | $R=2$ | $R=3$ | $R=1$ | $R=2$ | $R=3$ | $R=1$ | $R=2$ | $R=3$
$L=1$ | $-3.032_{(<0.001)}$ | $-2.997_{(0.001)}$ | $-2.995_{(0.002)}$ | $-1.085_{(<0.001)}$ | $\textbf{--1.065}_{(0.002)}$ | $-1.067_{(0.001)}$ | $-3.900_{(0.004)}$ | $-3.654_{(0.009)}$ | $\textbf{--3.625}_{(0.002)}$ | $-3.405_{(0.001)}$ | $-3.039_{(0.011)}$ | $\textbf{--2.980}_{(0.003)}$
$L=2$ | $-3.030_{(0.001)}$ | $-2.992_{(0.002)}$ | $\textbf{--2.987}_{(0.005)}$ | $-1.085_{(0.001)}$ | $-1.067_{(0.002)}$ | $\textbf{--1.065}_{(0.002)}$ | $-3.977_{(0.005)}$ | $-3.657_{(0.007)}$ | $-3.649_{(0.017)}$ | $-3.396_{(0.008)}$ | $-3.047_{(0.026)}$ | $-3.355_{(0.034)}$

  * •

*Numbers in parentheses are standard errors for three independent runs. Grey boxes indicate the choice in the original paper.

Table G4: Predictive metrics (Time MAE and Type KLD) for GraDK with different
kernel ranks.

| 16-node graph | 50-node graph | Wildfire | Theft
---|---|---|---|---
Kernel rank | $R=1$ | $R=2$ | $R=3$ | $R=1$ | $R=2$ | $R=3$ | $R=1$ | $R=2$ | $R=3$ | $R=1$ | $R=2$ | $R=3$
$L=1$ | $0.272/0.003$ | $0.153/0.001$ | 0.028/<0.001 | $1.912/0.004$ | 1.023/0.004 | $3.118/0.006$ | $0.943/0.107$ | $0.358/0.024$ | 0.207/0.006 | $2.787/0.097$ | $1.052/0.082$ | 0.640/0.079
$L=2$ | $0.226/0.002$ | $0.145/0.001$ | $0.064/0.002$ | $3.698/0.006$ | $1.024/0.004$ | $5.217/0.013$ | $0.945/0.109$ | $0.314/0.020$ | $0.946/0.110$ | $2.702/0.091$ | $1.001/0.082$ | $1.123/0.109$

  * •

*Grey boxes indicate the choice in the original paper.

Figure G2: Learned $\\{\alpha_{rl}\\}$ of GraDK with $L=2,R=3$ on wildfire
data.

##  Appendix H Experimental details and additional results

In this section, we give details regarding the experiments in Section 4,
including a description of the ground truth point process models for synthetic
data, latent graph structures for real data, experimental setup, and
additional results.

###  H.1 Data description

Our experiments are implemented using five synthetic and four real-world data
sets.

##### Synthetic data.

We provide a detailed description of the ground truth kernels and latent graph
structures we used for synthetic data generation:

  * (i)

3 nodes with a non-stationary temporal kernel and 1-hop (positive and
negative) graph influence:

$k(t^{\prime},t,v^{\prime},v)=1.5(0.5+0.5\cos(0.2t^{\prime}))e^{-2(t-t^{\prime})}\left(0.5B_{1}^{(0)}(v^{\prime},v)+0.2B_{2}^{(1)}(v^{\prime},v)\right),$

where $B_{1}^{(0)}=\text{diag}(0.5,0.7,0.5)$, $B_{2}^{(1)}(2,1)=-0.2$, and
$B_{2}^{(1)}(2,3)=0.4$. The true graph kernel is visualized in the first panel
of Figure H5.

  * (ii)

16 nodes with a non-stationary temporal kernel and 2-hop graph influence:

$k(t^{\prime},t,v^{\prime},v)=1.5(0.5+0.5\cos(0.2t^{\prime}))e^{-2(t-t^{\prime})}\left(0.2B_{1}^{(0)}(v^{\prime},v)-0.3B_{2}^{(1)}(v^{\prime},v)+0.1B_{3}^{(2)}(v^{\prime},v)\right),$

such that
$(0.2B_{1}^{(0)}-0.3B_{2}^{(1)}+0.1B_{3}^{(2)})=(0.2I-0.3\tilde{L}+0.1(2\tilde{L}^{2}-I))$.
Here $\tilde{L}$ is the scaled and normalized graph Laplacian defined in
Appendix C. This graph influence kernel is visualized in the top row of Figure
1 (first column).

  * (iii)

50 nodes with an exponentially decaying temporal kernel:
$g(t^{\prime},t)=2e^{-2(t-t^{\prime})}$. The graph kernel is constructed such
that the influence follows a Gaussian density (with 3 modes) along the
diagonal of the graph influence kernel with random noise in addition to off-
ring influence. The true graph kernel is visualized in the first panel of
Figure H6.

  * (iv)

50 nodes with a non-separable diffusion graph kernel:
$k(t^{\prime},t,v^{\prime},v)=e^{-2(t-t^{\prime})}\cdot\frac{1}{8\pi(t-t^{\prime})}e^{-\frac{d^{2}(v,v^{\prime})}{8(t-t^{\prime})}}$.
The true graph kernel is visualized in the first panel of Figure H7.

  * (v)

225 nodes with exponentially decaying temporal kernel:
$g(t^{\prime},t)=2e^{-2(t-t^{\prime})}$. The true graph kernel is visualized
in the first panel of Figure H8.

The latent graph structures for these five synthetic data experiments can be
found in Figure H3. Data sets are simulated using the thinning algorithm [10].
Each data set contains 1,000 sequences with an average length of 50.9, 105.8,
386.8, 558.1, and 498.3, respectively.

(a) 3-node graph

(b) 16-node graph

(c) 50-node graph

(d) 225-node graph

Figure H3: Latent graph structures for the synthetic data sets. From left to
right, the graphs correspond to the ones in the 3-node graph (synthetic data
1), the 16-node ring graph (synthetic data 2), the 50-node graph (synthetic
data 3 & 4), and the 225-node graph (synthetic data 5), respectively.

(a) Traffic graph

(b) Sepsis graph

(c) Wildfire graph

(d) Theft graph

Figure H4: Latent graph structures for the real data sets. From left to right,
the graphs correspond to the ones in Atlanta traffic congestion data, the
Sepsis data, the California wildfire data, and the Valencia sustraccion
(theft) data.

##### Real data.

We also apply the model to four real data sets across different real-world
settings.

  * (i)

Traffic congestion data. The Georgia Department of Transportation (GDOT)
provides traffic volume data for sensors embedded on roads and highways in the
state. We have access to such data for 5 sensors at the interchange of
interstates 75 and 85 in Midtown Atlanta from September 2018 to March 2019.
Traffic volume is measured in 15-minute intervals, and congestion events are
detected when the traffic count exceeds the third quartile of the daily
traffic volume. The result is 3,830 events which are split into 24-hour
trajectories (with an average of 24 events per day).

  * (ii)

Sepsis data. The Sepsis data set released by PhysioNet Challenges [20, 37] is
a physiological data set that records patient covariates consisting of
demographics, vital signs, and laboratory values for ICU patients from three
separate hospital systems. Based on the covariates, common and clinically
relevant Sepsis-Associated Derangements (SADs) can be identified through
expert opinions. A SAD is considered present if the patient’s covariates are
outside of normal limits [42]. In our Sepsis data, we extract the onsets of 12
SADs and the Sepsis (total of 13 medical indices, as shown in Table H5) as
discrete events, with a total number of 80,463 events. Each sequence contains
the events extracted from one patient during a 24-hour period. The event times
are measured in hours, with the average length of each sequence being 15.6.

  * (iii)

Wildfire data. The California Public Utilities Commission (CPUC) maintains a
large-scale multi-modal wildfire incident dataset. We extract a total of 2,428
wildfire occurrences in California from 2014 to 2019. The latitude-longitude
coordinates of incidents are bounded by the rectangular region [34.51,
-123.50] $\times$ [40.73, -118.49]. Note that the majority of the region has
no fire in the 5-year horizon due to the fact that fire incidents are likely
to cluster in space. Therefore, we apply the K-means algorithm to extract 25
clusters of wildfire incidents. The entire dataset is split into one-year
sequences with an average length of 436 events.

  * (iv)

Theft data. The proprietary crime data collected by the police department in
Valencia, Spain records the crime incidents that happened from 2015 to 2019,
including incident location, crime category, and distance to various landmarks
within the city. We analyze 9,372 sustraccions (smooth thefts) that happened
near 52 pubs within the Valencia city area. Each sustraccion is assigned to
the closest pub. We partition the data into quarter-year-long sequences with
an average length of 469 events.

The learned latent graph structures overlaid on the real-world geography are
displayed in Figure H4(a)(c)(d) for the Atlanta traffic congestion, California
wildfire, and Valencia theft data. The learned structure for the Sepsis data
is in Figure H4(b)

Table H5: Identified 13 medical indices in Sepsis data.

Number | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13
---|---|---|---|---|---|---|---|---|---|---|---|---|---
Abbreviation | RI | EI | OCD | Shock | DCO | CI | Chl | HI | OD(Lab) | Inf(lab) | OD(vs) | Inf(vs) | Sepsis
Full name |  Renal Injury |  Electrolyte Imbalance |  Oxygen Carrying Dysfunction | Shock |  Diminished Cardiac Output | Coagulopathy | Cholestasis |  Hepatocellular Injury |  Oxygenation Dysfunction (Lab) | Inflammation (Lab) |  Oxygenation Dysfunction (vital sign) |  Inflammation (vital sign) | Sepsis

Table H6: Training hyper-parameters for the baselines.

|  3-node synthetic |  16-node synthetic |  50-node synthetic |  50-node synthetic with non-separable kernel |  225-node synthetic | Traffic | Sepsis | Wildfire | Theft
---|---|---|---|---|---|---|---|---|---
Model |  Learning Rate |  Batch Size |  Learning Rate |  Batch Size |  Learning Rate |  Batch Size |  Learning Rate |  Batch Size |  Learning Rate |  Batch Size |  Learning Rate |  Batch Size |  Learning Rate |  Batch Size |  Learning Rate |  Batch Size |  Learning Rate |  Batch Size
RMTPP | $10^{-3}$ | $32$ | $10^{-3}$ | $32$ | $10^{-3}$ | $32$ | $10^{-3}$ | $32$ | $10^{-3}$ | $32$ | $10^{-3}$ | $32$ | $10^{-3}$ | $100$ | $10^{-3}$ | $2$ | $10^{-3}$ | $2$
FullyNN | $10^{-2}$ | $100$ | $10^{-2}$ | $100$ | $10^{-2}$ | $100$ | $10^{-2}$ | $100$ | $10^{-2}$ | $100$ | $10^{-3}$ | $50$ | $10^{-3}$ | $100$ | $10^{-3}$ | $20$ | $10^{-3}$ | $100$
DNSK | $10^{-1}$ | $32$ | $10^{-1}$ | $32$ | $10^{-1}$ | $32$ | $10^{-1}$ | $32$ | $10^{-1}$ | $32$ | $10^{-1}$ | $32$ | $10^{-1}$ | $100$ | $10^{-1}$ | $2$ | $10^{-1}$ | $2$
THP-S | $10^{-3}$ | $32$ | $10^{-3}$ | $32$ | $10^{-3}$ | $32$ | $10^{-3}$ | $32$ | $10^{-3}$ | $32$ | $10^{-3}$ | $64$ | $10^{-3}$ | $100$ | $10^{-3}$ | $2$ | $10^{-3}$ | $2$
SAHP-G | $10^{-4}$ | $32$ | $10^{-4}$ | $32$ | $10^{-4}$ | $32$ | $10^{-4}$ | $32$ | $10^{-4}$ | $32$ | $10^{-4}$ | $32$ | $10^{-3}$ | $100$ | $10^{-3}$ | $2$ | $10^{-3}$ | $2$

###  H.2 Detailed experimental setup

We choose our temporal basis functions to be fully connected neural networks
with two hidden layers of width 32. Each layer, except for the output layer,
is equipped with the softplus activation function. The GraDK method sets the
kernel rank to be $L=1$ and $R=3$ for the final learned model for each data
set (except $L=1,R=2$ on the data generated by the 50-node graph, $L=2,R=4$ on
the data generated by the 50-node graph with the non-separable kernel, and
$L=2,R=4$ on the Sepsis data). For each data set, all the models are trained
using 80% of the data and tested on the remaining 20%. Our model parameters
are estimated through objective functions in Section 3.4 using the Adam
optimizer [25] with a learning rate of $10^{-2}$ and batch size of $32$.

For the baseline of RMTPP, we test the dimension of {32, 64, 128} for the
hidden embedding in RNN and choose an embedding dimension of 32 in the
experiments. For FullyNN, we set the embedding dimension to be 64 and use a
fully-connected neural network with two hidden layers of width 64 for the
cumulative hazard function, as the default ones in the original paper. The
dimension of input embedding is set to 10. For DNSK, we adopt the structure
for the marked point processes (for spatio-temporal point processes in the
ablation study) and set the rank of temporal basis and mark basis (spatial
basis in the ablation study) to one and three. THP-S and SAHP-G use the
default parameters from their respective code implementations. Learning rate
and batch size parameters are provided for each baseline and experiment in
Table H6.

We note that the original paper of FullyNN [33] does not consider the modeling
of event marks, and the model can only be used for temporal point processes,
for it requires the definition of cumulative intensity and the calculation of
the derivative, which is inapplicable in event mark space. To compare with
this method in the setting of graph point processes, we extend FullyNN to
model the mark distribution. Specifically, an additional feedforward fully-
connected neural network $\Psi$ is introduced. For each index $i$, it takes
the hidden states of history up to $i$-th event $\bm{h}_{i}$ as input and
outputs its prediction for the likelihood of the next mark $v_{i+1}$ as
$\Psi(v_{i+1}|\bm{h}_{i})$. Given a sequence of $n$ events
$\\{(t_{i},v_{i})\\}_{i=1}^{n}$, the optimization objective of model log-
likelihood for the extended FullyNN is written as:

$\displaystyle\mathcal{L}$
$\displaystyle=\text{log}L(\\{t_{i}\\})+\text{log}L(\\{v_{i}\\})$
$\displaystyle=\sum_{i}\left[\text{log}\left\\{\frac{\partial}{\partial\tau}\Phi(\tau=t_{i+1}-t_{i}|\bm{h}_{i})\right\\}-\Phi(t_{i+1}-t_{i}|\bm{h}_{i})\right]+\sum_{i}\text{log}\left\\{\Psi(v_{i+1}|\bm{h}_{i})\right\\},$

where $\Phi$ is the cumulative hazard function (i.e., cumulative intensity
function) introduced in FullyNN, parametrized by a feedforward neural network.
Such an approach to modeling the distribution of event marks has been commonly
adopted in previous studies [14, 39].

###  H.3 Additional experimental results

Figure H5: Graph kernel, inter-event dependence, and conditional intensity
recovery for the synthetic data set generated by the 3-node graph with
negative (inhibitive) graph influence. The first column reflects the ground
truth, while the subsequent columns reflect the results obtained by GraDK (our
method), SAHP-G, and DNSK, respectively. Figure H6: Graph kernel, inter-event
dependence, and conditional intensity recovery for the synthetic data set
generated by the 50-node graph. Figure H7: Graph kernel, inter-event
dependence, and conditional intensity recovery for the synthetic data set
generated by the 50-node graph with a non-separable kernel.

##### Synthetic data.

Figure H5 presents the kernel and intensity recovery results for the synthetic
data set 1 generated by a non-stationary temporal kernel and a 3-node-graph
kernel with inhibitive influence. The recovery outcomes demonstrate the
efficacy of our proposed model in characterizing the temporal and graph-based
dependencies among events, as evident from the first and second rows of the
figure. Moreover, the results emphasize the capability of our model, GraDK, to
effectively capture the inhibitive effects of past events. In contrast, the
event dependencies represented by the normalized positive attention weights in
SAHP solely capture the triggering intensity of past events without accounting
for inhibitive influences. Similarly, the first row of Figure H6 displays the
true graph influence kernel in the 50-node synthetic data set and the learned
graph kernels by GraDK, SAHP-G, and DNSK. While SAHP-G exaggerates the self-
exciting influence of graph nodes and DNSK only learns some semblance of the
graph kernel behavior, our method accurately recovers both self- and inter-
node influence patterns, resulting in a faithful model of the true graph
kernel. The conditional intensity via each method for one testing trajectory
is displayed in the third row of Figure H6, which demonstrates the capability
of our model to capture the temporal dynamics of events. The same results are
presented in Figure H7 and Figure H8 for the synthetic data sets generated by
the 50-node graph with non-separable kernel and 225-node graph, respectively.

Figure H8: Graph kernel, inter-event dependence, and conditional intensity
recovery for the synthetic data set generated by the 225-node graph.

(a) Traffic data

(b) Sepsis data

(c) Wildfire data

Figure H9: Learned graph kernels for traffic, Sepsis, and wildfire data set
evaluated at $t^{\prime}=0$ and $t-t^{\prime}=1$. The columns present the
recovered kernels on each data set of GraDK+L3net, GraDK+GAT, SAHP-G, and
DNSK, respectively.

##### Real data.

We visualize the learned graph kernels by GraDK, SAHP, and DNSK on traffic,
Sepsis, and wildfire data in Figure H9. Our approach is able to learn intense
graph signals with complex patterns by taking advantage of the graph structure
and GNNs. While the attention structures adopted in SAHP contribute to
improved model prediction performance for future events, this approach suffers
from the limited model expressiveness and interpretability when attempting to
recover the underlying mechanism of the event generation process, indicated by
the weak and noisy graph signals. DNSK fails to uncover the intricate patterns
existing in graph kernels and only provides restrictive kernels for event
modeling without considering the latent graph structures among data.

(a) Traffic network structure

(b) Northbound sites

(c) Southbound sites

Figure H10: (a) Traffic network structure. Each traffic sensor is labeled with
a number. The blue and pink arrows indicate the monitored traffic directions
of sensors on northbound and southbound highways, respectively. (b)(c)
Conditional intensity of five sensors in a single day, which are categorized
into two groups according to the monitored traffic directions of the sensors.

Our model not only achieves exceptional interpretability but also holds
practical significance in the context of real-world point process data
modeling. This can be demonstrated through experimental results conducted on
traffic data. The five traffic sensors from which we collect data can be
categorized into two groups, northbound sites and southbound sites, according
to the directions of the highway they are monitoring. Figure H10(a) visualizes
the structure of the traffic network with sensors (graph nodes) labeled and
arrows indicating the highway directions. We then investigate the conditional
intensity learned by GraDK on each traffic sensor. Figure H10(b)(c) show the
conditional intensity functions of five sensors during one single day (i.e.,
computed given one sequence from the testing set) with two subgroups of
northbound and southbound sites. Note that similar temporal patterns can be
found within the same subgroup, which can be attributed to the fact that the
sensors in the same group are in the same direction and share the same traffic
flow. Also, all the intensity functions show a temporal pattern in which they
reach pinnacles during the morning (around 8:00) and evening (around 17:00)
rush hours, with a higher possibility for traffic congestion at southbound
sensors in the afternoon.