# Hawkes Processes on Graphons

Hongteng Xu  Dixin Luo  Hongyuan Zha

###### Abstract

We propose a novel framework for modeling multiple multivariate point
processes, each with heterogeneous event types that share an underlying space
and obey the same generative mechanism. Focusing on Hawkes processes and their
variants that are associated with Granger causality graphs, our model
leverages an uncountable event type space and samples the graphs with
different sizes from a nonparametric model called graphon. Given those graphs,
we can generate the corresponding Hawkes processes and simulate event
sequences. Learning this graphon-based Hawkes process model helps to 1) infer
the underlying relations shared by different Hawkes processes; and 2) simulate
event sequences with different event types but similar dynamics. We learn the
proposed model by minimizing the hierarchical optimal transport distance
between the generated event sequences and the observed ones, leading to a
novel reward-augmented maximum likelihood estimation method. We analyze the
properties of our model in-depth and demonstrate its rationality and
effectiveness in both theory and experiments.

Hawkes process, graphon, hierarchical optimal transport, heterogeneous event
sequences

##  1 Introduction

As a powerful statistical tool, Hawkes process (Hawkes, 1971) has been widely
used to model event sequences in the continuous-time domain. Suppose that we
have an event sequence
$\\{(t_{i},v_{i})\in[0,T]\times\mathcal{V}\\}_{i=1}^{N}$, where $[0,T]$ is the
observation time window, $\mathcal{V}$ is the set of event types, and
$(t_{i},v_{i})$ is the $i$-th event at time $t_{i}$ with type $v_{i}$.
Equivalently, we can represent the sequence by a counting process
$\bm{N}(t)=\\{N_{v}(t)\\}_{v\in\mathcal{V}}$, where $N_{v}(t)$ is the number
of the type-$v$ events till time $t$. A Hawkes process characterizes the
expected instantaneous rate of occurrence of the type-$v$ event at time $t$ by
a conditional intensity function (Liniger, 2009): for $v\in\mathcal{V}$ and
$t\in[0,T]$,

$\displaystyle\begin{aligned}
\lambda_{v}(t):=\frac{\mathbb{E}[\text{d}N_{v}(t)|\mathcal{H}_{t}]}{\text{d}t}=\mu_{v}+\sideset{}{{}_{t_{i}<t}}{\sum}\phi_{vv_{i}}(t,t_{i}).\end{aligned}$
(1)

Here, $\mathcal{H}_{t}=\\{(t_{i},v_{i})|t_{i}<t\\}$ contains the past events
till time $t$. $\mu_{v}\geq 0$ is the base rate of type-$v$ event.
$\\{\phi_{vv^{\prime}}(t,t^{\prime})\geq
0\\}_{v,v^{\prime}\in\mathcal{V},t^{\prime}<t}$ are the so called impact
functions, and $\phi_{vv^{\prime}}(t,t^{\prime})$ quantifies the influence of
the type-$v^{\prime}$ event at time $t^{\prime}$ on the type-$v$ event at time
$t$. Accordingly, $\sum_{t_{i}<t}\phi_{vv_{i}}(t,t_{i})$ accumulates the
impacts of the past events. The set of impact functions gives rise to the
Granger causality graph of the event types (Eichler et al., 2017; Xu et al.,
2016a), denoted as $G(\mathcal{V},\mathcal{E})$ — an edge
$v^{\prime}\rightarrow v\in\mathcal{E}$ means that a past type-$v^{\prime}$
event can trigger the occurrence of a type-$v$ event in the future, and
$v^{\prime}\rightarrow v\notin\mathcal{E}$ if and only if
$\phi_{vv^{\prime}}(t,t^{\prime})\equiv 0$.

Hawkes process, together with the corresponding Granger causality graph of
event types, has become instrumental for many applications involving event
sequences, such as social network modeling (Farajtabar et al., 2017) and
financial data analysis (Bacry et al., 2015). Interestingly, even with recent
models enhancing Hawkes processes with deep neural networks (Mei & Eisner,
2017; Zhang et al., 2020; Zuo et al., 2020), the work in (Tank et al., 2018)
shows that the group sparsity of their neural networks’ parameters can still
be interpreted by Granger causality of the event types.

Figure 1: An illustration of the Hawkes processes on a graphon.

Despite achieving many successes, the applications of the Hawkes-related
processes are limited for homogeneous scenarios in which all the event
sequences are generated by one point process defined on a known set of event
types. Although some methods consider learning multiple point processes for
the sequences in different clusters (Luo et al., 2015; Xu & Zha, 2017) or time
periods (Lin et al., 2016; Alaa et al., 2017), they still maintain a single
set of event types. This setting, however, is in conflict with the
heterogeneous nature of many real-world event sequences — the event types are
often sampled from an underlying event type space, and new sequences are
driven by the latent sampling process and are generated with event types
unobserved before. Accordingly, for different event sequences, their point
processes are defined with different event types, and thus, obey different
generative mechanisms. We illustrate this phenomenon via event sequences
defined on networks.

Social networks. Users of different networks, $e.g.$, Facebook and Twitter,
are actually sampled from the same underlying populations ($i.e.$, all the
Internet users in the world). When using Hawkes processes to model the user
behaviors on those different networks (Blundell et al., 2012; Zhou et al.,
2013; Zhao et al., 2015), the respective users are considered as event types
and their corresponding Granger causality graphs can be treated as different
subgraphs sampled from a large latent graph for the whole population.
Additionally, with the entering of new users and the exiting of old ones,
those networks are time-varying and their corresponding Hawkes processes at
different time are different.

Patient admissions. For a patient suffering from several diseases, his
admissions in a hospital over time are often assumed to be driven by the
Granger causality graph of his diseases ($i.e.$, disease graph), and thus,
modeled by a Hawkes process (Xu et al., 2017). For patients with different
profiles, even for the same patient in different age periods, their disease
graphs and the corresponding Hawkes processes can be very different. The
diseases appearing in each Hawkes process are sampled from the same set of
diseases, $e.g.$, the international classification of diseases (ICD), and each
specific disease graph is a subgraph of an unknown graph constructed by all
the diseases in the set. Moreover, with the development of biomedical science,
we may find new diseases and observe new admissions in the future.

Besides these two typical examples, the sequential shopping behaviors on
different platforms, the transactions of stocks in different markets, and the
diffusion of a virus in different cities, $etc.$, all these event sequences
are heterogeneous, whose event types can better be modeled as samples from an
underlying infinite even uncountable event type space. When modeling such
event sequences, we need to learn a generative model for their point processes
beyond just learning a single point process for each of them individually.

To this end, we propose a new graphon-based Hawkes process (GHP). Essentially,
our GHP is a hierarchical generative model for a collection of Hawkes
processes with heterogeneous types (and their variants). As illustrated in
Figure 1, it not only models the generative mechanisms of event sequences by
Hawkes processes but also designs a graphon model (Lovász, 2012) to generate
the event types of the different Hawkes processes from an uncountable event
type space. By sampling the graphon, we generate the parameters of various
Hawkes processes and simulate event sequences accordingly. Unlike existing
Hawkes-related processes, our GHP model is able to generate different Hawkes
processes with heterogeneous event types but similar dynamics. For more
complicated point processes, we can extend our GHP model by leveraging neural
networks and applying multi-dimensional graphons.

Our GHP model is theoretically grounded: with mild assumptions, we demonstrate
that for the generated Hawkes processes, the proposed model $i$) guarantees
their stationarity; $ii$) ensures their parameters to be Lipschitz continuous;
and $iii$) makes the difference between their corresponding event sequences
bounded. These properties guarantee the stability of our GHP model when
generating Hawkes processes and their event sequences.

Learning GHP from observed heterogeneous event sequences requires us to infer
and align the corresponding Hawkes processes with respect to the underlying
graphon, for which traditional methods like maximum likelihood estimation are
infeasible. To overcome this problem, we design a novel learning algorithm
based on the reward-augmented maximum likelihood (RAML) estimation (Norouzi et
al., 2016) and the hierarchical optimal transport (HOT) distance (Lee et al.,
2019; Yurochkin et al., 2019). In particular, given observed event sequences
and those generated by our GHP model, we calculate the HOT distance between
them and obtain an optimal transport matrix corresponding to their joint
probabilities. The probabilities work as the rewards modulating the log-
likelihood of each generated event sequence. Taking the reward-augment log-
likelihood as an objective, we estimate the parameters of GHP accordingly. We
verify the feasibility of our GHP model and its learning algorithm on both
synthetic and real-world data. When modeling sparse heterogeneous event
sequences that have many event types but small number of events, our GHP model
significantly mitigates the risk of over-fitting and thus outperforms other
state-of-the-art point process models.

##  2 Graphon-based Hawkes Processes

###  2.1 Generating Hawkes processes from a graphon

For a classic Hawkes process, we often parameterize its impact functions as
$\\{\phi_{vv^{\prime}}(t,t^{\prime})=a_{vv^{\prime}}\eta(t-t^{\prime})\\}_{v,v^{\prime}\in\mathcal{V}}$,
where the coefficient $a_{vv^{\prime}}\geq 0$ and the decay kernel
$\eta(t)\geq 0$. The decay kernel is predefined, and its integral is
$D=\int_{0}^{\infty}\eta(t)\text{d}t$. Such a Hawkes process is denoted as
$\text{HP}_{\mathcal{V}}(\bm{\mu},\bm{A})$, where $\mathcal{V}$ is the set of
event types, $\bm{\mu}=[\mu_{v}]\in\mathbb{R}^{|\mathcal{V}|}$ and
$\bm{A}=[a_{vv^{\prime}}]\in\mathbb{R}^{|\mathcal{V}|\times|\mathcal{V}|}$.
Here, $|\mathcal{V}|$ is the cardinality of $\mathcal{V}$. For
$\text{HP}_{\mathcal{V}}(\bm{\mu},\bm{A})$, $\bm{A}$ is the adjacency matrix
of the corresponding Granger causality graph.

A potential way to generate Hawkes processes is to first simulate their
Granger causality graphs. We apply this strategy based on a nonparametric
graph model called graphon (Lovász, 2012). A graphon is a two-dimensional
measurable function, denoted as $g:\Omega^{2}\mapsto[0,1]$, where $\Omega$ is
a measure space. Given a graphon, we can sample a matrix
$\bm{A}=[a_{vv^{\prime}}]\in[0,1]^{V\times V}$ with an arbitrary size $V$:

$\displaystyle\begin{aligned}
a_{vv^{\prime}}=g(x_{v},x_{v^{\prime}}),~{}x_{v}\sim\text{Uniform}(\Omega)~{}\text{for}~{}v=1,..,V.\end{aligned}$
(2)

Here, $\\{x_{v}\in\Omega\\}_{v=1}^{V}$ are $V$ independent variables sampled
from a uniform distribution. Accordingly, we generate a graph
$G(\mathcal{V},\mathcal{E})$ by setting $\mathcal{V}=\\{1,..,V\\}$ and
$v^{\prime}\rightarrow v\in\mathcal{E}\sim\text{Bernoulli}(a_{vv^{\prime}})$.
This graphon model is fundamental for modeling large-scale networks, which has
been widely used in network analysis (Gao & Caines, 2019).

Besides $g(x,y)$, we introduce a one-dimensional measurable function on
$\Omega$, $i.e.$, $f:\Omega\mapsto[0,+\infty)$, such that we can sample
$\bm{\mu}$ and $\bm{A}$ of a Hawkes process from $f(x)$ and $g(x,y)$,
respectively. Our graphon-based Hawkes process model consists of $f(x)$ and
$g(x,y)$, denoted as $\text{GHP}_{\Omega}(f,g)$. Here, we set $\Omega=[0,1]$
and implement the functions as

$\displaystyle\begin{aligned}
f(x)&=\text{softplus}(f_{1})(\exp(\sigma(f_{2})x)-1),\\\
g(x,y)&=\sigma\Bigl{(}\sideset{}{{}_{i,j\in\\{0,..,S\\}}}{\sum}(g_{ij}^{1}\sin
i\pi x+g_{ij}^{2}\cos i\pi x)\\\ &\hskip 62.59596pt\times(g_{ij}^{3}\sin j\pi
y+g_{ij}^{4}\cos j\pi y)\Bigr{)},\end{aligned}$ (3)

where $f(x)$ is an exponential function, $g(x,y)$ is designed based on the 2D
Fourier series, which has $4(S+1)^{2}$ coefficients, and $\sigma(\cdot)$ is
the sigmoid function. This implementation is simple and makes our model
satisfy some significant properties in theory, which will be shown in Section
2.3. Then the generative process defined by $\text{GHP}_{\Omega}(f,g)$ is

$\displaystyle\begin{aligned}
&\text{HP}_{\mathcal{V}}(\bm{\mu},\bm{A})\sim\text{GHP}_{\Omega}(f,g):\\\
&\quad 1)~{}V\sim\bm{\pi}=\\{\pi_{1},...,\pi_{V_{\max}}\\},\\\ &\quad
2)~{}\mathcal{V}=\\{1,..,V\\},~{}\text{and}~{}x_{v}\sim\text{Uniform}(\Omega),~{}\forall
v\in\mathcal{V}.\\\ &\quad
3)~{}\mu_{v}=f(x_{v}),~{}a_{vv^{\prime}}=\frac{1}{V_{\max}D}g(x_{v},x_{v^{\prime}}).\\\
&\bm{N}(t)\sim\text{HP}_{\mathcal{V}}(\bm{\mu},\bm{A}).\end{aligned}$ (4)

Here, $\bm{\pi}$ is a categorical distribution on $\\{1,...,V_{\max}\\}$,
which is often set as a uniform distribution, and $V_{\max}$ is the maximum
number of event types supported by our model. We treat $\Omega$ as an
uncountable event type space. In each trial, we sample $V$ latent event types
$\\{x_{v}\\}_{v=1}^{V}$ from $\Omega$, where the number of the event types $V$
is sampled from $\bm{\pi}$. Based on $\\{x_{v}\\}_{v=1}^{V}$, we sample
$\bm{\mu}$ and $\bm{A}$ from $f$ and $g$, respectively, and instantiate a
Hawkes process. Different from (2), we set
$a_{vv^{\prime}}=\frac{1}{V_{\max}D}g(x_{v},x_{v^{\prime}})$ in (4) to ensure
the Hawkes process is stationary.

######  Property 2.1 (Stationarity).

$\text{HP}_{\mathcal{V}}(\bm{\mu},\bm{A})\sim\text{GHP}_{\Omega}(f,g)$ is
asymptotically stationary as long as $|\mathcal{V}|\leq V_{\max}$.

Therefore, we can readily generate an event sequence $\bm{N}(t)$ from
$\text{HP}_{\mathcal{V}}(\bm{\mu},\bm{A})$ by various simulation methods,
$e.g.$, the branch processing (Møller & Rasmussen, 2006) and Ogata’s thinning
method (Ogata, 1981).

The key challenge in using GHP is that we cannot observe
$\\{x_{v}\\}_{v=1}^{V}$ because both the event type space $\Omega$ and the
sampled event types are latent. Accordingly, for the generated Hawkes
processes and their event sequences, we cannot directly match their event
types ($i.e.$, $\\{x_{v}\\}_{v=1}^{V}$) with the event types of real-world
sequences. To solve this problem, in Section 3.2 we will leverage optimal
transport (Villani, 2008; Peyré et al., 2019) to measure the distance between
heterogeneous event sequences. The learned optimal transport helps us to find
a soft alignment between the generated event types and the real ones, which
not only makes the generated event types and the corresponding point processes
semantically meaningful but also builds the foundation for the learning method
of our model (See Section 3).

###  2.2 Extensions

The proposed GHP provides us with a new framework to jointly model
heterogeneous event sequences. Beyond Hawkes processes, our GHP model can be
readily extended to generate more sophisticated types of point processes.

Nonlinear Hawkes process. For nonlinear Hawkes process (also called mutually-
correcting process) (Zhu, 2013; Xu et al., 2016b), its intensity function is
$\lambda_{v}(t)=\exp(\mu_{v}+\sum_{t_{i}<t}\phi_{vv_{i}}(t,t_{i}))$ and the
parameters can be negative. In this case, we can implement
$\text{GHP}_{\Omega}(f,g)$ with $f:\Omega\mapsto(-\infty,+\infty)$ and
$g:\Omega^{2}\mapsto(-\infty,+\infty)$, respectively.

Multi-kernel Hawkes process. The multi-kernel Hawkes process constructs its
impact functions by a set of decay kernels (Xu et al., 2016a), $i.e.$,
$\phi_{vv^{\prime}}(t)=\sum_{m=1}^{M}a_{vv^{\prime}m}\eta_{m}(t)$, where the
coefficients $a_{vv^{\prime}m}$’s are formulated as $M$ matrices
$\\{\bm{A}_{m}\\}_{m=1}^{M}$. In this case, we need to introduce several
graphons, $i.e.$, $\\{g_{1}(x,y),...,g_{M}(x,y)\\}$, to generate the $M$
matrices, and our GHP model becomes
$\text{GHP}_{\Omega}(f,\\{g_{m}\\}_{m=1}^{M})$.

Time-varying Hawkes process. The time-varying Hawkes process applies shift-
varying impact functions, $i.e.$,
$\phi_{vv^{\prime}}(t,t^{\prime})=a_{vv^{\prime}}(t)\eta(t-t^{\prime})$, where
the coefficient $a_{vv^{\prime}}(t)$ becomes a function of time. Similar to
the multi-kernel Hawkes process, when using a set of bases to represent the
coefficient function (Xu et al., 2017), $i.e.$,
$a_{vv^{\prime}}(t)=\sum_{m=1}^{M}a_{vv^{\prime}m}h_{m}(t)$, where $h_{m}(t)$
is the $m$-th base, we can still apply multiple graphons to generate impact
functions and rewrite our GHP model as
$\text{GHP}_{\Omega}(f,\\{g_{m}\\}_{m=1}^{M})$.

Neural Hawkes process Most existing neural network-based Hawkes processes
apply embedding layers to map the index of each event type to its latent code
(Mei & Eisner, 2017; Zhang et al., 2020; Zuo et al., 2020). For the neural
Hawkes process, we can replace the embedding layer with a function
$f(x):\Omega\mapsto\mathbb{R}^{M}$ such that we can generate $M$-dimensional
latent codes for uncountable event types in $\Omega$. If the neural Hawkes
process considers the interactions of different event types (Wang et al.,
2016), we can set the graphon as $g(x,y)=p(x)^{\top}q(y)$, where
$p(x):\Omega\mapsto\mathbb{R}^{M}$ and $q(y):\Omega\mapsto\mathbb{R}^{M}$,
respectively. Accordingly, the GHP becomes $\text{GHP}_{\Omega}(f,p,q)$.
Besides changing the point process model, we can also implement $f(x)$ and
$g(x,y)$ by deep neural networks, which is left for future work.

###  2.3 Theoretical analysis of the GHP model

In addition to verifying the stationarity of generated Hawkes processes, we
demonstrate two more properties of GHP based on the following mild
assumptions.

######  Assumption 2.2.

For $\text{GHP}_{\Omega}(f,g)$, we assume

  * A)

$f(x)$ is bi-Lipschitz continuous on $\Omega$, denoted as
$f\in\text{Lip}_{\Omega}(C_{1}^{f},C_{2}^{f})$: $\exists~{}0<C_{1}^{f}\leq
C_{2}^{f}<\infty$, $C_{1}^{f}\|x-x^{\prime}\|_{2}\leq|f(x)-f(x^{\prime})|\leq
C_{2}^{f}\|x-x^{\prime}\|_{2}$, $\forall x,x^{\prime}\in\Omega$.

  * B)

$f(x)$ has a unique zero point in $\Omega$, $i.e.$, $f(x_{0}^{f})=0$.

  * C)

$g(x,y)$ is strictly smaller than $1$, $i.e.$, $g:\Omega^{2}\mapsto[0,1)$.

  * D)

$g(x,y)$ is Lipschitz continuous on $\Omega^{2}$, denoted as
$g\in\text{Lip}_{\Omega^{2}}(C^{g})$: $\exists~{}0<C^{g}<\infty$,
$|g(x,y)-g(x^{\prime},y^{\prime})|\leq
C^{g}\|[x;y]-[x^{\prime};y^{\prime}]\|_{2}$,
$\forall[x;y],[x^{\prime};y^{\prime}]\in\Omega^{2}$.

Clearly, GHP defined in (3) satisfies the assumptions. Based on the above
assumptions, we prove that the parameters of the Hawkes process generated by
our GHP model is Lipschitz continuous.

######  Property 2.3 (Lipschitz Continuity).

For $\text{HP}_{\mathcal{V}}(\bm{\mu}_{1},\bm{A}_{1})$ and
$\text{HP}_{\mathcal{U}}(\bm{\mu}_{2},\bm{A}_{2})\sim\text{GHP}_{\Omega}(f,g)$,
where $\text{GHP}_{\Omega}(f,g)$ satisfies Assumption 2.2, their parameters
satisfy

$\displaystyle\begin{aligned}
&C_{1}^{f}d_{\text{w}}(\bm{x}_{1},\bm{x}_{2})\leq
d_{\text{w}}(\bm{\mu}_{1},\bm{\mu}_{2})\leq
C_{2}^{f}d_{\text{w}}(\bm{x}_{1},\bm{x}_{2}),\\\
&d_{\text{w}}(\bm{A}_{1},\bm{A}_{2})\leq
C^{g}d_{\text{w}}(\bm{x}_{1}^{\times},\bm{x}_{2}^{\times}),\\\
&d_{\text{gw}}(\bm{A}_{1},\bm{A}_{2})\leq
C^{g}d_{\text{gw}}(\bm{x}_{1}^{\times},\bm{x}_{2}^{\times}),\end{aligned}$ (5)

where $\bm{x}_{1}=\\{x_{v,1}\\}_{v=1}^{|\mathcal{V}|}$ and
$\bm{x}_{2}=\\{x_{u,2}\\}_{u=1}^{|\mathcal{U}|}$ are the latent event types,
and
$\bm{x}_{1}^{\times}=\\{[x_{v,1};x_{v^{\prime},1}]\\}_{v,v^{\prime}=1}^{|\mathcal{V}|}$
and
$\bm{x}_{2}^{\times}=\\{[x_{u,2};x_{u^{\prime},2}]\\}_{u,u^{\prime}=1}^{|\mathcal{U}|}$
enumerate the pairs of the latent event types. $d_{\text{w}}$ is the discrete
Wasserstein distance (or called the earth mover’s distance) and the
$d_{\text{gw}}$ is the discrete Gromov-Wasserstein distance.111The definitions
of $d_{\text{w}}$ and $d_{\text{gw}}$ are given in Appendix A.

Property 2.3 shows that $i$) for the generated Hawkes processes, the
difference between their parameters is bounded by the difference between their
latent event types; and $ii$) the parameters of each generated Hawkes process
are robust to the perturbations of the latent event types.

Because the difference between generated Hawkes processes is bounded, the
difference between the corresponding event sequences is bounded as well.
Specifically, for a point process, its average intensity vector, defined as
$\bar{\bm{\lambda}}:=\frac{\mathbb{E}[\text{d}\bm{N}(t)]}{\text{d}t}$,
reflects the dynamics of its event sequences (Chiu et al., 2013). For this key
statistics, we have

######  Property 2.4.

For $\text{HP}_{\mathcal{V}}(\bm{\mu}_{1},\bm{A}_{1})$ and
$\text{HP}_{\mathcal{U}}(\bm{\mu}_{2},\bm{A}_{2})\sim\text{GHP}_{\Omega}(f,g)$,
where $\text{GHP}_{\Omega}(f,g)$ satisfies Assumption 2.2 and
$|\mathcal{V}|\leq|\mathcal{U}|$, their average intensity vectors, $i.e.$,
$\bar{\bm{\lambda}}_{1}$ and $\bar{\bm{\lambda}}_{2}$, satisfy

$\displaystyle\begin{aligned}
\frac{d_{\text{w}}(\bar{\bm{\lambda}}_{1},\bar{\bm{\lambda}}_{2})}{\|\bar{\bm{\lambda}}_{1}\|_{2}}\leq&\frac{\frac{\sqrt{2U}C^{g}}{C_{1}^{f}\|\bm{I}_{V}-D\bm{A}_{1}\|_{2}}+\frac{1}{\|\bm{\mu}_{1}\|_{2}}}{1-D\|\bm{A}_{1}\|_{2}}\Bigl{(}d_{\text{w}}(\bm{\mu}_{1},\bm{\mu}_{2})\\\
&+\sqrt{\frac{U-V}{V}}\|\bm{\mu}_{1}\|_{2}\Bigr{)}+\sqrt{\frac{U-V}{UV}},\end{aligned}$
(6)

where $\|\cdot\|_{2}$ is the $\ell_{2}$-norm for vectors and the spectral norm
for matrices, $U=|\mathcal{U}|$, $V=|\mathcal{V}|$,
$D=\int_{0}^{\infty}\eta(t)\text{d}t$ is the integral of the decay kernel used
in the Hawkes processes, and $C_{1}^{f}$ and $C^{g}$ are the constants defined
in Assumption 2.2.

Furthermore, if $|\mathcal{V}|=|\mathcal{U}|$, we can simplify Property 2.4 as

######  Corollary 2.5.

For $\text{HP}_{\mathcal{V}}(\bm{\mu}_{1},\bm{A}_{1})$ and
$\text{HP}_{\mathcal{U}}(\bm{\mu}_{2},\bm{A}_{2})\sim\text{GHP}_{\Omega}(f,g)$,
where $\text{GHP}_{\Omega}(f,g)$ satisfies Assumption 2.2 and
$|\mathcal{V}|=|\mathcal{U}|=V$, we have

$\displaystyle\begin{aligned}
\frac{d_{\text{w}}(\bar{\bm{\lambda}}_{1},\bar{\bm{\lambda}}_{2})}{\|\bar{\bm{\lambda}}_{1}\|_{2}}\leq\frac{d_{\text{w}}(\bm{\mu}_{1},\bm{\mu}_{2})}{1-D\|\bm{A}_{1}\|_{2}}\left(\frac{\sqrt{2V}{C^{g}}/{C_{1}^{f}}}{\|\bm{I}_{V}-D\bm{A}_{1}\|_{2}}+\frac{1}{\|\bm{\mu}_{1}\|_{2}}\right).\end{aligned}$

##  3 Learning Algorithm

###  3.1 A reward-augmented maximum likelihood

We propose a novel method to learn GHP model from observed heterogeneous event
sequences. Denote $\mathcal{N}=\\{\bm{N}_{l}(t)\\}_{l=1}^{L}$ as the set of
real-world event sequences and
$\widehat{\mathcal{N}}=\\{\widehat{\bm{N}}_{k}(t)\\}_{k=1}^{K}$ the set of the
event sequences generated by our model. Because the correspondence of real-
world event types in the latent event type space is unknown, as mentioned in
Section 2.1, we need to simultaneously learn the underlying graphon of our
model and align the event types of the generated Hawkes processes with the
real ones. To achieve this aim, we formulate the following optimization
problem,

$\displaystyle\begin{aligned}
\sideset{}{{}_{\theta}}{\min}-\sideset{}{{}_{\widehat{\bm{N}}_{k}\in\widehat{\mathcal{N}}}}{\sum}\sideset{}{{}_{\bm{N}_{l}\in\mathcal{N}}}{\max}q(\widehat{\bm{N}}_{k}|\bm{N}_{l})\log
p(\widehat{\bm{N}}_{k};\theta).\end{aligned}$ (7)

where
$p(\widehat{\bm{N}}_{k};\theta)=\frac{\prod_{(t_{i},v_{i})}\lambda_{v_{i}}^{k}(t_{i};\theta)}{\exp(\sum_{v\in\mathcal{V}}\int_{0}^{T}\lambda_{v}^{k}(t;\theta)\text{d}t)}$
is the likelihood of the $k$-th generated event sequence, $\theta$ represents
the model parameter $\\{f_{1},f_{2},\\{g_{ij}^{m}\\}\\}$, and
$q(\widehat{\bm{N}}_{k}|\bm{N}_{l})$ is the probability of
$\widehat{\bm{N}}_{k}$ conditioned on the $l$-th real sequence $\bm{N}_{l}$.
Essentially, the conditional probability $q(\widehat{\bm{N}}|\bm{N})$ measures
the similarity between the generated sequence and the real one. When the two
sequences yield the same generative mechanism and have similar dynamics, the
real sequence provides useful prior information, and thus, the occurrence of
the generated sequence is with a high probability.

In (7), the log-likelihood of each generated sequence is weighted by its
maximum conditional probability with respect to the real sequences, $i.e.$,
$\max_{\bm{N}_{l}\in\mathcal{N}}q(\widehat{\bm{N}}_{k}|\bm{N}_{l}))$. The
weight measures the overall similarity between the the generated sequence
$\widehat{\bm{N}}_{k}$ and the real ones. A large weight indicates that the
generated sequence is informative for our learning problem because it is
similar to at least one real sequence. Otherwise, the sequence is less useful.
Additionally, assuming the empirical distribution of the real sequences to be
uniform, we have $q(\widehat{\bm{N}}|\bm{N})\propto
q(\widehat{\bm{N}},\bm{N})$, and the optimization problem becomes,

$\displaystyle\begin{aligned}
\sideset{}{{}_{\theta}}{\min}-\sideset{}{{}_{\widehat{\bm{N}}_{k}\in\widehat{\mathcal{N}}}}{\sum}\sideset{}{{}_{\bm{N}_{l}\in\mathcal{N}}}{\max}q(\widehat{\bm{N}}_{k},\bm{N}_{l})\log
p(\widehat{\bm{N}}_{k};\theta).\end{aligned}$ (8)

Remark. the above formulation (8) can be considered as a variant of the
reward-augmented maximum likelihood (RAML) estimation method (Norouzi et al.,
2016) (see also Section 3.3). For sequence $\widehat{\bm{N}}_{k}$, the weight
$\max_{\bm{N}_{l}\in\mathcal{N}}q(\widehat{\bm{N}}_{k},\bm{N}_{l})$ plays the
role of its reward and is assigned to its log-likelihood. The higher reward
the log-likelihood obtains, the more significant it is in learning.

###  3.2 Hierarchical optimal transport between heterogeneous event sequences

The key of our learning algorithm, which is also its main novelty, is
computing the joint distribution $q(\widehat{\bm{N}},\bm{N})$ based on the
hierarchical optimal transport (HOT) model (Lee et al., 2019; Yurochkin et
al., 2019). In particular, the HOT model not only captures the optimal
transport between the generated event sequences and the real ones but also
captures the optimal transport between their event types. Given
$\widehat{\mathcal{N}}=\\{\widehat{\bm{N}}_{k}\\}_{k=1}^{K}$ and
$\mathcal{N}=\\{\bm{N}_{l}\\}_{l=1}^{L}$, we compute the optimal transport
distance between them as $d_{\text{ot}}(\widehat{\mathcal{N}},\mathcal{N})$

$\displaystyle\begin{aligned}
:=&\sideset{}{{}_{\bm{Q}\in\Pi\left(\frac{1}{K}\bm{1}_{K},\frac{1}{L}\bm{1}_{L}\right)}}{\min}\sideset{}{{}_{k,l}}{\sum}q(\widehat{\bm{N}}_{k},\bm{N}_{l})d(\widehat{\bm{N}}_{k},\bm{N}_{l})\\\
=&\sideset{}{{}_{\bm{Q}\in\Pi\left(\frac{1}{K}\bm{1}_{K},\frac{1}{L}\bm{1}_{L}\right)}}{\min}\langle\bm{D},\bm{Q}\rangle,\end{aligned}$
(9)

where the polytope
$\Pi\left(\frac{1}{K}\bm{1}_{K},\frac{1}{L}\bm{1}_{L}\right)=\\{\bm{Q}\geq\bm{0}~{}|~{}\bm{Q}\bm{1}_{L}=\frac{1}{K}\bm{1}_{K},\bm{Q}^{\top}\bm{1}_{K}=\frac{1}{L}\bm{1}_{L}\\}$
is the set of the doubly-stochastic matrices having marginals
$\frac{1}{K}\bm{1}_{K}$ and $\frac{1}{L}\bm{1}_{L}$,
$\bm{D}=[d(\widehat{\bm{N}}_{k},\bm{N}_{l})]\in\mathbb{R}^{K\times L}$ is a
distance matrix, whose element measures the distance between the sequences.

$\bm{Q}^{*}=\arg\min_{\bm{Q}\in\Pi(\frac{1}{K}\bm{1}_{K},\frac{1}{L}\bm{1}_{L})}\langle\bm{D},\bm{Q}\rangle$,
the optimizer of (9), is the optimal transport matrix between the two sets of
event sequences. When $\widehat{\mathcal{N}}$ and $\mathcal{N}$ correspond to
the sets of generated event sequences and the real ones, this matrix is the
desired joint distribution, $i.e.$,
$\bm{Q}^{*}=[q^{*}(\widehat{\bm{N}}_{k},\bm{N}_{l})]$. This optimization
problem can be solved by many efficient methods, $e.g.$, the Sinkhorn scaling
method (Cuturi, 2013) and the proximal point method (Xie et al., 2020).

Figure 2: An illustration of the hierarchical optimal transport distance
between two sets of event sequences.

For $\widehat{\bm{N}}_{k}=\\{\widehat{N}_{u}^{k}\\}_{u\in\mathcal{V}_{k}}$ and
$\bm{N}_{l}=\\{N_{v}^{l}\\}_{v\in\mathcal{V}_{l}}$, where $\mathcal{V}_{k}$
and $\mathcal{V}_{l}$ are the sets of their event types, we also implement
their distance as an optimal transport distance:
$d(\widehat{\bm{N}}_{k},\bm{N}_{l})$

$\displaystyle\begin{aligned}
:=&\sideset{}{{}_{\bm{T}\in\Pi\left(\frac{1}{|\mathcal{V}_{k}|}\bm{1}_{|\mathcal{V}_{k}|},\frac{1}{|\mathcal{V}_{l}|}\bm{1}_{|\mathcal{V}_{l}|}\right)}}{\min}\sideset{}{{}_{u,v}}{\sum}T_{uv}d(\widehat{N}_{u}^{k},N_{v}^{l})\\\
=&\sideset{}{{}_{\bm{T}\in\Pi\left(\frac{1}{|\mathcal{V}_{k}|}\bm{1}_{|\mathcal{V}_{k}|},\frac{1}{|\mathcal{V}_{l}|}\bm{1}_{|\mathcal{V}_{l}|}\right)}}{\min}\langle\bm{D}_{kl},\bm{T}\rangle,\end{aligned}$
(10)

where
$\bm{D}_{kl}=[d(\widehat{N}_{u}^{k},N_{v}^{l})]\in\mathbb{R}^{|\mathcal{V}_{k}|\times|\mathcal{V}_{l}|}$
is the distance matrix for $\widehat{\bm{N}}_{k}$ and $\bm{N}_{l}$, and
$d(\widehat{N}_{u}^{k},N_{v}^{l})=\frac{1}{T}\int_{0}^{T}|\widehat{N}_{u}^{k}(t)-N_{v}^{l}(t)|\text{d}t$
measures the difference between the sequence of the type-$u$ events and that
of the type-$v$ events in $[0,T]$.

Plugging (10) into (9), we measure the difference between two sets of
heterogeneous event sequences by a hierarchical optimal transport distance, in
which the ground distance used in (9) is also an optimal transport distance.
Figure 2 illustrates the hierarchical optimal transport distance. In the
proposed HOT distance, the optimal transport matrix $\bm{Q}^{*}$ derived by
(9) achieves a soft alignment between the generated sequences and the real
ones, which corresponds to the joint distribution in (8). Additionally, the
optimal transport matrix $\bm{T}^{*}$ derived by (10) aligns the event types
of a generated sequence with those of a real one, which indicates the
correspondence of real-world event types in the latent event type space. In
Section 5, we will show that based on $\bm{T}^{*}$ GHP can generate
semantically-meaningful Hawkes processes and their event sequences.

###  3.3 Further analysis

Our HOT-based RAML method (denoted as RAML-HOT) has two advantages over the
original RAML in (Norouzi et al., 2016). Firstly, the reward used in the
original RAML is the sum of the conditional probabilities, $i.e.$,
$\sum_{\bm{N}_{l}\in\mathcal{N}}q(\hat{\bm{N}}_{k}|\bm{N}_{l})$. Accordingly,
a generated sequence earns a high reward only when it is similar to most real
sequences. This setting is unnecessary even unreasonable in our problem: a
generated sequence is likely to close to a small number of real sequences
because the real sequences are heterogeneous and yield different Hawkes
processes. Secondly, the original RAML implements $q(\widehat{\bm{N}}|\bm{N})$
empirically as an exponential pay-off distribution, $i.e.$,
$q(\widehat{\bm{N}}|\bm{N})=\frac{1}{Z}\exp(\frac{r(\widehat{\bm{N}},\bm{N})}{\tau})$,
where $Z$ is the normalizing constant, $\tau$ is the hyperparameter, and
$r(\widehat{\bm{N}},\bm{N})$ is a predefined reward function. Different from
the original RAML, our RAML-HOT method computes the joint distribution
$q(\widehat{\bm{N}},\bm{N})$ based on the HOT distance and the reward
$\max_{\bm{N}_{l}\in\mathcal{N}}q(\widehat{\bm{N}}_{k},\bm{N}_{l})$ is more
reasonable and interpretable.

Algorithm 1 shows the steps of our learning method and the original RAML when
learning a GHP model. Given $L$ real sequences, for each of them we denote
$\mathcal{O}(V)$ as the number of its event types and $\mathcal{O}(I)$ the
number of events per event type. When learning GHP, RAML-HOT generates a batch
of sequences and computes its HOT distance to a batch of real sequences.
Because of solving $B^{2}+1$ optimal transport problems, its computational
complexity is $\mathcal{O}(B^{2}IV^{2})$, where $B$ is batch size. Regarding
computational cost, GHP is suitable for modeling sparse heterogeneous event
sequences, in which both $V$ and $I$ are small and thus our RAML-HOT method is
efficient. Such sequences are common in real-world applications: $i$) The
admissions of different patients in a hospital cover many kinds of diseases,
but each patient often has a limited number of diseases and admissions. $ii$)
The Linkedin users cover many types of jobs, but each user has few job-hopping
behaviors among a small number of jobs. In such situations, GHP captures the
point process per sequence, whose number of event types ($i.e.$, $V$) is
limited. Compared to modeling a large point process model for all the
sequences, applying our GHP model can mitigate the risk of over-fitting.

##  4 Related Work

Hawkes processes. Because of its quantitative power and good interpretability,
Hawkes process has been a significant tool for event sequence analysis and
achieved encouraging performance in many applications like social network
analysis (Zhou et al., 2013; Farajtabar et al., 2017) and financial
engineering (Bacry et al., 2015). These years, many efforts have been made to
develop the variants of Hawkes process, $e.g.$, the mixture model of Hawkes
processes (Xu & Zha, 2017), the recurrent neural networks in the continuous
time (Du et al., 2016; Mei & Eisner, 2017) and the Hawkes processes with
attention mechanisms (Zhang et al., 2020; Zuo et al., 2020). Most existing
models are learned by the maximum likelihood estimation. Recently, more
cutting-edge techniques are applied, $e.g.$, Wasserstein generative
adversarial network (Xiao et al., 2017), reinforcement learning (Li et al.,
2018), and noisy contrastive estimation (Mei et al., 2020). However, most
existing methods cannot learn multiple Hawkes processes with different event
types.

Algorithm 1 Learning a GHP model

1: Input Real event sequences $\mathcal{N}$.

2: Initialize the model parameter $\theta$ randomly.

3: for each epoch

4: for each batch of real sequences
$\\{\bm{N}_{b}\\}_{b=1}^{B}\subset\mathcal{N}$

5: Generate $B$ sequences $\\{\widehat{\bm{N}}_{b}\\}_{b=1}^{B}$ via (4).

6: Calculate $d(\widehat{\bm{N}}_{b},\bm{N}_{b^{\prime}})$ by (10) and obtain
the matrix $\bm{D}$.

7: RAML: Set the reward function
$r(\widehat{\bm{N}}_{b},\bm{N}_{b^{\prime}})=-d(\widehat{\bm{N}}_{b},\bm{N}_{b^{\prime}})$
and $q(\widehat{\bm{N}}_{b}|\bm{N}_{b^{\prime}})$ an exponential pay-off
distribution.

8: Our RAML-HOT: Solve (9) and obtain
$\bm{Q}^{*}=[q^{*}(\widehat{\bm{N}}_{b},\bm{N}_{b^{\prime}})]$.

9: Calculate the loss function in (8).

10: Update $\theta$ by the Adam algorithm (Kingma & Ba, 2014).

Graphons. Graphon is a nonparametric graph model generating arbitrary-size
graphs in an infinite dimensional space (Lovász, 2012). Given observed graphs,
most existing methods learn graphons as stochastic block models (Channarond et
al., 2012; Airoldi et al., 2013; Chan & Airoldi, 2014), low-rank matrices
(Keshavan et al., 2010; Chatterjee et al., 2015; Xu, 2018) or Gromov-
Wasserstein barycenters (Xu et al., 2020), which approximate graphons by 2D
step functions based on the weak regularity lemma (Frieze & Kannan, 1999).

Optimal transport. The theory of optimal transport (Villani, 2008) has been
widely used in distribution estimation (Boissard et al., 2015) and matching
(Courty et al., 2017), and data generation (Arjovsky et al., 2017). Because of
its usefulness, many methods have been proposed to compute the optimal
transport efficiently, $e.g.$, the Sinkhorn scaling algorithm (Cuturi, 2013)
and its stochastic variant (Altschuler et al., 2017), the Bregman ADMM
algorithm (Wang & Banerjee, 2014), the proximal point method (Xie et al.,
2020), and the sliced Wasserstein distance (Kolouri et al., 2018). Recently,
hierarchical optimal transport (HOT) models are proposed in (Lee et al., 2019;
Yurochkin et al., 2019), which achieve encouraging performance on data
clustering. Our work makes the first attempt to introduce the HOT model into
event sequence analysis.

(a) $d_{\text{fgw}}(\hat{\theta},\theta)$

(b) $d_{\text{ot}}(\widehat{\mathcal{N}},\mathcal{N})$

(c) The influence of $B$

(d) The influence of $V_{\max}$

Figure 3: The experimental results on synthetic data. The runtime in (c, d) is
derived by running on a CPU.

##  5 Experiments

###  5.1 Experiments on synthetic data

To test our learning method, we first learn GHP models from synthetic
heterogeneous event sequences. The synthetic sequences are generated by a
predefined GHP model. For the predefined model, we set $V_{\max}=20$, the
decay kernel $\kappa(t)=\exp(-t)$, the number of Fourier bases ($i.e.$, the
$S$ in (3)) of $g(x,y)$ as $5$, and sampled the model parameter $\theta$ from
a multivariate normal distribution. Given the predefined model, we simulate
$120$ event sequences in the time window $[0,50]$ by the steps in (4), in
which we apply the Ogata’s thinning method (Ogata, 1981). We select $100$
sequences for learning a new GHP model, $10$ sequences for validation, and the
remaining $10$ sequences for testing the learned model. We evaluate the
learned model based on two criteria. Firstly, we compute the Fused Gromov-
Wasseserstein (FGW) distance (Vayer et al., 2018) between the estimated model
parameter $\hat{\theta}$ and the ground truth $\theta$:

$\displaystyle\begin{aligned}
d_{\text{fgw}}(\hat{\theta},\theta):=&\sideset{}{{}_{\pi\in\Pi}}{\inf}\mathbb{E}_{x,x^{\prime}\sim\pi}[|\hat{f}(x)-f(x^{\prime})|^{2}]+\\\
&\mathbb{E}_{x,x^{\prime},y,y^{\prime}\sim\pi\otimes\pi}[|\hat{g}(x,y)-g(x^{\prime},y^{\prime})|^{2}].\end{aligned}$
(11)

The FGW distance minimizes the expected error between the model parameters by
finding finds an optimal transport $\pi$, whose implementation is in Appendix
B. Secondly, we simulate a set of sequences based on the learned model and
calculate its HOT distance to the testing set, $i.e.$,
$d_{\text{ot}}(\widehat{\mathcal{N}},\mathcal{N})$.

Setting the number of training sequences from $10$ to $100$, we test our
learning method (RAML-HOT) and compare it with the original RAML (Norouzi et
al., 2016). For each method, we set the number of epochs to be $20$ and the
learning rate to be $0.01$. For our RAML-HOT method, we apply the Sinkhorn
scaling method (Cuturi, 2013) to compute the HOT distance. Figure 3(a) and
Figure 3(b) show the averaged performance of the two learning methods in $10$
trials. With the increase of training data, both our RAML-HOT and the RAML
improve their learning results consistently, achieving smaller
$d_{\text{fgw}}(\hat{\theta},\theta)$ and
$d_{\text{ot}}(\widehat{\mathcal{N}},\mathcal{N})$ with smaller standard
deviation. Moreover, we can find that our RAML-HOT method outperforms the RAML
method on the two measurements. This result verifies the feasibility of our
RAML-HOT method and demonstrates its advantages claimed in Section 3.3 — the
reward used in (7) is suitable for our problem, and leveraging the HOT
distance works better than using the exponential pay-off distribution.

For our RAML-HOT method, the batch size $B$ is a key hyperparameter.
Generally, using a large batch size may improve learning results. However, for
our method, whose computational complexity is quadratic to the batch size, we
need to carefully set the batch size to achieve a trade-off between
performance and efficiency. Figure 3(c) visualize the runtime per batch and
the $d_{\text{fgw}}(\hat{\theta},\theta)$ achieved by our method with respect
to different batch sizes. We find that the
$d_{\text{fgw}}(\hat{\theta},\theta)$ is relatively stable but the runtime
increases quadratically with respect to the batch size. According to the
result, we set $B=10$ in our experiments.

Besides the batch size, the maximum number of event types $V_{\max}$ is also
significant. According to (4), for the event sequences generated by our GHP
model, the expected number of their event types is
$\mathbb{E}[V]=\frac{\hat{V}_{\max}}{2}$. In the training phase, the maximum
number of event types used to learn the GHP model, denoted as
$\hat{V}_{\max}$, may be different from the ground truth $V_{\max}$. Setting
$\hat{V}_{\max}$ too large or too small may lead to the model misspecification
problem. As shown in Figure 3(d), the runtime of our method increases
quadratically with respect to $\hat{V}_{\max}$, which verifies the
computational complexity in Section 3.3. The best
$d_{\text{fgw}}(\hat{\theta},\theta)$ is achieved when the
$\hat{V}_{\max}=V_{\max}$. In practice, given a set of training sequences, we
calculate the averaged number of event types per sequence, denoted as
$\bar{V}$, and set $\hat{V}_{\max}=2\bar{V}$.

###  5.2 Modeling sparse heterogeneous event sequences

As aforementioned, our GHP model is suitable for modeling sparse heterogeneous
event sequences. We demonstrate the usefulness of our GHP model on two
representative real-world datasets. The first is the Linkedin dataset, which
contains the job-hopping and promotion behaviors of 2,439 Linkedin users (Xu
et al., 2017). The dataset has 3,730 kinds of jobs ($i.e.$, the event types).
However, most users seldom change their jobs, and each of their event
sequences contains 1 - 6 events in general. The second is the MIMIC-III
dataset. It contains 2,371 patients, each with more than two admissions in a
hospital (Johnson et al., 2016). The dataset covers 2,789 kinds of diseases,
but each patient suffers from extremely few of them and has a limited number
of admissions. Given these two datasets, we apply our RAML-HOT method to learn
GHP models and compare the models with state-of-the-art point process models.
Specifically, we consider six baselines: the classic Hawkes process (HP) (Zhou
et al., 2013), the time-varying Hawkes process (TVHP) (Xu et al., 2017), the
recurrent marked temporal point process (RMTPP) (Du et al., 2016), the neural
Hawkes process (NHP) (Mei & Eisner, 2017), the self-attentive Hawkes process
(SAHP) (Zhang et al., 2020), and the transformer Hawkes process (THP) (Zuo et
al., 2020). For our GHP model, we implement it to generate classic Hawkes
processes ($i.e.$, GHP${}_{\text{HP}}$) and extend it to generate time-varying
Hawkes processes ($i.e.$, GHP${}_{\text{TVHP}}$).

Table 1: Comparisons on real-world data Method | LinkedIn | MIMIC-III
---|---|---
NLL | $d_{ot}(\widehat{\mathcal{N}},\mathcal{N})$ | NLL | $d_{ot}(\widehat{\mathcal{N}},\mathcal{N})$
HP |  144.45${}_{\pm\text{20.70}}$ |  9.29${}_{\pm\text{1.38}}$ |  87.72${}_{\pm\text{7.73}}$ |  10.30${}_{\pm\text{0.69}}$
TVHP |  113.82${}_{\pm\text{8.34}}$ |  8.66${}_{\pm\text{1.57}}$ |  63.25${}_{\pm\text{3.08}}$ |  10.06${}_{\pm\text{0.63}}$
RMTPP |  127.39${}_{\pm\text{13.44}}$ |  8.83${}_{\pm\text{1.49}}$ |  82.46${}_{\pm\text{6.18}}$ |  11.76${}_{\pm\text{0.54}}$
NHP |  52.58${}_{\pm\text{14.52}}$ |  7.47${}_{\pm\text{1.26}}$ |  60.05${}_{\pm\text{5.27}}$ |  9.98${}_{\pm\text{0.71}}$
SAHP |  38.91${}_{\pm\text{10.33}}$ |  7.09${}_{\pm\text{0.80}}$ |  54.45${}_{\pm\text{3.12}}$ |  10.01${}_{\pm\text{0.95}}$
THP |  30.64${}_{\pm\text{7.03}}$ |  6.44${}_{\pm\text{0.61}}$ |  42.08${}_{\pm\text{5.26}}$ |  9.85${}_{\pm\text{0.88}}$
GHP${}_{\text{HP}}$ |  19.36${}_{\pm\text{2.97}}$ |  5.23${}_{\pm\text{0.28}}$ |  33.79${}_{\pm\text{6.54}}$ |  9.36${}_{\pm\text{2.45}}$
GHP${}_{\text{TVHP}}$ | 17.55${}_{\pm\text{2.61}}$ | 4.71${}_{\pm\text{0.15}}$ | 31.63${}_{\pm\text{5.83}}$ | 8.96${}_{\pm\text{2.29}}$

For each dataset, we train the models above based on 80% sequences and test
them on the remaining 20% sequences based on two measurements. Firstly, for
each model we can simulate a set of event sequences and calculate their
optimal transport distance to the testing set, $i.e.$,
$d_{ot}(\widehat{\mathcal{N}},\mathcal{N})$. Secondly, given the learned
method, we can calculate the negative log-likelihood (NLL) of the testing
sequences. When calculating $d_{ot}(\widehat{\mathcal{N}},\mathcal{N})$, our
GHP models apply the HOT distance based on (9, 10). The optimal transport
$\bm{Q}^{*}=[q^{*}(\bm{N}_{k},\bm{N}_{l})]$ derived by (9) helps match the
simulated sequences with the testing ones. For each pair of the sequence, the
optimal transport $\bm{T}^{*}=[T^{*}_{uv}]$ derived by (10) indicates the
correspondence between the event types of the testing sequence in the latent
event type space, $i.e.$, the latent event type
$\\{x_{1},...,x_{|\mathcal{V}|}\\}\in\Omega$ for the real-world event types
$\mathcal{V}$. For the $v$-th event type of the $l$-th testing sequence
$\bm{N}_{l}$, we first estimate the probability that it matches with the
$u$-th latent event type of the $k$-th generated sequence $\hat{\bm{N}}_{k}$
as $p(x_{u}^{k}|v)\propto T^{*}_{uv}q^{*}(\bm{N}_{k},\bm{N}_{l})$. Then, we
take $\\{x_{u}^{k}\\}$ as landmarks on $\Omega$ and approximate the
probability density $p(x|v)$ by the kernel density estimation, $i.e.$,
$p(x|v)=\frac{1}{Z}\sum_{u,k}p(x_{u}^{k}|v)\exp(-\frac{|x-x_{u}^{k}|^{2}}{2\sigma^{2}})$,
where $Z$ is the normalizing constant and $\sigma$ is the bandwidth of the
Gaussian kernel. For each event type in the testing sequence, we select its
latent event type corresponding to the largest $p(x|v)$, $i.e.$,
$x^{*}=\max_{x}p(x|v)$. Given the latent event types, we obtain the Hawkes
process from our GHP model and calculate the NLL of the testing sequence.
Table 1 shows the performance of various models in $10$ trials. In particular,
the baselines are learned as a single point process with a huge number of
event types from sparse event sequences, which have a high risk of over-
fitting. Our GHP models, on the contrary, describe each sparse event sequence
by a small point process sampled from an underlying graphon and learn the
point processes jointly. As a result, we can find that our GHP models
outperform the baselines consistently.

(a) Linkedin

(b) MIMIC-III

Figure 4: The graphons of the two real-world datasets.

In Figure 4, we show the probability densities of some representative real-
world event types in the latent space and check their triggering patterns on
the graphons. The graphons are visualized the resolution $200\times 200$. For
the Linkedin dataset, we find the pairs of “UCB, graduate student” and
“Google, research scientist” in the graphon according to their probability
density. The values of the pairs indicate that a graduate student at UC
Berkeley is likely to be a researcher at Google, while a researcher at Google
may not be go back to school, which we think is reasonable in practice. For
the MIMIC-III dataset, we find the pairs of “Hypertensive kidney disease” and
“Atherosclerosis” in the graphon. The values of the pairs reflect the
following facts: if a patient has atherosclerosis, his risk of having the
kidney disease caused by hypertensive will increase; however, the reasons for
the hypertensive kidney disease are complicated, and a patient having this
disease may not have atherosclerosis.

##  6 Conclusions

In this work, we propose a graphon-based Hawkes process model, capturing the
generative mechanism of multiple point processes with graph structures from
heterogeneous event sequences. Our GHP model is a new member of hierarchical
generative models for event sequence analysis. To our knowledge, it makes the
first attempt to combine graphon models with point processes. In the future,
we will improve GHP and its learning algorithm, $e.g.$, developing efficient
algorithms to compute the HOT distance with lower complexity and building the
GHP model based on deep neural networks. Additionally, the HOT distance used
in our model provides a potential solution to align event types of
heterogeneous event sequences. Combining GHP with existing event sequence
alignment methods (Xu et al., 2018; Trouleau et al., 2019; Luo et al., 2019),
we plan to develop a new framework of data fusion and augmentation for large-
scale point processes.

## References

  * Airoldi et al. (2013) Airoldi, E. M., Costa, T. B., and Chan, S. H.  Stochastic blockmodel approximation of a graphon: Theory and consistent estimation.  In _Advances in Neural Information Processing Systems_ , pp. 692–700, 2013.
  * Alaa et al. (2017) Alaa, A. M., Hu, S., and Schaar, M.  Learning from clinical judgments: Semi-markov-modulated marked hawkes processes for risk prognosis.  In _International Conference on Machine Learning_ , pp. 60–69, 2017\.
  * Altschuler et al. (2017) Altschuler, J., Weed, J., and Rigollet, P.  Near-linear time approximation algorithms for optimal transport via Sinkhorn iteration.  In _Advances in Neural Information Processing Systems_ , pp. 1964–1974, 2017.
  * Arjovsky et al. (2017) Arjovsky, M., Chintala, S., and Bottou, L.  Wasserstein GAN.  _arXiv preprint arXiv:1701.07875_ , 2017.
  * Bacry et al. (2015) Bacry, E., Mastromatteo, I., and Muzy, J.-F.  Hawkes processes in finance.  _Market Microstructure and Liquidity_ , 1(01):1550005, 2015.
  * Blundell et al. (2012) Blundell, C., Beck, J., and Heller, K. A.  Modelling reciprocating relationships with hawkes processes.  In _Advances in Neural Information Processing Systems_ , pp. 2600–2608, 2012.
  * Boissard et al. (2015) Boissard, E., Le Gouic, T., Loubes, J.-M., et al.  Distribution’s template estimate with Wasserstein metrics.  _Bernoulli_ , 21(2):740–759, 2015.
  * Chan & Airoldi (2014) Chan, S. and Airoldi, E.  A consistent histogram estimator for exchangeable graph models.  In _International Conference on Machine Learning_ , pp. 208–216, 2014.
  * Channarond et al. (2012) Channarond, A., Daudin, J.-J., Robin, S., et al.  Classification and estimation in the stochastic blockmodel based on the empirical degrees.  _Electronic Journal of Statistics_ , 6:2574–2601, 2012\.
  * Chatterjee et al. (2015) Chatterjee, S. et al.  Matrix estimation by universal singular value thresholding.  _The Annals of Statistics_ , 43(1):177–214, 2015\.
  * Chiu et al. (2013) Chiu, S. N., Stoyan, D., Kendall, W. S., and Mecke, J.  _Stochastic geometry and its applications_.  John Wiley & Sons, 2013.
  * Courty et al. (2017) Courty, N., Flamary, R., and Ducoffe, M.  Learning Wasserstein embeddings.  _arXiv preprint arXiv:1710.07457_ , 2017.
  * Cuturi (2013) Cuturi, M.  Sinkhorn distances: Lightspeed computation of optimal transport.  In _Advances in neural information processing systems_ , pp. 2292–2300, 2013.
  * Du et al. (2016) Du, N., Dai, H., Trivedi, R., Upadhyay, U., Gomez-Rodriguez, M., and Song, L.  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_ , pp. 1555–1564, 2016.
  * Eichler et al. (2017) Eichler, M., Dahlhaus, R., and Dueck, J.  Graphical modeling for multivariate hawkes processes with nonparametric link functions.  _Journal of Time Series Analysis_ , 38(2):225–242, 2017.
  * Farajtabar et al. (2017) Farajtabar, M., Wang, Y., Gomez-Rodriguez, M., Li, S., Zha, H., and Song, L.  Coevolve: A joint point process model for information diffusion and network evolution.  _The Journal of Machine Learning Research_ , 18(1):1305–1353, 2017.
  * Frieze & Kannan (1999) Frieze, A. and Kannan, R.  Quick approximation to matrices and applications.  _Combinatorica_ , 19(2):175–220, 1999.
  * Gao & Caines (2019) Gao, S. and Caines, P. E.  Graphon control of large-scale networks of linear systems.  _IEEE Transactions on Automatic Control_ , 2019.
  * Hawkes (1971) Hawkes, A. G.  Spectra of some self-exciting and mutually exciting point processes.  _Biometrika_ , pp. 83–90, 1971.
  * Johnson et al. (2016) Johnson, A. E., Pollard, T. J., Shen, L., Li-wei, H. L., Feng, M., Ghassemi, M., Moody, B., Szolovits, P., Celi, L. A., and Mark, R. G.  MIMIC-III, a freely accessible critical care database.  _Scientific data_ , 3:160035, 2016.
  * Keshavan et al. (2010) Keshavan, R. H., Montanari, A., and Oh, S.  Matrix completion from a few entries.  _IEEE transactions on information theory_ , 56(6):2980–2998, 2010.
  * Kingma & Ba (2014) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kolouri et al. (2018) Kolouri, S., Pope, P. E., Martin, C. E., and Rohde, G. K.  Sliced-wasserstein autoencoder: an embarrassingly simple generative model.  _arXiv preprint arXiv:1804.01947_ , 2018.
  * Lee et al. (2019) Lee, J., Dabagia, M., Dyer, E., and Rozell, C.  Hierarchical optimal transport for multimodal distribution alignment.  In _Advances in Neural Information Processing Systems_ , pp. 13453–13463, 2019.
  * Li et al. (2018) Li, S., Xiao, S., Zhu, S., Du, N., Xie, Y., and Song, L.  Learning temporal point processes via reinforcement learning.  In _Advances in neural information processing systems_ , pp. 10781–10791, 2018.
  * Lin et al. (2016) Lin, P., Guo, T., Wang, Y., Chen, F., et al.  Infinite hidden semi-markov modulated interaction point process.  _Advances in Neural Information Processing Systems_ , 29:3900–3908, 2016.
  * Liniger (2009) Liniger, T. J.  _Multivariate hawkes processes_.  PhD thesis, ETH Zurich, 2009.
  * Lovász (2012) Lovász, L.  _Large networks and graph limits_ , volume 60.  American Mathematical Soc., 2012.
  * Luo et al. (2015) Luo, D., Xu, H., Zhen, Y., Ning, X., Zha, H., Yang, X., and Zhang, W.  Multi-task multi-dimensional hawkes processes for modeling event sequences.  2015\.
  * Luo et al. (2019) Luo, D., Xu, H., and Carin, L.  Fused gromov-wasserstein alignment for hawkes processes.  _arXiv preprint arXiv:1910.02096_ , 2019.
  * Mei & Eisner (2017) Mei, H. and Eisner, J. M.  The neural hawkes process: A neurally self-modulating multivariate point process.  _Advances in Neural Information Processing Systems_ , 30:6754–6764, 2017.
  * Mei et al. (2020) Mei, H., Wan, T., and Eisner, J.  Noise-contrastive estimation for multivariate point processes.  _arXiv preprint arXiv:2011.00717_ , 2020.
  * Mémoli (2011) Mémoli, F.  Gromov-Wasserstein distances and the metric approach to object matching.  _Foundations of computational mathematics_ , 11(4):417–487, 2011.
  * Møller & Rasmussen (2006) Møller, J. and Rasmussen, J. G.  Approximate simulation of hawkes processes.  _Methodology and Computing in Applied Probability_ , 8(1):53–64, 2006.
  * Norouzi et al. (2016) Norouzi, M., Bengio, S., Jaitly, N., Schuster, M., Wu, Y., Schuurmans, D., et al.  Reward augmented maximum likelihood for neural structured prediction.  In _Advances In Neural Information Processing Systems_ , pp. 1723–1731, 2016.
  * Ogata (1981) Ogata, Y.  On lewis’ simulation method for point processes.  _IEEE transactions on information theory_ , 27(1):23–31, 1981.
  * Peyré et al. (2019) Peyré, G., Cuturi, M., et al.  Computational optimal transport.  _Foundations and Trends® in Machine Learning_ , 11(5-6):355–607, 2019.
  * Rabin et al. (2011) Rabin, J., Peyré, G., Delon, J., and Bernot, M.  Wasserstein barycenter and its application to texture mixing.  In _International Conference on Scale Space and Variational Methods in Computer Vision_ , pp. 435–446. Springer, 2011.
  * Tank et al. (2018) Tank, A., Covert, I., Foti, N., Shojaie, A., and Fox, E.  Neural granger causality for nonlinear time series.  _arXiv preprint arXiv:1802.05842_ , 2018.
  * Trouleau et al. (2019) Trouleau, W., Etesami, J., Grossglauser, M., Kiyavash, N., and Thiran, P.  Learning hawkes processes under synchronization noise.  In _International Conference on Machine Learning_ , pp. 6325–6334. PMLR, 2019.
  * Van Loan & Golub (1983) Van Loan, C. F. and Golub, G. H.  _Matrix computations_.  Johns Hopkins University Press Baltimore, 1983.
  * Vayer et al. (2018) Vayer, T., Chapel, L., Flamary, R., Tavenard, R., and Courty, N.  Fused Gromov-Wasserstein distance for structured objects: theoretical foundations and mathematical properties.  _arXiv preprint arXiv:1811.02834_ , 2018.
  * Villani (2008) Villani, C.  _Optimal transport: old and new_ , volume 338.  Springer Science & Business Media, 2008.
  * Wang & Banerjee (2014) Wang, H. and Banerjee, A.  Bregman alternating direction method of multipliers.  In _NeurIPS_ , 2014.
  * Wang et al. (2016) Wang, Y., Du, N., Trivedi, R., and Song, L.  Coevolutionary latent feature processes for continuous-time user-item interactions.  _Advances in neural information processing systems_ , 29:4547–4555, 2016.
  * Xiao et al. (2017) Xiao, S., Farajtabar, M., Ye, X., Yan, J., Song, L., and Zha, H.  Wasserstein learning of deep generative point process models.  In _Advances in neural information processing systems_ , pp. 3247–3257, 2017.
  * Xie et al. (2020) Xie, Y., Wang, X., Wang, R., and Zha, H.  A fast proximal point method for computing exact Wasserstein distance.  In _Uncertainty in Artificial Intelligence_ , pp. 433–453. PMLR, 2020.
  * Xu & Zha (2017) Xu, H. and Zha, H.  A dirichlet mixture model of hawkes processes for event sequence clustering.  _Advances in Neural Information Processing Systems_ , 30:1354–1363, 2017.
  * Xu et al. (2016a) Xu, H., Farajtabar, M., and Zha, H.  Learning granger causality for hawkes processes.  In _International conference on machine learning_ , pp. 1717–1726, 2016a.
  * Xu et al. (2016b) Xu, H., Wu, W., Nemati, S., and Zha, H.  Patient flow prediction via discriminative learning of mutually-correcting processes.  _IEEE transactions on Knowledge and Data Engineering_ , 29(1):157–171, 2016b.
  * Xu et al. (2017) Xu, H., Luo, D., and Zha, H.  Learning hawkes processes from short doubly-censored event sequences.  In _International Conference on Machine Learning_ , pp. 3831–3840, 2017.
  * Xu et al. (2018) Xu, H., Carin, L., and Zha, H.  Learning registered point processes from idiosyncratic observations.  In _International Conference on Machine Learning_ , pp. 5443–5452. PMLR, 2018.
  * Xu et al. (2019) Xu, H., Luo, D., Zha, H., and Carin, L.  Gromov-Wasserstein learning for graph matching and node embedding.  In _International Conference on Machine Learning_ , pp. 6932–6941, 2019.
  * Xu et al. (2020) Xu, H., Luo, D., Carin, L., and Zha, H.  Learning graphons via structured gromov-wasserstein barycenters.  _arXiv preprint arXiv:2012.05644_ , 2020.
  * Xu (2018) Xu, J.  Rates of convergence of spectral methods for graphon estimation.  In _International Conference on Machine Learning_ , pp. 5433–5442, 2018.
  * Yurochkin et al. (2019) Yurochkin, M., Claici, S., Chien, E., Mirzazadeh, F., and Solomon, J.  Hierarchical optimal transport for document representation.  _arXiv preprint arXiv:1906.10827_ , 2019.
  * Zhang et al. (2020) Zhang, Q., Lipani, A., Kirnap, O., and Yilmaz, E.  Self-attentive hawkes process.  In _International Conference on Machine Learning_ , pp. 11183–11193. PMLR, 2020.
  * Zhao et al. (2015) Zhao, Q., Erdogdu, M. A., He, H. Y., Rajaraman, A., and Leskovec, J.  Seismic: A self-exciting point process model for predicting tweet popularity.  In _Proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining_ , pp. 1513–1522, 2015.
  * Zhou et al. (2013) Zhou, K., Zha, H., and Song, L.  Learning triggering kernels for multi-dimensional hawkes processes.  In _International Conference on Machine Learning_ , pp. 1301–1309, 2013.
  * Zhu (2013) Zhu, L.  _Nonlinear Hawkes Processes_.  PhD thesis, New York University, 2013.
  * Zuo et al. (2020) Zuo, S., Jiang, H., Li, Z., Zhao, T., and Zha, H.  Transformer hawkes process.  _arXiv preprint arXiv:2002.09291_ , 2020.

##  Appendix A The Properties of Our GHP model

###  A.1 The proof of Property 2.1

Property 2.1
_$\text{HP}_{\mathcal{V}}(\bm{\mu},\bm{A})\sim\text{GHP}_{\Omega}(f,g)$ is
asymptotically stationary as long as $|\mathcal{V}|\leq V_{\max}$._

###### Proof.

For a classical shift-invariant Hawkes process, its intensity function is

$\displaystyle\begin{aligned}
\lambda_{v}(t)=\mu_{v}+\sideset{}{{}_{t_{i}<t}}{\sum}\phi_{vv^{\prime}}(t-t_{i})=\mu_{v}+\sideset{}{{}_{t_{i}<t}}{\sum}a_{vv^{\prime}}\eta(t-t_{i})~{}\text{for}~{}v\in\mathcal{V}.\end{aligned}$
(12)

We can construct a matrix
$\bm{\Phi}=[\phi_{vv^{\prime}}]\in\mathbb{R}^{V\times V}$, whose element is
$\phi_{vv^{\prime}}=\int_{0}^{\infty}\phi_{vv^{\prime}}(t)\text{d}t$.
According to the Proposition 1 in (Bacry et al., 2015), the Hawkes process is
asymptotically stationary if the impact functions satisfy:

$\displaystyle\begin{aligned}
1)~{}\forall~{}v,v^{\prime}\in\mathcal{V},~{}\phi_{vv^{\prime}}(t)\begin{cases}\geq
0,&t\geq 0,\\\ =0,&t<0.\end{cases}\quad
2)~{}\|\bm{\Phi}\|_{2}<1.\end{aligned}$ (13)

When setting $\phi_{vv^{\prime}}(t)=a_{vv^{\prime}}\kappa(t)$, as shown in
(12), we have

$\displaystyle\begin{aligned}
\bm{\Phi}=D\bm{A},~{}\text{where}~{}D=\int_{0}^{\infty}\eta(t)\text{d}t.\end{aligned}$
(14)

In the generative process of our GHP model, we have

$\displaystyle\begin{aligned}
a_{vv^{\prime}}=\frac{1}{V_{\max}D}g(x_{v},x_{v^{\prime}}),~{}\text{where}~{}V_{\max}\geq
V~{}\text{and}~{}g:\Omega\mapsto[0,1).\end{aligned}$ (15)

Accordingly, we have

$\displaystyle\begin{aligned}
\|\bm{A}\|_{2}\leq\|\bm{A}\|_{F}=\Bigl{(}\sideset{}{{}_{v,v^{\prime}}}{\sum}a_{vv^{\prime}}^{2}\Bigr{)}^{\frac{1}{2}}=\frac{1}{V_{\max}D}\Bigl{(}\sideset{}{{}_{v,v^{\prime}}}{\sum}g^{2}(x_{v},x_{v^{\prime}})\Bigr{)}^{\frac{1}{2}}<\frac{V}{V_{\max}D}\leq\frac{1}{D}.\end{aligned}$
(16)

Here, the first inequality is based on the relationship between a matrix’s
spectral norm and its Frobenius norm. The second (strict) inequality is based
on Assumption 2.2 ($i.e.$, $0\leq g(x,y)<1$ $\forall~{}(x,y)\in\Omega^{2}$).
Plugging (16) into (14), we have $\|\bm{\Phi}\|_{2}<1$, thus the stability
condition in (13) is satisfied. ∎

###  A.2 The proof of Property 2.3

Before proving Property 2.3, we first introduce the definition of the discrete
Wasserstein distance (the earth mover’s distance) and that of the discrete
Gromov-Wasserstein distance.

######  Definition A.1 (Earth Mover’s Distance).

Given $\bm{a}=\\{\bm{a}_{m}\in\mathbb{R}^{D}\\}_{m=1}^{M}$ and
$\bm{b}=\\{\bm{b}_{n}\in\mathbb{R}^{D}\\}_{n=1}^{N}$, the discrete Wasserstein
distance between them is

$\displaystyle\begin{aligned}
d_{\text{w}}(\bm{a},\bm{b}):=\sideset{}{{}_{\bm{T}\in\Pi\left(\frac{1}{M}\bm{1}_{M},\frac{1}{N}\bm{1}_{N}\right)}}{\min}\Bigl{(}\sideset{}{{}_{m,n}}{\sum}t_{mn}\|\bm{a}_{m}-\bm{b}_{n}\|_{2}^{2}\Bigr{)}^{\frac{1}{2}}=\sideset{}{{}_{\bm{T}\in\Pi\left(\frac{1}{M}\bm{1}_{M},\frac{1}{N}\bm{1}_{N}\right)}}{\min}\langle\bm{D},\bm{T}\rangle^{\frac{1}{2}}.\end{aligned}$
(17)

Here, $\bm{T}=[t_{mn}]$ is a doubly stochastic matrix in the set

$\Pi\left(\frac{1}{M}\bm{1}_{M},\frac{1}{N}\bm{1}_{N}\right)=\left\\{\bm{T}=[t_{mn}]\;\bigg{|}\;t_{mn}\geq
0,\bm{T}\bm{1}_{N}=\frac{1}{M}\bm{1}_{M},\bm{T}^{\top}\bm{1}_{M}=\frac{1}{N}\bm{1}_{N}\right\\},$

where $\bm{1}_{N}$ is the $N$-dimensional all-one vector. $\bm{D}=[d_{mn}]$ is
a distance matrix, whose element $d_{mn}=\|\bm{a}_{m}-\bm{b}_{n}\|_{2}^{2}$.
The optimal $\bm{T}$ corresponding to the distance, $i.e.$,
$\bm{T}^{*}=\arg\min_{\bm{T}\in\Pi(\frac{1}{M}\bm{1}_{M},\frac{1}{N}\bm{1}_{N})}\langle\bm{D},\bm{T}\rangle^{\frac{1}{2}}$,
is the so-called optimal transport matrix.

######  Definition A.2 (Discrete Gromov-Wasserstein Distance).

Given $\bm{a}=\\{\bm{a}_{m}\in\mathbb{R}^{D}\\}_{m=1}^{M}$ and
$\bm{b}=\\{\bm{b}_{n}\in\mathbb{R}^{D^{\prime}}\\}_{n=1}^{N}$, where $D$ can
be different from $D^{\prime}$, the discrete Gromov-Wasserstein distance
between them is

$\displaystyle\begin{aligned}
d_{\text{gw}}(\bm{a},\bm{b}):=\sideset{}{{}_{\bm{T}\in\Pi(\frac{1}{M}\bm{1}_{M},\frac{1}{N}\bm{1}_{N})}}{\min}\Bigl{(}\sideset{}{{}_{m,m^{\prime},n,n^{\prime}}}{\sum}t_{mn}t_{m^{\prime}n^{\prime}}|\|\bm{a}_{m}-\bm{a}_{m^{\prime}}\|_{2}-\|\bm{b}_{n}-\bm{b}_{n^{\prime}}\|_{2}|^{2}\Bigr{)}^{\frac{1}{2}}.\end{aligned}$
(18)

Similar to the discrete Wasserstein distance, $\bm{T}=[t_{mn}]$ is a doubly
stochastic matrix in the set
$\Pi\left(\frac{1}{M}\bm{1}_{M},\frac{1}{N}\bm{1}_{N}\right)$ and the optimal
$\bm{T}$ corresponding to the distance is the optimal transport matrix. When
two distance matrices $\bm{A}=\\{a_{mm^{\prime}}\\}_{m,m^{\prime}=1}^{M}$ and
$\bm{B}=\\{b_{nn^{\prime}}\\}_{n,n^{\prime}=1}^{N}$ are provided directly,
where $a_{mm^{\prime}}=\|\bm{a}_{m}-\bm{a}_{m^{\prime}}\|_{2}$ and
$b_{nn^{\prime}}=\|\bm{b}_{n}-\bm{b}_{n^{\prime}}\|_{2}$, we can rewrite the
discrete Gromov-Wasserstein distance equivalently as

$\displaystyle\begin{aligned}
d_{\text{gw}}(\bm{A},\bm{B}):=\sideset{}{{}_{\bm{T}\in\Pi\left(\frac{1}{M}\bm{1}_{M},\frac{1}{N}\bm{1}_{N}\right)}}{\min}\Bigl{(}\sideset{}{{}_{m,m^{\prime},n,n^{\prime}}}{\sum}t_{mn}t_{m^{\prime}n^{\prime}}|a_{mm^{\prime}}-b_{nn^{\prime}}|^{2}\Bigr{)}^{\frac{1}{2}}.\end{aligned}$
(19)

Property 2.3 _For $\text{HP}_{\mathcal{V}}(\bm{\mu}_{1},\bm{A}_{1})$ and
$\text{HP}_{\mathcal{U}}(\bm{\mu}_{2},\bm{A}_{2})\sim\text{GHP}_{\Omega}(f,g)$,
where $\text{GHP}_{\Omega}(f,g)$ satisfies Assumption 2.2, their parameters
satisfy_

$\displaystyle\begin{aligned}
&C_{1}^{f}d_{\text{w}}(\bm{x}_{1},\bm{x}_{2})\leq
d_{\text{w}}(\bm{\mu}_{1},\bm{\mu}_{2})\leq
C_{2}^{f}d_{\text{w}}(\bm{x}_{1},\bm{x}_{2}),\\\
&d_{\text{w}}(\bm{A}_{1},\bm{A}_{2})\leq
C^{g}d_{\text{w}}(\bm{x}_{1}^{\times},\bm{x}_{2}^{\times}),\\\
&d_{\text{gw}}(\bm{A}_{1},\bm{A}_{2})\leq
C^{g}d_{\text{gw}}(\bm{x}_{1}^{\times},\bm{x}_{2}^{\times}),\end{aligned}$

_where $\bm{x}_{1}=\\{x_{v,1}\\}_{v=1}^{|\mathcal{V}|}$ and
$\bm{x}_{2}=\\{x_{u,2}\\}_{u=1}^{|\mathcal{U}|}$ are the latent event types,
and
$\bm{x}_{1}^{\times}=\\{[x_{v,1};x_{v^{\prime},1}]\\}_{v,v^{\prime}=1}^{|\mathcal{V}|}$
and
$\bm{x}_{2}^{\times}=\\{[x_{u,2};x_{u^{\prime},2}]\\}_{u,u^{\prime}=1}^{|\mathcal{U}|}$
enumerate the pairs of the latent event types._

###### Proof.

Denote $|\mathcal{V}|=V$ and $|\mathcal{U}|=U$. Denote $\bm{T}^{x}$ as the
optimal transport matrix corresponding to
$d_{\text{w}}(\bm{x}_{1},\bm{x}_{2})$, $i.e.$,
$\bm{T}^{x}=\arg\min_{\bm{T}\in\Pi\left(\frac{1}{V}\bm{1}_{V},\frac{1}{U}\bm{1}_{U}\right)}(\sum_{v,u}t_{vu}|x_{v,1}-x_{u,2}|^{2})^{\frac{1}{2}}=\arg\min_{\bm{T}\Pi\left(\frac{1}{V}\bm{1}_{V},\frac{1}{U}\bm{1}_{U}\right)}\langle\bm{D}_{x},\bm{T}\rangle^{\frac{1}{2}}$,
where $\bm{D}_{x}=[d^{x}_{vu}]$ and $d^{x}_{vu}=|x_{v,1}-x_{u,2}|^{2}$.
Similarly, denote $\bm{T}^{\mu}$ as the optimal transport matrix corresponding
to $d_{\text{w}}(\bm{\mu}_{1},\bm{\mu}_{2})$, $i.e.$,
$d_{\text{w}}(\bm{\mu}_{1},\bm{\mu}_{2})=\langle\bm{D}_{\mu},\bm{T}^{\mu}\rangle$,
where $\bm{D}_{\mu}=[d^{\mu}_{vu}]$ and
$d^{\mu}_{vu}=|\mu_{v,1}-\mu_{u,2}|^{2}$. Obviously, we have

$\displaystyle\forall~{}\bm{T}\in\Pi\left(\frac{1}{V}\bm{1}_{V},\frac{1}{U}\bm{1}_{U}\right),\quad\langle\bm{D}_{x},\bm{T}\rangle^{\frac{1}{2}}\geq\langle\bm{D}_{x},\bm{T}^{x}\rangle^{\frac{1}{2}},\quad\text{and}\quad\langle\bm{D}_{\mu},\bm{T}\rangle^{\frac{1}{2}}\geq\langle\bm{D}_{\mu},\bm{T}^{\mu}\rangle^{\frac{1}{2}}.$
(20)

Based on (20), we have

$\displaystyle\begin{aligned}
d_{\text{w}}(\bm{\mu}_{1},\bm{\mu}_{2})&=\Bigl{(}\sideset{}{{}_{m,n}}{\sum}t_{vu}^{\mu}|\mu_{v,1}-\mu_{u,2}|^{2}\Bigr{)}^{\frac{1}{2}}\leq\Bigl{(}\sideset{}{{}_{v,u}}{\sum}t_{vu}^{x}|\mu_{v,1}-\mu_{u,2}|^{2}\Bigr{)}^{\frac{1}{2}}\\\
&\leq
C_{2}^{f}\Bigl{(}\sideset{}{{}_{v,u}}{\sum}t_{vu}^{x}|x_{v,1}-x_{u,2}|^{2}\Bigr{)}^{\frac{1}{2}}=C_{2}^{f}d_{\text{w}}(\bm{x}_{1},\bm{x}_{2}).\\\
d_{\text{w}}(\bm{\mu}_{1},\bm{\mu}_{2})&=\Bigl{(}\sideset{}{{}_{v,u}}{\sum}t_{vu}^{\mu}|\mu_{v,1}-\mu_{u,2}|^{2}\Bigr{)}^{\frac{1}{2}}\geq
C_{1}^{f}\Bigl{(}\sideset{}{{}_{v,u}}{\sum}t_{vu}^{\mu}|x_{v,1}-x_{u,2}|^{2}\Bigr{)}^{\frac{1}{2}}\\\
&\geq
C_{1}^{f}\Bigl{(}\sideset{}{{}_{v,u}}{\sum}t_{vu}^{x}|x_{v,1}-x_{u,2}|^{2}\Bigr{)}^{\frac{1}{2}}=C_{1}^{f}d_{\text{w}}(\bm{x}_{1},\bm{x}_{2}).\end{aligned}$
(21)

For $\bm{A}_{1}=[a_{vv^{\prime}}^{1}]$ and $\bm{A}_{2}=[a_{uu^{\prime}}^{2}]$,
their discrete Wasserstein distance is

$\displaystyle\begin{aligned}
d_{\text{w}}(\bm{A}_{1},\bm{A}_{2})=\sideset{}{{}_{\bm{T}\in\Pi(\frac{1}{V^{2}}\bm{1}_{V^{2}},\frac{1}{U^{2}}\bm{1}_{U^{2}})}}{\min}\langle\bm{D}_{A},\bm{T}\rangle^{\frac{1}{2}}=\langle\bm{D}_{A},\bm{T}^{A}\rangle^{\frac{1}{2}}\end{aligned}$
(22)

where $\bm{D}_{A}=[d_{vv^{\prime},uu^{\prime}}]\in\mathbb{R}^{V^{2}\times
U^{2}}$,
$d_{vv^{\prime},uu^{\prime}}=|a_{vv^{\prime}}^{1}-a_{uu^{\prime}}^{2}|^{2}$,
and $\bm{T}^{A}$ is the optimal transport matrix. Their discrete Gromov-
Wasserstein distance is

$\displaystyle\begin{aligned}
d_{\text{gw}}(\bm{A}_{1},\bm{A}_{2})=\sideset{}{{}_{\bm{T}\in\Pi(\frac{1}{V}\bm{1}_{V},\frac{1}{U}\bm{1}_{U})}}{\min}\langle\bm{D}_{A},\bm{T}\otimes\bm{T}\rangle^{\frac{1}{2}}=\langle\bm{D}_{A},\bm{T}_{\text{gw}}^{A}\otimes\bm{T}_{\text{gw}}^{A}\rangle^{\frac{1}{2}},\end{aligned}$
(23)

where $\bm{T}_{\text{gw}}^{A}$ is the optimal transport matrix and $\otimes$
represents the Kronecker multiplication between two matrices.

Similarly, we represent

$d_{\text{w}}(\bm{x}_{1}^{\times},\bm{x}_{2}^{\times})=\min_{\bm{T}\in\Pi(\frac{1}{V^{2}}\bm{1}_{V^{2}},\frac{1}{U^{2}}\bm{1}_{U^{2}})}\langle\bm{D}_{X},\bm{T}\rangle^{\frac{1}{2}},\quad
d_{\text{gw}}(\bm{x}_{1}^{\times},\bm{x}_{2}^{\times})=\min_{\bm{T}\in\Pi(\frac{1}{V}\bm{1}_{V},\frac{1}{U}\bm{1}_{U})}\langle\bm{D}_{X},\bm{T}\otimes\bm{T}\rangle^{\frac{1}{2}},$

respectively, where $\bm{D}_{X}=[D_{vv^{\prime},uu^{\prime}}]$ and
$D_{vv^{\prime},uu^{\prime}}=\|[x_{v,1};x_{v^{\prime},1}]-[x_{u,2};x_{u^{\prime},2}]\|_{2}^{2}$.
Accordingly, we denote $\bm{T}^{X}$ as the optimal transport matrix of
$d_{\text{w}}(\bm{x}_{1}^{\times},\bm{x}_{2}^{\times})$ and
$\bm{T}_{\text{gw}}^{X}$ as the optimal transport matrix of
$d_{\text{gw}}(\bm{x}_{1}^{\times},\bm{x}_{2}^{\times})$.

Similar to the derivation shown in (21), we have

$\displaystyle\begin{aligned}
&\langle\underbrace{\bm{D}_{A},\bm{T}^{A}\rangle^{\frac{1}{2}}}_{d_{\text{w}}(\bm{A}_{1},\bm{A}_{2})}\leq\langle\bm{D}_{A},\bm{T}^{X}\rangle^{\frac{1}{2}}\leq
C^{g}\underbrace{\langle\bm{D}_{X},\bm{T}^{X}\rangle^{\frac{1}{2}}}_{d_{\text{w}}(\bm{x}_{1}^{\times},\bm{x}_{2}^{\times})},\\\
&\langle\underbrace{\bm{D}_{A},\bm{T}_{\text{gw}}^{A}\otimes\bm{T}_{\text{gw}}^{A}\rangle^{\frac{1}{2}}}_{d_{\text{gw}}(\bm{A}_{1},\bm{A}_{2})}\leq\langle\bm{D}_{A},\bm{T}_{\text{gw}}^{X}\otimes\bm{T}_{\text{gw}}^{X}\rangle^{\frac{1}{2}}\leq
C^{g}\underbrace{\langle\bm{D}_{X},\bm{T}_{\text{gw}}^{X}\otimes\bm{T}_{\text{gw}}^{X}\rangle^{\frac{1}{2}}}_{d_{\text{gw}}(\bm{x}_{1}^{\times},\bm{x}_{2}^{\times})}.\end{aligned}$
(24)

∎

###  A.3 The proof of Property 2.4

######  Definition A.3 (Average Intensity (Bacry et al., 2015)).

For the stationary Hawkes process defined in (12), its counting process is
denoted as $\bm{N}(t)=\\{N_{v}(t)\\}_{v\in\mathcal{V},t\in[0,T]}$, where
$N_{v}(t)$ is the number of the type-$v$ events till time $t$, and its average
intensity is

$\displaystyle\bar{\bm{\lambda}}:=\frac{\mathbb{E}[\text{d}\bm{N}(t)]}{\text{d}t}=(\bm{I}_{V}-\bm{\Phi})^{-1}\bm{\mu}=(\bm{I}_{V}-D\bm{A})^{-1}\bm{\mu}.$
(25)

According to Campbell’s theorem (Chiu et al., 2013), given
$\bar{\bm{\lambda}}=[\bar{\lambda_{v}}]$, we have

$\displaystyle\mathbb{E}[\text{d}N_{v}(t)]=T\bar{\lambda}_{v}=\int_{0}^{T}\lambda_{v}(t)\text{d}t,~{}\forall~{}v\in\mathcal{V}.$
(26)

In other words, the average intensity reflects the overall dynamics of
different event types.

The proof of Property 2.4 is based on the theory of optimal transport and some
well-known theorems.

######  Property A.4 (Triangle inequality (Villani, 2008)).

For arbitrary $\bm{a}=\\{\bm{a}_{l}\in\mathbb{R}^{D}\\}_{l=1}^{L}$,
$\bm{b}=\\{\bm{b}_{m}\in\mathbb{R}^{D}\\}_{m=1}^{M}$, and
$\bm{c}=\\{\bm{c}_{n}\in\mathbb{R}^{D}\\}_{n=1}^{N}$, we have

$\displaystyle d_{\text{w}}(\bm{a},\bm{c})\leq
d_{\text{w}}(\bm{a},\bm{b})+d_{\text{w}}(\bm{b},\bm{c}).$ (27)

######  Theorem A.5 (One-dimensional Earth Mover’s Distance (Rabin et al.,
2011)).

For two sets of 1D points, $i.e.$, $\bm{a}=\\{a_{n}\in\mathbb{R}\\}_{n=1}^{N}$
and $\bm{b}=\\{b_{n}\in\mathbb{R}\\}_{n=1}^{N}$, their earth mover’s distance
has a closed form solution with complexity $\mathcal{O}(N\log N)$.

$\displaystyle\begin{aligned}
d_{\text{w}}(\bm{a},\bm{b})=\frac{1}{\sqrt{N}}\|\text{\rm
sort}(\bm{a})-\text{\rm
sort}(\bm{b})\|_{2}=\frac{1}{\sqrt{N}}\|\bm{a}-\bm{P}\bm{b}\|_{2},\end{aligned}$
(28)

where $\text{\rm sort}(\cdot)$ sorts the elements of a vector in a descending
order, and $\bm{P}\in\\{\bm{P}\in\\{0,1\\}^{N\times
N}~{}|~{}\bm{P}\bm{1}_{N}=\bm{1}_{N},\bm{P}^{\top}\bm{1}_{N}=\bm{1}_{N}\\}$ is
a permutation matrix, mapping the $n$-th largest element of $\bm{b}$ to the
$n$-th largest element of $\bm{a}$ for $n=1,...,N$. Obviously,
$\frac{1}{N}\bm{P}$ is the optimal transport matrix.

Additionally, according to the definition of the earth mover’s distance, we
have the following theorem:

######  Theorem A.6.

For a set of 1D points, $i.e.$, $\bm{a}=\\{a_{m}\in\mathbb{R}\\}_{m=1}^{M}$,
if we pad $N-M$ zeros to $\bm{a}$ and obtain
$\bm{a}^{\prime}=[\bm{a};\bm{0}_{N-M}]$, we have

$\displaystyle
d_{\text{w}}(\bm{a},\bm{a}^{\prime})\leq\sqrt{\frac{N-M}{MN}}\|\bm{a}\|_{2}$
(29)

###### Proof.

For $\bm{a}$ and $\bm{a}^{\prime}$, we obtain a distance matrix
$\bm{D}=[\bm{D}_{a},\bm{D}_{0}]\in\mathbb{R}^{M\times N}$. Here,
$\bm{D}_{a}=[d_{mm^{\prime}}]\in\mathbb{R}^{M\times M}$ and
$d_{mm^{\prime}}=|a_{m}-a_{m^{\prime}}|^{2}$. Obviously, the diagonal element
$d_{mm}=0$ for $m=1,..,M$. $\bm{D}_{0}=[d_{mn}]\in\mathbb{R}^{M\times(N-M)}$
and $d_{mn}=|a_{m}|^{2}$ for all $n=1,...,N-M$. Accordingly, we can design a
valid transport matrix
$\hat{\bm{T}}=[\frac{1}{N}\bm{I}_{M},\frac{1}{MN}\bm{1}_{M\times(N-M)}]$, such
that
$\hat{\bm{T}}\in\Pi\left(\frac{1}{M}\bm{1}_{M},\frac{1}{N}\bm{1}_{N}\right)$.
Therefore, we have

$\displaystyle\begin{aligned}
d_{\text{w}}(\bm{a},\bm{a}^{\prime})=\min_{\bm{T}\in\Pi\left(\frac{1}{M}\bm{1}_{M},\frac{1}{N}\bm{1}_{N}\right)}\langle\bm{D},\bm{T}\rangle^{\frac{1}{2}}\leq\langle\bm{D},\hat{\bm{T}}\rangle^{\frac{1}{2}}=\left(\underbrace{\sideset{}{{}_{m}}{\sum}d_{mm}}_{=0}+\frac{1}{MN}\sideset{}{{}_{m=1}^{M}}{\sum}\sideset{}{{}_{n=1}^{N-M}}{\sum}|a_{m}|^{2}\right)^{\frac{1}{2}}=\sqrt{\frac{N-M}{MN}}\|\bm{a}\|_{2}.\end{aligned}$

∎

######  Theorem A.7 (The perturbation theory of linear system (Van Loan &
Golub, 1983)).

Suppose that we have a linear system $\bm{A}\bm{x}=\bm{b}$, where
$\bm{A}\in\mathbb{R}^{N\times N}$, $\bm{x}\in\mathbb{R}^{N}$, and
$\bm{b}\in\mathbb{R}^{N}$. Given $\bm{A}^{\prime}=\bm{A}+\bm{E}$ and
$\bm{b}^{\prime}=\bm{b}+\bm{\epsilon}$, where $\bm{E}$ and $\bm{\epsilon}$ are
noise in the system, we denote
$\bm{x}^{\prime}=\bm{A}^{\prime-1}\bm{b}^{\prime}=\bm{x}+\bm{\epsilon}_{x}$,
where the relative error of $\bm{x}$, $i.e.$,
$\frac{\|\bm{\epsilon}_{x}\|_{2}}{\|\bm{x}\|_{2}}$, is bounded as

$\displaystyle\frac{\|\bm{\epsilon}_{x}\|_{2}}{\|\bm{x}\|_{2}}\leq\tau_{\bm{A}}\left(\frac{\|\bm{E}\|_{2}}{\|\bm{A}\|_{2}}+\frac{\|\bm{\epsilon}\|_{2}}{\|\bm{\mu}\|_{2}}\right),$
(30)

where $\tau_{\bm{A}}$ is the condition number of $\bm{A}$.

Based on the properties and the theorems above, we can proof Property 2.4 as
follows.

Property 2.4 _For $\text{HP}_{\mathcal{V}}(\bm{\mu}_{1},\bm{A}_{1})$ and
$\text{HP}_{\mathcal{U}}(\bm{\mu}_{2},\bm{A}_{2})\sim\text{GHP}_{\Omega}(f,g)$,
where $\text{GHP}_{\Omega}(f,g)$ satisfies Assumption 2.2 and
$V=|\mathcal{V}|\leq|\mathcal{U}|=U$, their average intensity vectors, $i.e.$,
$\bar{\bm{\lambda}}_{1}$ and $\bar{\bm{\lambda}}_{2}$, satisfy_

$\frac{d_{\text{w}}(\bar{\bm{\lambda}}_{1},\bar{\bm{\lambda}}_{2})}{\|\bar{\bm{\lambda}}_{1}\|_{2}}\leq\frac{1}{1-D\|\bm{A}_{1}\|_{2}}\left(\frac{\sqrt{2U}C^{g}}{C_{1}^{f}\|\bm{I}_{V}-D\bm{A}_{1}\|_{2}}+\frac{1}{\|\bm{\mu}_{1}\|_{2}}\right)\left(d_{\text{w}}(\bm{\mu}_{1},\bm{\mu}_{2})+\sqrt{\frac{U-V}{V}}\|\bm{\mu}_{1}\|_{2}\right)+\sqrt{\frac{U-V}{VU}}.$

###### Proof.

According to (4), our $\text{GHP}(f,g)$ model generates a Hawkes process
$\text{HP}_{\mathcal{V}}(\bm{\mu},\bm{A})$ from the samples in $\Omega$.
Denote $\bm{x}_{1}=\\{x_{v,1}\\}_{v=1}^{V}$ as the samples for
$\text{HP}_{\mathcal{V}}(\bm{\mu}_{1},\bm{A}_{1})$ and
$\bm{x}_{2}=\\{x_{u,2}\\}_{u=1}^{U}$ the samples for
$\text{HP}_{\mathcal{U}}(\bm{\mu}_{2},\bm{A}_{2})$. We have

$\displaystyle\begin{aligned}
&\mu_{v,1}=f(x_{v,1}),&&a_{vv^{\prime},1}=g(x_{v,1},x_{v^{\prime},1})&&\forall~{}x_{v,1}\in\bm{x}_{1},~{}\mu_{v,1}\in\bm{\mu}_{1},~{}a_{vv^{\prime},1}\in\bm{A}_{1},\\\
&\mu_{u,2}=f(x_{u,2}),&&a_{uu^{\prime},2}=g(x_{u,2},x_{u^{\prime},2})&&\forall~{}x_{u,2}\in\bm{x}_{2},~{}\mu_{u,2}\in\bm{\mu}_{2},~{}a_{uu^{\prime},2}\in\bm{A}_{2}.\end{aligned}$
(31)

Because $V\leq U$, for $\text{HP}_{\mathcal{V}}(\bm{\mu}_{1},\bm{A}_{1})$ we
pad $\bm{\mu}_{1}$ and $\bm{A}_{1}$ with zeros, $i.e.$,
$\tilde{\bm{\mu}}_{1}=[\bm{\mu}_{1};\bm{0}_{U-V}]\in\mathbb{R}^{U}$ and
$\tilde{\bm{A}}_{1}=\begin{bmatrix}\bm{A}_{1}&\bm{0}_{V\times(U-V)}\\\
\bm{0}_{V\times(U-V)}^{\top}&\bm{0}_{(U-V)\times(U-V)}\end{bmatrix}\in\mathbb{R}^{U\times
U}$, such that $\tilde{\bm{\mu}_{1}}$ ($\tilde{\bm{A}_{1}}$) is as large as
$\bm{\mu}_{2}$ ($\bm{A}_{2}$). Accordingly, in $\Omega$ we denote the samples
corresponding to $\tilde{\bm{\mu}_{1}}$ and $\tilde{\bm{A}_{1}}$ as
$\tilde{\bm{x}}_{1}=\\{\tilde{x}_{v,1}\\}_{v=1}^{U}$, which is constructed by
padding $\bm{x}_{1}$ with $U-V$ zero points of $f(x)$, $i.e.$,

$\displaystyle\tilde{\bm{x}}_{1}=\\{\tilde{x}_{v,1}\\}_{v=1}^{U}=\\{x_{1,1},..,x_{V,1},\underbrace{x_{0}^{f},...,x_{0}^{f}}_{U-V}\\},$
(32)

where $x_{0}^{f}$ is the unique zero point of $f(x)$ (Assumption 2.2).

Because the Hawkes process generated by our GHP model is stationary (Property
2.1), according to (25) we have

$\displaystyle(\bm{I}_{V}-D\bm{A}_{1})\bar{\bm{\lambda}}_{1}=\bm{\mu}_{1},\quad(\bm{I}_{U}-D\bm{A}_{2})\bar{\bm{\lambda}}_{2}=\bm{\mu}_{2},\quad\text{and}\quad(\bm{I}_{U}-D\tilde{\bm{A}}_{1})\tilde{\bm{\lambda}}_{1}=\tilde{\bm{\mu}}_{1},$
(33)

where $\tilde{\bm{\lambda}}_{1}=[\bar{\bm{\lambda}}_{1};\bm{0}_{U-V}]$ is the
average intensity $\bar{\bm{\lambda}}_{1}$ with padded zeros.

Following the notations used in the proof of Property 2.3, we denote
$\bm{T}^{\mu}$ as the optimal transport matrix for
$d_{\text{w}}(\tilde{\bm{\mu}}_{1},\bm{\mu}_{2})$ and $\bm{T}^{\lambda}$ the
optimal transport matrix for
$d_{\text{w}}(\tilde{\bm{\lambda}}_{1},\bar{\bm{\lambda}}_{2})$. According to
Theorem A.5, these two matrices are normalized permutation matrices, $i.e.$,
$\bm{T}^{\mu}=\frac{1}{U}\bm{P}^{\mu}$ and
$\bm{T}^{\lambda}=\frac{1}{U}\bm{P}^{\lambda}$. Then, we have

$\displaystyle\begin{aligned}
d_{\text{w}}(\bar{\bm{\lambda}}_{1},\bar{\bm{\lambda}}_{2})&\leq
d_{\text{w}}(\bar{\bm{\lambda}}_{1},\tilde{\bm{\lambda}}_{1})+d_{\text{w}}(\tilde{\bm{\lambda}}_{1},\bar{\bm{\lambda}}_{2})&(\text{Property~{}\ref{prop:tri}})\\\
&\leq\sqrt{\frac{U-V}{VU}}\|\bar{\bm{\lambda}}_{1}\|_{2}+\frac{1}{\sqrt{U}}\|\tilde{\bm{\lambda}}_{1}-\bm{P}^{\lambda}\bar{\bm{\lambda}}_{2}\|_{2}.&(\text{Theorem~{}\ref{the:padzero}
+ Theorem~{}\ref{the:emd1d}})\\\
&\leq\sqrt{\frac{U-V}{VU}}\|\bar{\bm{\lambda}}_{1}\|_{2}+\frac{1}{\sqrt{U}}\|\tilde{\bm{\lambda}}_{1}-\bm{P}^{\mu}\bar{\bm{\lambda}}_{2}\|_{2}.&(\text{Based
on}~{}(\ref{eq:cond_t}))\end{aligned}$ (34)

Because the permutation matrix $\bm{P}^{\mu}$ satisfies
$\bm{P}^{\mu}(\bm{P}^{\mu})^{\top}=\bm{I}_{U}$, we have

$\displaystyle\bm{P}^{\mu}\bm{\mu}_{2}=\bm{P}^{\mu}(\bm{I}_{U}-D\bm{A}_{2})\bar{\bm{\lambda}}_{2}=\bm{P}^{\mu}(\bm{I}_{U}-D\bm{A}_{2})(\bm{P}^{\mu})^{\top}\bm{P}^{\mu}\bar{\bm{\lambda}}_{2}.$
(35)

We can treat (35) as a perturbed version of the linear system
$(\bm{I}_{U}-D\tilde{\bm{A}}_{1})\tilde{\bm{\lambda}}_{1}=\tilde{\bm{\mu}}_{1}$
and obtain

$\displaystyle\begin{aligned}
&\|\tilde{\bm{\lambda}}_{1}-\bm{P}^{\mu}\bar{\bm{\lambda}}_{2}\|_{2}&\\\
&\leq\|\tilde{\bm{\lambda}}_{1}\|_{2}\tau_{\bm{I}_{U}-D\tilde{\bm{A}}_{1}}\left(\frac{\|(\bm{I}_{U}-D\tilde{\bm{A}}_{1})-\bm{P}^{\mu}(\bm{I}_{U}-D\bm{A}_{2})(\bm{P}^{\mu})^{\top}\|_{2}}{\|\bm{I}_{U}-D\tilde{\bm{A}}_{1}\|_{2}}+\frac{\|\tilde{\bm{\mu}}_{1}-\bm{P}^{\mu}\bm{\mu}_{2}\|_{2}}{\|\tilde{\bm{\mu}}_{1}\|_{2}}\right)&(\text{Theorem}~{}\ref{the:noise})\\\
&=\|\tilde{\bm{\lambda}}_{1}\|_{2}\tau_{\bm{I}_{U}-D\tilde{\bm{A}}_{1}}\left(\frac{D\|\tilde{\bm{A}}_{1}-\bm{P}^{\mu}\bm{A}_{2}(\bm{P}^{\mu})^{\top}\|_{2}}{\|\bm{I}_{U}-D\tilde{\bm{A}}_{1}\|_{2}}+\frac{\sqrt{U}d_{\text{w}}(\tilde{\bm{\mu}}_{1},\bm{\mu}_{2})}{\|\tilde{\bm{\mu}}_{1}\|_{2}}\right)&(\text{Theorem}~{}\ref{the:emd1d})\\\
&\leq\|\bar{\bm{\lambda}}_{1}\|_{2}\frac{1}{1-D\|\bm{A}_{1}\|_{2}}\left(\frac{D\|\tilde{\bm{A}}_{1}-\bm{P}^{\mu}\bm{A}_{2}(\bm{P}^{\mu})^{\top}\|_{2}}{\|\bm{I}_{V}-D\bm{A}_{1}\|_{2}}+\frac{\sqrt{U}d_{\text{w}}(\tilde{\bm{\mu}}_{1},\bm{\mu}_{2})}{\|\bm{\mu}_{1}\|_{2}}\right).\end{aligned}$
(36)

The second inequality in (36) is because 1)
$\|\tilde{\bm{\lambda}}_{1}\|_{2}=\|\bar{\bm{\lambda}}_{1}\|_{2}$; 2)
$\|\tilde{\bm{\mu}}_{1}\|_{2}=\|\bm{\mu}_{1}\|_{2}$; 3)
$\tau_{\bm{I}_{U}-D\tilde{\bm{A}}_{1}}=\frac{\sigma_{\max}(\bm{I}_{U}-D\tilde{\bm{A}}_{1})}{\sigma_{\min}(\bm{I}_{U}-D\tilde{\bm{A}}_{1})}=\frac{1-D\sigma_{\min}(\tilde{\bm{A}}_{1})}{1-D\sigma_{\max}(\tilde{\bm{A}}_{1})}\leq\frac{1}{1-D\|\tilde{\bm{A}}_{1}\|_{2}}=\frac{1}{1-D\|\bm{A}_{1}\|_{2}}$;
and 4)
$\|\bm{I}_{U}-D\tilde{\bm{A}}_{1}\|_{2}=1-D\sigma_{\min}(\tilde{\bm{A}}_{1})\geq
1-D\sigma_{\min}(\bm{A}_{1})=\|\bm{I}_{V}-D\bm{A}_{1}\|_{2}$, where
$\sigma_{\min}$ ($\sigma_{\max}$) represents the minimum (the maximum)
eigenvalue of a matrix.

For the
$\|\tilde{\bm{A}}_{1}-\bm{P}^{\mu}\bm{A}_{2}(\bm{P}^{\mu})^{\top}\|_{2}$ in
(36), we have

$\displaystyle\begin{aligned}
&\|\tilde{\bm{A}}_{1}-\bm{P}^{\mu}\bm{A}_{2}(\bm{P}^{\mu})^{\top}\|_{2}&\\\
&=\Bigl{(}\sideset{}{{}_{v,v^{\prime}=1}^{U}}{\sum}\sideset{}{{}_{u,u^{\prime}=1}^{U}}{\sum}|\tilde{a}_{vv^{\prime},1}-a_{uu^{\prime},2}|^{2}p_{vu}^{\mu}p_{v^{\prime}u^{\prime}}^{\mu}\Bigr{)}^{\frac{1}{2}}&\\\
&\leq
C^{g}\Bigl{(}\sideset{}{{}_{v,v^{\prime}=1}^{U}}{\sum}\sideset{}{{}_{u,u^{\prime}=1}^{U}}{\sum}\|[\tilde{x}_{v,1};\tilde{x}_{v^{\prime},1}]-[x_{u,2};x_{u^{\prime},2}]\|_{2}^{2}p_{vu}^{\mu}p_{v^{\prime}u^{\prime}}^{\mu}\Bigr{)}^{\frac{1}{2}}&(\text{Property}~{}\ref{prop:lip})\\\
&=C^{g}\Bigl{(}\sideset{}{{}_{v,v^{\prime}=1}^{U}}{\sum}\sideset{}{{}_{u,u^{\prime}=1}^{U}}{\sum}(|\tilde{x}_{v,1}-x_{u,2}|^{2}+|\tilde{x}_{v^{\prime},1}-x_{u^{\prime},2}|^{2})p_{vu}^{\mu}p_{v^{\prime}u^{\prime}}^{\mu}\Bigr{)}^{\frac{1}{2}}&\\\
&=C^{g}\left(\sum_{v,u=1}^{U}p_{vu}\sum_{v^{\prime},u^{\prime}=1}^{U}|\tilde{x}_{v^{\prime},1}-x_{u^{\prime},2}|^{2}p_{v^{\prime}u^{\prime}}^{\mu}+C^{g}\sum_{v^{\prime},u^{\prime}=1}^{N}p_{v^{\prime}u^{\prime}}\sum_{v,u=1}^{N}|\tilde{x}_{v,1}-x_{u,2}|^{2}p_{vu}^{\mu}\right)^{\frac{1}{2}}&\\\
&=\sqrt{2}C^{g}\Bigl{(}U\sideset{}{{}_{v,u=1}^{U}}{\sum}|\tilde{x}_{v,1}-x_{u,2}|^{2}p_{vu}^{\mu}\Bigr{)}^{\frac{1}{2}}&\\\
&\leq\frac{\sqrt{2}C^{g}}{C_{1}^{f}}\Bigl{(}U\sideset{}{{}_{v,u=1}^{U}}{\sum}|\tilde{\mu}_{v,1}-\mu_{u,2}|^{2}p_{vu}^{\mu}\Bigr{)}^{\frac{1}{2}}&(\text{Property}~{}\ref{prop:lip})\\\
&=\frac{\sqrt{2U}C^{g}}{C_{1}^{f}}\|\tilde{\bm{\mu}}_{1}-\bm{P}^{\mu}\bm{\mu}_{2}\|_{2}&(\text{Theorem}~{}\ref{the:emd1d})\\\
&=\frac{\sqrt{2}UC^{g}}{C_{1}^{f}}d_{\text{w}}(\tilde{\bm{\mu}}_{1},\bm{\mu}_{2}).&\end{aligned}$
(37)

Plugging (37) into (36), we have

$\displaystyle\begin{aligned}
&\frac{\|\tilde{\bm{\lambda}}_{1}-\bm{P}^{\mu}\bar{\bm{\lambda}}_{2}\|_{2}}{\|\bar{\bm{\lambda}}_{1}\|_{2}}&\\\
&\leq\frac{d_{\text{w}}(\tilde{\bm{\mu}}_{1},\bm{\mu}_{2})}{1-D\|\bm{A}_{1}\|_{2}}\left(\frac{\sqrt{2}UC^{g}}{C_{1}^{f}\|\bm{I}_{V}-K\bm{A}_{1}\|_{2}}+\frac{\sqrt{U}}{\|\bm{\mu}_{1}\|_{2}}\right)&\\\
&\leq\frac{1}{1-D\|\bm{A}_{1}\|_{2}}\left(\frac{\sqrt{2}UC^{g}}{C_{1}^{f}\|\bm{I}_{V}-D\bm{A}_{1}\|_{2}}+\frac{\sqrt{U}}{\|\bm{\mu}_{1}\|_{2}}\right)(d_{\text{w}}(\tilde{\bm{\mu}}_{1},\bm{\mu}_{1})+d_{\text{w}}(\bm{\mu}_{1},\bm{\mu}_{2}))&(\text{Property}~{}\ref{prop:tri})\\\
&\leq\frac{1}{1-D\|\bm{A}_{1}\|_{2}}\left(\frac{\sqrt{2}UC^{g}}{C_{1}^{f}\|\bm{I}_{V}-D\bm{A}_{1}\|_{2}}+\frac{\sqrt{U}}{\|\bm{\mu}_{1}\|_{2}}\right)\left(\sqrt{\frac{U-V}{V}}\|\bm{\mu}_{1}\|_{2}+d_{\text{w}}(\bm{\mu}_{1},\bm{\mu}_{2})\right).&(\text{Theorem}~{}\ref{the:padzero})\end{aligned}$
(38)

Finally, plugging (38) into (34), we have

$\displaystyle\begin{aligned}
\frac{d_{\text{w}}(\bar{\bm{\lambda}}_{1},\bar{\bm{\lambda}}_{2})}{\|\bar{\bm{\lambda}}_{1}\|_{2}}&\leq\sqrt{\frac{U-V}{VU}}+\frac{1}{1-D\|\bm{A}_{1}\|_{2}}\left(\frac{\sqrt{2U}C^{g}}{C_{1}^{f}\|\bm{I}_{V}-K\bm{A}_{1}\|_{2}}+\frac{1}{\|\bm{\mu}_{1}\|_{2}}\right)\left(\sqrt{\frac{U-V}{V}}\|\bm{\mu}_{1}\|_{2}+d_{\text{w}}(\bm{\mu}_{1},\bm{\mu}_{2})\right).\end{aligned}$
(39)

∎

###  A.4 The proof of Corollary 2.5

Plugging $U=V$ into (39), we obtain

$\displaystyle\begin{aligned}
\frac{d_{\text{w}}(\bar{\bm{\lambda}}_{1},\bar{\bm{\lambda}}_{2})}{\|\bar{\bm{\lambda}}_{1}\|_{2}}&\leq\frac{d_{\text{w}}(\bm{\mu}_{1},\bm{\mu}_{2})}{1-D\|\bm{A}_{1}\|_{2}}\left(\frac{\sqrt{2V}C^{g}}{C_{1}^{f}\|\bm{I}_{V}-D\bm{A}_{1}\|_{2}}+\frac{1}{\|\bm{\mu}_{1}\|_{2}}\right).\end{aligned}$
(40)

##  Appendix B Details of Our Algorithm and Experiments

###  B.1 The implementation of $d(\widehat{N}_{u}^{k},N_{v}^{l})$ in (10)

Denote $\\{t_{1}^{u},...,t_{I}^{u}\\}\subset[0,T]$ as the sequence
corresponding to $\widehat{N}_{u}^{k}(t)$ and
$\\{t_{1}^{v},...,t_{J}^{v}\\}\subset[0,T]$ the sequence corresponding to
$N_{v}^{l}(t)$. Without loss of generality, we assume $I\leq J$ and calculate
$d(\widehat{N}_{u}^{k},N_{v}^{l})$ as

$\displaystyle\begin{aligned}
d(\widehat{N}_{u}^{k},N_{v}^{l})=\frac{1}{T}\int_{0}^{T}|\widehat{N}_{u}^{k}(t)-N_{v}^{l}(t)|\text{d}t=\frac{1}{T}\sideset{}{{}_{i=1}^{I}}{\sum}|t_{i}^{u}-t_{i}^{v}|+\sideset{}{{}_{i=I+1}^{J}}{\sum}|T-t_{i}^{v}|.\end{aligned}$
(41)

As proven in (Xiao et al., 2017), the distance in (41) is a valid metric for
the event sequences with a single event type.

###  B.2 The significance of $\bm{V}_{\max}$ in theory

Suppose that we have two GHP models, whose maximum numbers of event types are
$V_{\max}$ and $\hat{V}_{\max}$, respectively. Based on (4), we know that the
expected numbers of event types of their event sequences are
$\frac{V_{\max}}{2}$ and $\frac{\hat{V}_{\max}}{2}$. For the two sequences
having $\frac{\hat{V}_{\max}}{2}$ and $\frac{V_{\max}}{2}$ event types,
respectively, Property 2.4 indicates that the difference between their average
intensity vectors is
$\mathcal{O}\left(\sqrt{\frac{|V_{\max}-\hat{V}_{\max}|}{\min\\{V_{\max},\hat{V}_{\max}\\}}}\right)$.
Therefore, when training our GHP model, we need to carefully set
$\bm{V}_{\max}$ based on the training data. Empirically, we calculate the
averaged number of event types per sequence, denoted as $\bar{V}$ and set
$\bm{V}_{\max}=2\bar{V}$.

###  B.3 The Sinkhorn scaling algorithm

When calculating the HOT distance, we need to solve a series of optimal
transport problem. All the problems can be written in the following matrix
format:

$\displaystyle\min_{\bm{T}\in\Pi(\bm{p},\bm{q})}\langle\bm{D},\bm{T}\rangle.$
(42)

We apply the Sinkhorn scaling algorithm (Cuturi, 2013) to solve this problem
approximately. In particular, we add an entropic regularizer with a weight
$\beta$ into (42) and rewrite it as

$\displaystyle\min_{\bm{T}\in\Pi(\bm{p},\bm{q})}\langle\bm{D},\bm{T}\rangle+\beta\langle\log\bm{T},\bm{T}\rangle.$
(43)

Then, we can solve (43) by the following algorithm.

Algorithm 2
$\min_{\bm{T}\in\Pi(\bm{p},\bm{q})}\langle\bm{D},\bm{T}\rangle+\beta\langle\log\bm{T},\bm{T}\rangle$

1: Initialize $\bm{T}^{(0)}=\bm{p}\bm{q}^{\top}$, $\bm{a}=\bm{p}$,
$\bm{C}=\exp(-\frac{1}{\beta}\bm{D})$.

2: for $j=0,...,J-1$

3: Sinkhorn iteration: $\bm{b}=\frac{\bm{q}}{\bm{C}^{\top}\bm{a}}$,
$\bm{a}=\frac{\bm{p}}{\bm{C}\bm{b}}$,

4: $\bm{T}=\text{diag}(\bm{a}){C}\text{diag}(\bm{b})$.

###  B.4 The implementation of $d_{\text{fgw}}(\hat{\theta},\theta)$

Given two GHP modeled on $\Omega$, denoted as
$\text{GHP}_{\Omega}(f_{a},g_{a})$ and $\text{GHP}_{\Omega}(f_{b},g_{b})$.
Denote $\mathcal{S}_{\Omega}$ as the set of measure-preserving mappings from
$\Omega$ to $\Omega$. Based on the theory of graphon (Lovász, 2012), the
distance between these two GHP models can be measured by the $\delta_{2}$
distance between their parameters:

$\displaystyle
d(\text{GHP}_{\Omega}(f_{a},g_{a}),\text{GHP}_{\Omega}(f_{b},g_{b}))=\inf_{s\in\mathcal{S}_{\Omega}}\|f_{a}-f_{b}^{s}\|_{L_{2}}+\|g_{a}-g_{b}^{s}\|_{L_{2}},$
(44)

where $f_{b}^{s}(x)=f_{b}(s(x))$ and $g_{b}^{s}(x,y)=g_{b}(s(x),s(y))$. We say
the two GHP models are equivalent if there exists at least one map
$s\in\mathcal{S}_{\Omega}$ making
$d(\text{GHP}_{\Omega}(f_{a},g_{a}),\text{GHP}_{\Omega}(f_{b},g_{b}))=0$.

According to the theory of optimal transport (Villani, 2008), the first term
of (44) can be implemented as the Wasserstein distance between $f_{a}$ and
$f_{b}$. For the second term, the work in (Lovász, 2012) implies that we can
rewrite them as the Gromov-Wasserstein distance (Mémoli, 2011) between $g_{a}$
and $g_{b}$. Combining these two distance together leads to the Fused Gromov-
Wasseserstein (FGW) distance (Vayer et al., 2018):

$\displaystyle\begin{aligned}
&d_{\text{fgw}}(\text{GHP}_{\Omega}(f_{a},g_{a}),\text{GHP}_{\Omega}(f_{b},g_{b}))\\\
:=&\sideset{}{{}_{\pi\in\Pi(p,q)}}{\inf}\mathbb{E}_{x,x^{\prime}\sim\pi}\left[|f_{a}(x)-f_{b}(x^{\prime})|^{2}\right]+\mathbb{E}_{x,x^{\prime},y,y^{\prime}\sim\pi\otimes\pi}\left[|g_{a}(x,y)-g_{b}(x^{\prime},y^{\prime})|^{2}\right].\end{aligned}$
(45)

Here, we assume $p$ and $q$ are two uniform distribution on $\Omega$. In our
experiment, $\text{GHP}_{\Omega}(f_{a},g_{a})$ and
$\text{GHP}_{\Omega}(f_{b},g_{b})$ correspond to the ground truth model and
the learning result, respectively, and we use (45) as the measurement of the
estimation error.

In practice, we set $\Omega=[0,1]$ and uniformly $N$ samples from it, $i.e.$,
$\left\\{0,\frac{1}{N},...,\frac{N-1}{N}\right\\}$. Accordingly, we obtain the
discrete representation of each function, $i.e.$,
$\bm{f}_{a}=[f_{i}^{a}]\in\bm{R}^{N}$, $\bm{f}_{b}=[f_{i}^{b}]\in\bm{R}^{N}$,
$\bm{G}_{a}=[g_{ij}^{a}]\in\bm{R}^{N\times N}$, and
$\bm{G}_{b}=[g_{ij}^{b}]\in\bm{R}^{N\times N}$. Then, we obtain the discrete
version of (45)

$\displaystyle\begin{aligned}
&\sideset{}{{}_{\bm{T}\in\Pi(\bm{p},\bm{q})}}{\min}\sum_{i,j=1}^{N}T_{ij}|f_{i}^{a}-f_{j}^{b}|^{2}+\sum_{i,i^{\prime},j,j^{\prime}=1}^{N}T_{ij}T_{i^{\prime}j^{\prime}}|g_{ii^{\prime}}^{a}-g_{jj^{\prime}}^{b}|^{2}\\\
=&\sideset{}{{}_{\bm{T}\in\Pi(\bm{p},\bm{q})}}{\min}\langle\bm{D}_{f},\bm{T}\rangle+\langle\bm{D}_{g}-2\bm{G}_{a}\bm{T}\bm{G}_{b}^{\top},\bm{T}\rangle,\end{aligned}$
(46)

where $\bm{D}_{f}=[|f_{i}^{a}-f_{j}^{b}|^{2}]$,
$\bm{D}_{g}=\frac{1}{N}(\bm{G}_{a}\odot\bm{G}_{a}+\bm{G}_{b}\odot\bm{G}_{b})$,
and $\odot$ is the Hadamard product. This problem can be solved by the
proximal gradient method in (Xu et al., 2019).

Algorithm 3
$\min_{\bm{T}\in\Pi(\bm{p},\bm{q})}\langle\bm{D}_{f},\bm{T}\rangle+\langle\bm{D}_{g}-2\bm{G}_{a}\bm{T}\bm{G}_{b}^{\top},\bm{T}\rangle$

1: Initialize $\bm{T}^{(0)}=\bm{p}\bm{q}^{\top}$, $\bm{a}=\bm{p}$

2: for $j=0,...,J-1$

3:
$\bm{C}=\exp\left(-\frac{1}{\alpha}(\bm{D}_{f}+\bm{D}_{g}-2\bm{G}_{a}\bm{T}^{(j)}\bm{G}_{b}^{\top})\right)\odot\bm{T}^{(j)}$.

4: Sinkhorn iteration: $\bm{b}=\frac{\bm{q}}{\bm{C}^{\top}\bm{a}}$,
$\bm{a}=\frac{\bm{p}}{\bm{C}\bm{b}}$,

5: $\bm{T}^{(j+1)}=\text{diag}(\bm{a})\bm{C}\text{diag}(\bm{b})$.

6: Return $\bm{T}^{(J)}$

###  B.5 The $d_{ot}(\widehat{\mathcal{N}},\mathcal{N})$ of baselines

The baselines also calculate $d_{ot}(\widehat{\mathcal{N}},\mathcal{N})$ by
(9). However, because the event types of their generated sequences perfectly
correspond to those of the testing sequences, they can calculate the distance
between each pair of sequences as
$d(\hat{\bm{N}}_{k},\bm{N}_{l})=\frac{1}{|\mathcal{V}|}\sum_{v\in\mathcal{V}}d(\widehat{N}_{v}^{k},N_{v}^{l})$
rather than using (10).