Variational Inference for Sparse Gaussian Process Modulated Hawkes Process

TL;DR

This paper introduces a scalable variational inference method for a non-parametric Bayesian Hawkes process using a sparse Gaussian process to model triggering kernels, improving efficiency and model fitting.

Contribution

It proposes a novel variational inference schema with a tighter ELBO and linear time complexity for Bayesian Hawkes processes, enhancing scalability and accuracy.

Findings

Outperforms state-of-the-art methods on synthetic and social media data.

Achieves linear time inference leveraging stationarity.

Tighter ELBO improves model selection accuracy.

Abstract

The Hawkes process (HP) has been widely applied to modeling self-exciting events including neuron spikes, earthquakes and tweets. To avoid designing parametric triggering kernel and to be able to quantify the prediction confidence, the non-parametric Bayesian HP has been proposed. However, the inference of such models suffers from unscalability or slow convergence. In this paper, we aim to solve both problems. Specifically, first, we propose a new non-parametric Bayesian HP in which the triggering kernel is modeled as a squared sparse Gaussian process. Then, we propose a novel variational inference schema for model optimization. We employ the branching structure of the HP so that maximization of evidence lower bound (ELBO) is tractable by the expectation-maximization algorithm. We propose a tighter ELBO which improves the fitting performance. Further, we accelerate the novel variational inference schema to linear time complexity by leveraging the stationarity of the triggering kernel. Different from prior acceleration methods, ours enjoys higher efficiency. Finally, we exploit synthetic data and two large social media datasets to evaluate our method. We show that our approach outperforms state-of-the-art non-parametric frequentist and Bayesian methods. We validate the efficiency of our accelerated variational inference schema and practical utility of our tighter ELBO for model selection. We observe that the tighter ELBO exceeds the common one in model selection.

Tables

Table 1. Table 1: Notations

VBHP	VBPP	VI
$𝒟\equiv {\left\{{𝒙}_n\right\}}_{n=1}^N$	$𝒟\equiv {\left\{{𝒙}_n\right\}}_{n=1}^N$	$𝒙$
$B,\mu ,f,𝒖$	$f,𝒖$	$𝒛$
$k_0,c_0,{\left\{{\alpha }_i\right\}}_{i=1}^R,\gamma$	${\left\{{\alpha }_i\right\}}_{i=1}^R,\gamma$	$𝜽$
$k,c,𝒎,𝑺,{\left\{{\alpha }_i\right\}}_{i=1}^R,\gamma ,$	$𝒎,𝑺,{\left\{{\alpha }_i\right\}}_{i=1}^R,\gamma$	${𝜽}^{\prime }$
${\left\{{\left\{q_{ij}\right\}}_{j=0}^{i-1}\right\}}_{i=1}^N$

Table 2. Table 2: Results on Synthetic and Real World Data (Mean ± plus-or-minus \pm One Standard Variance). VBHP (C) and (T) use the CELBO and the TELBO to update the hyper-parameters respectively.

Measure	Data	SumExp	ODE	WH	Gibbs Hawkes	VBHP (C)	VBHP (T)
${\text{L}}_2$	Sin	$\varphi$:${0.693}_{\pm 0.028}$	${0.665}_{\pm 0.121}$	${2.463}_{\pm 0.145}$	${0.408}_{\pm 0.198}$	${0.152}_{\pm 0.091}$	${0.183}_{\pm 0.076}$
		$\mu$:${2.968}_{\pm 1.640}$	${4.514}_{\pm 3.808}$	${6.794}_{\pm 5.054}$	${4.108}_{\pm 3.949}$	${0.640}_{\pm 0.528}$	${0.579}_{\pm 0.523}$
	Cos	$\varphi$:${0.473}_{\pm 0.102}$	${0.697}_{\pm 0.065}$	${1.743}_{\pm 0.083}$	${0.667}_{\pm 0.686}$	${0.325}_{\pm 0.073}$	${0.292}_{\pm 0.096}$
		$\mu$:${2.751}_{\pm 1.902}$	${7.030}_{\pm 5.662}$	${6.099}_{\pm 4.613}$	${4.685}_{\pm 4.421}$	${0.555}_{\pm 0.294}$	${0.515}_{\pm 0.293}$
	Exp	$\varphi$:${0.133}_{\pm 0.138}$	${1.835}_{\pm 0.539}$	${2.254}_{\pm 2.042}$	${0.676}_{\pm 0.233}$	${0.257}_{\pm 0.086}$	${0.235}_{\pm 0.102}$
		$\mu$:${3.290}_{\pm 1.991}$	${8.969}_{\pm 8.604}$	${16.66}_{\pm 20.95}$	${7.648}_{\pm 9.647}$	${0.471}_{\pm 0.432}$	${0.486}_{\pm 0.418}$
HLL	Sin	${3.490}_{\pm 0.400}$	${3.489}_{\pm 0.413}$	${3.233}_{\pm 0.273}$	${3.492}_{\pm 0.406}$	${3.488}_{\pm 0.400}$	${3.497}_{\pm 0.406}$
	Cos	${3.874}_{\pm 0.544}$	${3.872}_{\pm 0.552}$	${3.613}_{\pm 0.373}$	${3.871}_{\pm 0.562}$	${3.876}_{\pm 0.541}$	${3.878}_{\pm 0.548}$
	Exp	${2.825}_{\pm 0.481}$	${2.822}_{\pm 0.496}$	${2.782}_{\pm 0.490}$	${2.826}_{\pm 0.492}$	${2.826}_{\pm 0.491}$	${2.829}_{\pm 0.487}$
	ACTIVE	${1.692}_{\pm 1.371}$	${0.880}_{\pm 2.716}$	${0.710}_{\pm 0.943}$	${1.323}_{\pm 2.160}$	${1.824}_{\pm 1.159}$	${1.867}_{\pm 1.181}$
	SEISMIC	${2.943}_{\pm 0.959}$	${2.582}_{\pm 1.665}$	${1.489}_{\pm 1.796}$	${3.110}_{\pm 1.251}$	${3.143}_{\pm 0.895}$	${3.164}_{\pm 0.843}$

Equations

$$\lambda (\mathbf{x})=\mu +\sum _{{\mathbf{x}}_i<\mathbf{x}}\varphi (\mathbf{x}-{\mathbf{x}}_i),$$

$$\log p(\mathcal{D},B|\mu ,\varphi )=$$

$$\displaystyle\sum_{i=1}^{N}\Big{(}\sum_{j=1}^{i-1}b_{ij}\log\phi_{ij}+b_{i0}\log\mu\Big{)}-\sum_{i=1}^{N}\int_{\mathcal{T}_{i}}\phi-\mu|\mathcal{T}|,~{}~{}$$

$$\log p(\mathbf{x}|\mathbf{\theta })$$

$$=\log \frac{p(\mathbf{x},\mathbf{z}|\mathbf{\theta })}{q(\mathbf{z}|{\mathbf{\theta }}^{\prime })}-\log \frac{p(\mathbf{z}|\mathbf{x},\mathbf{\theta })}{q(\mathbf{z}|{\mathbf{\theta }}^{\prime })}$$

$$\displaystyle=\underbrace{\underbrace{\operatorname{\mathbb{E}}_{q(\bm{z}|\bm{\theta}^{\prime})}\big{[}\log p(\bm{x}|\bm{z},\bm{\theta})\big{]}}_{\text{reconstruction term}}-\underbrace{\text{KL}(q(\bm{z}|\bm{\theta}^{\prime})||p(\bm{z}|\bm{\theta}))}_{\text{regularisation term}}}_{\equiv\text{CELBO}(q(\bm{z}),p(\bm{x}|\bm{z}),p(\bm{z}))}$$

$$+\underset{\begin{array}{c}\text{intractable (non-negative) gap}\end{array}}{\underbrace{\text{KL}(q(\mathbf{z}|{\mathbf{\theta }}^{\prime })||p(\mathbf{z}|\mathbf{x},\mathbf{\theta }))}}$$

$$\ge \text{CELBO}(q(\mathbf{z}),p(\mathbf{x}|\mathbf{z}),p(\mathbf{z})),$$

$$K(\bm{x},\bm{x}^{\prime})\equiv\gamma\prod_{r=1}^{R}\exp\Big{(}-\dfrac{(x_{r}-x_{r}^{\prime})^{2}}{2\alpha_{r}}\Big{)}.$$

$$q(f,\mathbf{u})\equiv p(f|\mathbf{u})q(\mathbf{u})$$

$$q(f)$$

$$\nu (\mathbf{x})$$

$$\Sigma (\mathbf{x},{\mathbf{x}}^{\prime })$$

$$\text{CELBO}(q(f,\mathbf{u}),p(\mathcal{D}|f,\mathbf{u}),p(f,\mathbf{u}))$$

$$=E_{q(f)}[\log p(\mathcal{D}|f)]-\text{KL}(q(\mathbf{u})||p(\mathbf{u})).$$

$$q(B)=\prod _{i=1}^N\prod _{j=0}^{i-1}{q}_{ij}^{b_{ij}}.$$

$$q(B,\mu ,f,\mathbf{u})\equiv q(B)q(\mu )p(f|\mathbf{u})q(\mathbf{u}),$$

$$\text{CELBO}(q(B,\mu ,f,\mathbf{u}),p(\mathcal{D}|B,\mu ,f,\mathbf{u}),p(B,\mu ,f,\mathbf{u}))$$

$$\displaystyle=\underbrace{\operatorname{\mathbb{E}}_{q(B,\mu,f)}\Big{[}\log p(\mathcal{D},B|f,\mu)\Big{]}}_{\text{Data Dependent Expectation (DDE)}}+H_{B}$$

$$-\text{KL}(q(\mu )||p(\mu ))-\text{KL}(q(\mathbf{u})||p(\mathbf{u})).$$

$$E_{q(B)}[\log p(\mathcal{D},B|f,\mu )]$$

$$\displaystyle=\sum_{i=1}^{N}\Big{(}\sum_{j=1}^{i-1}q_{ij}\log f^{2}_{ij}+q_{i0}\log\mu\Big{)}-\sum_{i=1}^{N}\int_{\mathcal{T}_{i}}f^{2}-\mu|\mathcal{T}|$$

$$\displaystyle=\sum_{i=1}^{N}\Big{[}\sum_{j=1}^{i-1}q_{ij}\operatorname{\mathbb{E}}_{q(f)}(\log f_{ij}^{2})+q_{i0}\big{(}\psi(k)+\log c\big{)}$$

$$\displaystyle\quad-\int_{\mathcal{T}_{i}}\operatorname{\mathbb{E}}_{q(f)}^{2}(f)-\int_{\mathcal{T}_{i}}\operatorname{\text{Var}}_{q(f)}(f)\Big{]}-kc|\mathcal{T}|,$$

$$\underbrace{\operatorname{\mathbb{E}}_{q(B,\mu,f)}\Big{[}\log p(\mathcal{D},B|f,\mu)\Big{]}+H_{B}}_{\equiv\text{TELBO}}\leq\log p(\mathcal{D}).$$

$$\log p(\mathcal{D})$$

$$\underset{n\to \infty }{\lim }\frac{1}{n}\sum _{i=1}^n\log p(\mathcal{D},B_i)$$

$$=\sum _{B}q(B)\log p(\mathcal{D},B),$$

$$\sum _{i=1}^{\infty }\log p(\mathcal{D},B_i)=\log p(\{(\mathcal{D},B_i){\}}_{i=1}^{\infty }).$$

$$\log p(\{(\mathcal{D},B_i){\}}_{i=1}^n)$$

$$-\text{KL}$$

$$=\underset{n\to \infty }{\lim }\frac{1}{n}\log p(\{(\mathcal{D},B_i){\}}_{i=1}^n)$$

$$\displaystyle\overset{(a)}{\geq}\lim_{n\rightarrow\infty}\dfrac{1}{n}E_{q(f,\mu)}\big{[}\log p(\{(\mathcal{D},B_{i})\}_{i=1}^{n}|f,\mu)]$$

$$\displaystyle\quad\overset{(b)}{=}\lim_{n\rightarrow\infty}\sum_{n_{B}}\dfrac{n_{B}}{n}E_{q(f,\mu)}\big{[}\log p(\mathcal{D},B|f,\mu)\big{]}$$

$$\displaystyle\quad=E_{q(f,\mu,B)}\big{[}\log p(\mathcal{D},B|f,\mu)\big{]}$$

$$\{{\alpha }_i^{\ast }{\}}_{i=1}^R,{\gamma }^{\ast }={\text{argmax}}_{\{{\alpha }_i{\}}_{i=1}^R,\gamma }\text{TELBO}.$$

$$q_{ij}$$

$${\varphi }_{\sin }(x)=0.9[\sin (3x)+1],x\in [0,\pi /2];\text{otherwise},0;$$

$${\varphi }_{\cos }(x)=\cos (2x)+1,x\in [0,\pi /2];\text{otherwise},0;$$

Full text

1 Introduction

The Hawkes process (HP) (?) is particularly useful to model self-exciting point data – i.e., when the occurrence of a point increases the likelihood of occurrence of new points. The process is parameterized using a background intensity $\mu$, and a triggering kernel $\phi$. The Hawkes process can be alternatively represented as a cluster of Poisson processes (PPes) (?). In the cluster, a PP with an intensity $\mu$ (denoted as PP($\mu$)) generates immigrant points which are considered to arrive in the system from the outside, and every existing point triggers offspring points, which are generated internally through the self-excitement, following a PP($\phi$). Points can therefore be structured into clusters where each cluster contains either a point and its direct offspring or the background process (an example is shown in Fig. 1a). Connecting all points using the triggering relations yields a tree structure, which is called the branching structure (an example is shown in Fig. 1b corresponding to Fig.1a). With the branching structure, we can decompose the HP into a cluster of PPes. The triggering kernel $\phi$ is shared among all cluster Poisson processes relating to a HP, and it determines the overall behavior of the process. Consequently, designing the kernel functions is of utmost importance for employing the HP to a new application, and its study has attracted much attention.

Prior non-parametric frequentist solutions. In the case in which the optimal triggering kernel for a particular application is unknown, a typical solution is to express it using a non-parametric form, such as the work of ? (?); ? (?); ? (?); ? (?). These are all frequentist methods and among them, the Wiener-Hoef equation based method (?) enjoys linear time complexity. Our method has a similar advantage of linear time complexity per iteration. The Euler-Lagrange equation based solutions (?; ?) require discretizing the input domain so they face a problem of poorly scaling with the dimension of the domain. The same problem is also faced by ? (?)’s discretization based method. In contrast, our method requires no dixscretization so enjoys scalability with the dimension of the domain.

Prior Bayesian solutions. The Bayesian inference for the HP has also been studied, including the work of ? (?); ? (?); ? (?). These work require either constructing a parametric triggering kernel (?; ?) or discretizing the input domain to scale with the data size (?). The shortcoming of discretization is just mentioned and to overcome it, ? (?) propose a continuous non-parametric Bayesian HP and resort to an unscalable Markov chain Monte Carlo (MCMC) estimator to the posterior distribution. A more recent solution based on Gibbs sampling (?) obtains linear time complexity per iteration by exploiting the HP’s branching structure and the stationarity of the triggering kernel similar to the work of ? (?). This approach considers only high-probability triggering relationships in computations for acceleration. However, those relationships are updated in each iteration, which is less efficient than our pre-computing them. Besides, our variational inference schema enjoys faster convergence than that of the Gibbs sampling based method.

Contributions. In this paper, we propose the first sparse Gaussian process modulated HP which employs a novel variational inference schema, enjoys a linear time complexity per iteration and scales to large real world data. Our method is inspired by the variational Bayesian PP (VBPP) (?) which provides the Bayesian non-parametric inference only for the whole intensity of the HP without for its components: the background intensity $\mu$ and the triggering kernel $\phi$. Thus, the VBPP loses the internal reactions between points, and developing the variational Bayesian non-parametric inference for the HP is non-trivial and more challenging than the VBPP. In this paper, we adapt the VBPP for the HP and term the new approach the variational Bayesian Hawkes process (VBHP). The contributions are summarized:

(1) (Sec.3) We introduce a new Bayesian non-parametric HP which employs a sparse Gaussian process modulated triggering kernel and a Gamma distributed background intensity.

(2) (Sec.3&4) We propose a new variational inference schema for such a model. Specifically, we employ the branching structure of the HP so that maximization of the evidence lower bound (ELBO) is tractable by the expectation-maximization algorithm, and we contribute a tighter ELBO which improves the fitting performance of our model.

(3) (Sec.5) We propose a new acceleration trick based on the stationarity of the triggering kernel. The new trick enjoys higher efficiency than prior methods and accelerates the variational inference schema to linear time complexity per iteration.

(4) (Sec.6) We empirically show that VBHP provides more accurate predictions than state-of-the-art methods on synthetic data and on two large online diffusion datasets. We validate the linear time complexity and faster convergence of our accelerated variational inference schema compared to the Gibbs sampling method, and the practical utility of our tighter ELBO for model selection, which outperforms the common one in model selection.

2 Prerequisite

In this section, we review the Hawkes process, the variational inference and its application to the Poisson process.

2.1 Hawkes Process (HP)

The HP (?) is a self-exciting point process, in which the occurrence of a point increases the arrival rate $\lambda(\cdot)$, a.k.a. the (conditional) intensity, of new points. Given a set of time-ordered points $\mathcal{D}=\{\bm{x}_{i}\}_{i=1}^{N}$, $\bm{x}_{i}\in\mathbb{R}^{R}$, the intensity at $\bm{x}$ conditioned on given points is written as:

[TABLE]

where $\mu>0$ is a constant background intensity, and $\phi:\mathbb{R}^{R}\rightarrow[0,\infty)$ is the triggering kernel.

We are particularly interested in the branching structure of the HP. As introduced in Sec.1, each point $\bm{x}_{i}$ has a parent that we represent by a one-hot vector $\bm{b}_{i}=[b_{i0},b_{i1},\cdots,b_{i,i-1}]^{T}$: each element $b_{ij}$ is binary, $b_{ij}=1$ represents that $\bm{x}_{i}$ is triggered by $\bm{x}_{j}$ ($0\leq j\leq i-1$, $\bm{x}_{0}$: the background), and $\sum_{j=0}^{i-1}b_{ij}=1$. A branching structure $B$ is a set of $\bm{b}_{i}$, namely $B=\{\bm{b}_{i}\}_{i=1}^{N}$, and if the probability of $b_{ij}=1$ is defined as $p_{ij}\equiv p(b_{ij}=1)$, the probability of $B$ can be expressed as $p(B)=\prod_{i=1}^{N}\prod_{j=0}^{i-1}p_{ij}^{b_{ij}}$. Since $\bm{b}_{i}$ is a one-hot vector, there is $\sum_{j=0}^{i-1}p_{ij}=1$ for all $i$.

Given both $\mathcal{D}$ and a branching structure $B$, the log likelihood of $\mu$ and $\phi$ becomes a sum of log likelihoods of PPes. With $B$, the HP $\mathcal{D}$ can be decomposed into PP($\mu$) and $\{\text{PP}(\phi(\bm{x}-\bm{x}_{i}))\}_{i=1}^{N}$. The data domain of PP($\mu$) equals that of the HP, which we denote $\mathcal{T}$, and we denote the data domain of PP($\phi(\bm{x}-\bm{x}_{i})$) by $\mathcal{T}_{i}\subset\mathcal{T}$. As a result, the log likelihood $\log p(D,B|\mu,\phi)$ is expressed as:

[TABLE]

where $|\mathcal{T}|\equiv\int_{\mathcal{T}}1\text{ d}\bm{x}$, $\phi_{ij}\equiv\phi(\bm{x}_{i}-\bm{x}_{j})$ and $\int_{\mathcal{T}_{i}}\phi\equiv\int_{\mathcal{T}_{i}}\phi(\bm{x})\text{ d}\bm{x}$. Throughout this paper, we simplify the integration by omitting the integration variable, which is $\bm{x}$ unless otherwise specified. This work is developed based on the univariate HP. However, it can be extended to be multivariate following the same procedure.

2.2 Variational Inference (VI)

Consider the latent variable model $p(\bm{x},\bm{z}|\bm{\theta})$ where $\bm{x}$ and $\bm{z}$ are the data and the latent variables respectively. The variational approach introduces a variational distribution to approximate the posterior distribution $q(\bm{z}|\bm{\theta}^{\prime})\approx p(\bm{z}|\bm{x},\bm{\theta})$ and maximizes a lower bound of the log-likelihood, which can be derived from the non-negative gap perspective:

[TABLE]

where we omit $\bm{\theta}$ and $\bm{\theta}^{\prime}$ in conditions. We term this the Common ELBO (CELBO) to differentiate with our tighter ELBO of Sec.4. For notational convenience, we will often omit conditioning on $\bm{\theta}$ and $\bm{\theta}^{\prime}$ hereinafter. Optimizing the CELBO w.r.t. $\bm{\theta}^{\prime}$ balances between the reconstruction error and the Kullback-Leibler (KL) divergence from the prior. Generally, the conditional $p(\bm{x}|\bm{z})$ is known, so is the prior. Thus, for an appropriate choice of $q$, it is easier to work with this lower bound than with the intractable posterior $p(\bm{z}|\bm{x})$. We also see that, due to the form of the intractable gap, if $q$ is from a distribution family containing elements close to the true unknown posterior, then $q$ will be close to the true posterior when the CELBO is close to the true likelihood. An alternative derivation applies Jensen’s inequality (?).

2.3 Variational Bayesian Poisson Process (VBPP)

VBPP (?) applies the VI to the Bayesian Poisson process, which exploits the sparse Gaussian process (GP) to model the Poisson intensity. Specifically, VBPP uses a squared link function to map a sparse GP distributed function $f$ to the Poisson intensity $\lambda(\bm{x})=f^{2}(\bm{x})$. The sparse GP employs the ARD kernel:

[TABLE]

where $\gamma$ and $\{\alpha_{r}\}_{r=1}^{R}$ are GP hyper-parameters. Let $\bm{u}\equiv(f(\bm{z}_{1}),f(\bm{z}_{2}),\cdots,f(\bm{z}_{M}))$ where $\bm{z}_{i}$ are inducing points. The prior and the approximate posterior distributions of $\bm{u}$ are Gaussian distributions $p(\bm{u})=\mathcal{N}(\bm{u}|\bm{0},\bm{K}_{zz})$ and $q(\bm{u})=\mathcal{N}(\bm{u}|\bm{m},\bm{S})$ where $\bm{m}$ and $\mathbf{S}$ are the mean vector and the covariance matrix respectively. Note both $\bm{u}$ and $f$ employ zero mean priors. Notations of VBPP are connected with those of VI (Sec.2.2) in Table1.

Importantly, the variational joint distribution of $f$ and $\bm{u}$ uses the exact conditional distribution $p(f|\bm{u})$, i.e.,

[TABLE]

which in turn leads to the posterior GP:

[TABLE]

Then, the CELBO is obtained by using Eqn.(8):

[TABLE]

Note that the second term is the KL divergence between two multivariate Gaussian distributions, so is available in closed form. The first term turns out to be the expectation w.r.t. $q(f)$ of the log-likelihood $\log p(\mathcal{D}|f)=\sum_{i=1}^{N}\log f^{2}(\bm{x}_{i})-\int_{\mathcal{T}}f^{2}$. The expectation of the integral part is relatively straight-forward to compute and the expectation of the other (data-dependent) part is available in almost closed-form with a hyper-geometric function.

3 Variational Bayesian Hawkes Process

3.1 Notations

To extend VBPP to HP, we introduce two more variables: the background intensity $\mu$ and the branching structure $B$, defined in Sec.2.1. We assume that the prior distribution of $\mu$ is a Gamma distribution $p(\mu)=\text{Gamma}(\mu|k_{0},c_{0})$111$\text{Gamma}(\mu|k_{0},c_{0})=\dfrac{1}{\Gamma(k_{0})c_{0}^{k}}\mu^{k_{0}-1}e^{-x/c_{0}}$ and the posterior distribution is approximated by another Gamma distribution $q(\mu)=\text{Gamma}(\mu|k,c)$. For $B=\{\bm{b}_{i}\}_{i=1}^{N}$ given $\mathcal{D}=\{\bm{x}_{i}\}_{i=1}^{N}$, we assume the variational posterior $b_{ij}$, $j=0,\cdots,i-1$, have a categorical distribution: $q_{ij}=q(b_{ij}=1)$ and $\sum_{j=0}^{i-1}q_{ij}=1$, and thus, the variational posterior probability of $B$ is expressed as:

[TABLE]

The same squared link function is adopted for the triggering kernel $\phi(\bm{x})=f^{2}(\bm{x})$, so are the priors for $f$ and $\bm{u}$, namely $\mathcal{N}(f|\bm{0},\mathbf{K}_{xx^{\prime}})$ and $\mathcal{N}(\bm{u}|\bm{0},\mathbf{K}_{zz^{\prime}})$. More link functions such as $\exp(\cdot)$ are discussed by ? (?). Moreover, we use the same variational joint posterior on $f$ and $\bm{u}$ as Eqn.(10). Consequently, we complete the variational joint distribution on all latent variables as below:

[TABLE]

and notations of VBHP are summarized in Table 1.

3.2 CELBO

Based on Eqn.(8)&(17), we obtain the CELBO for VBHP (see details in A.1 of the appendix (?)):

[TABLE]

where $H_{B}=-\sum_{i=1}^{N}\sum_{j=0}^{i-1}q_{ij}\log q_{ij}$ is the entropy of the variational posterior $B$. The KL terms are between gamma and Gaussian distributions for which closed forms are provided in A.2 of the appendix (?).

3.3 Data Dependent Expectation

Now, we are left with the problem of computing the data dependent expectation (DDE) in Eqn.(19). The DDE is w.r.t. the variational posterior probability $q(B,\mu,f)$. From Eqn.(17), $q(B,\mu,f)=\int q(B,\mu,f,\mathbf{u})\text{ d}\mathbf{u}=q(B)q(\mu)q(f)$ and $q(f)$ is identical to Eqn.(11). As a result, we can compute the DDE w.r.t. $q(B)$ first, and then w.r.t. $q(\mu)$ and $q(f)$.

Expectation w.r.t. $q(B)$. From Eqn.(3), we easily obtain $\log p(\mathcal{D},B|f,\mu)$ by replacing $\phi$ with $f^{2}$, whereupon it is clear that only $b_{ij}$ in $\log p(\mathcal{D},B|f,\mu)$ is dependent on $B$. Therefore, $\operatorname{\mathbb{E}}_{q(B)}[\log p(\mathcal{D},B|f,\mu)]$ is computed as:

[TABLE]

where $f_{ij}\equiv f(\bm{x}_{i}-\bm{x}_{j})$.

Expectation w.r.t. $q(f)$ and $q(\mu)$. We compute the expectation w.r.t. $q(f)$ and $q(\mu)$ by exploiting the expectation and the log expectation of the Gamma distribution: $\operatorname{\mathbb{E}}_{q(\mu)}(\mu)=kc$ and $\operatorname{\mathbb{E}}_{q(\mu)}(\log(\mu))=\psi(k)+\log c$, and also the property $\operatorname{\mathbb{E}}(x^{2})=\operatorname{\mathbb{E}}(x)^{2}+\operatorname{\text{Var}}(x)$:

[TABLE]

where $\psi$ is the Digamma function. We provide closed form expressions for $\int_{\mathcal{T}_{i}}\operatorname{\mathbb{E}}_{q(f)}^{2}(f)$ and $\int_{\mathcal{T}_{i}}\operatorname{\text{Var}}_{q(f)}(f)$ in A.3 of the appendix (?). As in VBPP, $\operatorname{\mathbb{E}}_{q(f)}(\log f_{ij}^{2})=-\widetilde{G}(-\nu_{ij}^{2}/(2\Sigma_{ij}))+\log(\Sigma_{ij}/2)-C$ is available in the closed form with a hyper-geometric function, where $\nu_{ij}=\nu(\bm{x}_{i}-\bm{x}_{j})$, $\Sigma_{ij}=\Sigma(\bm{x}_{i}-\bm{x}_{j},\bm{x}_{i}-\bm{x}_{j})$ ($\nu$ and $\Sigma$ are defined in Eqn.(11)), $C\approx 0.57721566$ is the Euler-Mascheroni constant and $\widetilde{G}$ is defined as $\widetilde{G}(z)=\prescript{}{1}{F}_{1}^{(1,0,0)}(0,1/2,z)$, i.e., the partial derivative of the confluent hyper-geometric function $\prescript{}{1}{F}_{1}$ w.r.t. the first argument. We compute $\widetilde{G}$ using the method of ? (?) and implement $\tilde{G}$ and $\tilde{G}^{\prime}$ by linear interpolation of a lookup table (see a demo in the supplementary material).

3.4 Predictive Distribution of $\phi$

The predictive distribution of $f(\bm{x})$ depends on the posterior $\mathbf{u}$. We assume that the optimal variational distribution of $\mathbf{u}$ approximates the true posterior distribution, namely $q(\mathbf{u}|\mathcal{D},{\bm{\theta}^{\prime}}^{*})=\mathcal{N}(\mathbf{u}|\bm{m}^{*},\mathbf{S}^{*})\approx p(\mathbf{u}|\mathcal{D},\bm{\theta})$. Therefore, there is $q(f|\mathcal{D},{\bm{\theta}^{\prime}}^{*})\approx p(f|\mathcal{D},\bm{\theta})$, i.e., the approximate predictive $f(\tilde{\bm{x}})\sim\mathcal{N}(\bm{K}_{\tilde{x}z}\bm{K}_{zz}^{-1}\bm{m}^{*},\bm{K}_{\tilde{x}\tilde{x}}-\bm{K}_{\tilde{x}z}\bm{K}_{zz^{\prime}}^{-1}\bm{K}_{z\tilde{x}}+\bm{K}_{\tilde{x}z}\bm{K}_{zz^{\prime}}^{-1}\mathbf{S}^{*}\bm{K}_{zz^{\prime}}^{-1}\bm{K}_{z\tilde{x}})\equiv\mathcal{N}(\tilde{\nu},\tilde{\sigma}^{2})$. Given the relation $\phi=f^{2}$, it is straightforward to derive the corresponding $\phi(\tilde{\bm{x}})\sim\text{Gamma}(\tilde{k},\tilde{c})$ where the shape $\tilde{k}=(\tilde{\nu}^{2}+\tilde{\sigma}^{2})^{2}/[2\tilde{\sigma}^{2}(2\tilde{\nu}^{2}+\tilde{\sigma}^{2})]$ and the scale $\tilde{c}=2\tilde{\sigma}^{2}(2\tilde{\nu}^{2}+\tilde{\sigma}^{2})/(\tilde{\nu}^{2}+\tilde{\sigma}^{2})$.

4 New Variational Inference Schema

We now propose a new variational inference (VI) schema which uses a tighter ELBO than the common one, i.e.

Theorem 1**.**

For VBHP, there is a tighter ELBO

[TABLE]

Remark*.*

TELBO is tighter because it is equivalent to the CELBO (Eqn.(19)) except without subtracting non-negative KL divergences over $\mu$ and $\mathbf{u}$. Other graphical models such as the variational Gaussian mixture model (?) have a similar TELBO. Later on, we propose a new VI schema based on the TELBO.

Proof.

With the variational posterior probability of the branching structure $q(B)$ defined in Eqn.(16) and through the Jensen’s inequality, we have:

[TABLE]

where $H_{B}$ is the entropy of $B$ defined in Eqn.(19). The term $\sum_{B}q(B)\log p(\mathcal{D},B)$ can be understood as follows. Consider that infinite branching structures are drawn from $q(B)$ independently, say $\{B_{i}\}_{i=1}^{\infty}$. Given a branching structure $B_{i}$, the Hawkes process can be decomposed into a cluster of Poisson processes, denoted as $(\mathcal{D},B_{i})$, and the corresponding log-likelihood is $\log p(\mathcal{D},B_{i})$. Then, $\sum_{B}q(B)\log p(\mathcal{D},B)$ is the mean of all log likelihoods $\{\log p(\mathcal{D},B_{i})\}_{i=1}^{\infty}$,

[TABLE]

where $n_{B}$ is the number of occurences of branching structure $B$. Since all branching structures $\{B_{i}\}_{i=1}^{\infty}$ are i.i.d., the clusters of Poisson processes generated over $\{B_{i}\}_{i=1}^{\infty}$ should also be independent, i.e., $\{(\mathcal{D},B_{i})\}_{i=1}^{\infty}$ are i.i.d.. It follows that

[TABLE]

We compute the CELBO of $\log p(\{(\mathcal{D},B_{i})\}_{i=1}^{\infty})$ by making $\bm{z}=(\mu,f,\mathbf{u})$ and $\bm{x}=\{(\mathcal{D},B_{i})\}_{i=1}^{n}$ in Eqn.(8):

[TABLE]

Further, we plug Eqn.(28) and Eqn.(30) into Eqn.(27):

[TABLE]

where (a) is because the finite values of KL terms are divided by infinitely large $n$, and (b) is due to i.i.d. $(\mathcal{D},B_{i})$ and the variational posterior $B$ being independent of $f$ and $\mu$. Finally, we plug the above inequality into Eqn.(25) and obtain the TELBO. ∎

New Optimization Schema for VBHP

To optimize the model parameters, we employ the expectation-maximization algorithm. Specifically, in the E step, all $q_{ij}$ are optimized to maximize the CELBO, and in the M step, $\bm{m}$, $\bm{S}$, $k$ and $c$ are updated to increase the CELBO. We don’t use the TELBO to optimize the variational distributions because it doesn’t guarantee minimizing the KL divergence between variational and true posterior distributions. Instead, the TELBO is employed to select GP hyper-parameters:

[TABLE]

The TELBO bounds the marginal likelihood more tightly than CELBO, and is therefore expected to lead to a better predictive performance — an intuition which we empirically validate in Sec.6.

The updating equations for $q_{ij}$ are derived through maximization of Eqn.(19) under the constraints $\sum_{j=0}^{i-1}q_{ij}=1$ for all $i$. This maximization problem is dealt with the Lagrange multiplier method, and yields the below updating equations:

[TABLE]

where $A_{i}=\theta\exp(\psi(k))+\sum_{j=1}^{i-1}\exp(\operatorname{\mathbb{E}}_{q(f)}(\log f_{ij}^{2}))$ is the normalizer.

Furthermore, and similarly to VBPP, we fix the inducing points on a regular grid over $\mathcal{T}$. Despite the observation that more inducing points lead to better fitting accuracy (?; ?), in the case of our more complex VBHP, more inducing points may cause slow convergence (Fig.5a (?)) for some hyper-parameters, and therefore lead to poor performance in limited iterations. Generally, more inducing points improve accuracy at the expense of longer fitting time.

5 Acceleration Trick

Time Complexity Without Acceleration

In the E step of model optimization, updating $q_{ij}$ requires computing the mean and the variance of all $f_{ij}$, which both take $O(M^{3}+M^{2}N^{2})$ with $N$ points in the HP and $M$ inducing points. Here, we omit the dimension of data $R$ since normally $M>R$ for a regular grid of inducing points. Similarly, in the M step, computing the hyper-geometric term requires the means and variances of all the $f_{ij}$. Finally, computation of the integral terms takes $O(M^{3}N)$. Thus, the total time complexity per iteration is $O(M^{3}N+M^{2}N^{2})$.

Acceleration to Linear Time Complexity

To accelerate our VBHP, similarly to ? (?) we exploit the stationarity of the triggering kernel, assuming the kernel has negligible values for sufficiently large inputs. As a result, sufficiently distant pairs of points do not enter into the computations. This trick reduces possible parents of a point from all prior points to a set of neighbors. The number of relevant neighbors is bounded by a constant $C$ and as a result the total time complexity is reduced to $O(CM^{3}N)$.

Specifically, we introduce a compact region $\mathcal{S}=\bigtimes_{r=1}^{R}[\mathcal{S}_{r}^{\min},\mathcal{S}_{r}^{\max}]\subseteq\mathcal{T}$ so that $\phi(\bm{x}_{i}-\bm{x}_{j})=0$ and $q_{ij}=0$ if $\bm{x}_{i}-\bm{x}_{j}\not\in\mathcal{S}$. As a result, all terms related to $\bm{x}_{i}-\bm{x}_{j}\not\in\mathcal{S}$ vanish. To choose a suitable $\mathcal{S}$, we again use the TELBO, taking the smallest $\bm{S}$ for which the TELBO doesn’t drop significantly; we optimize $S_{r}^{\min}$ and $S_{r}^{\max}$ by grid search with other dimensions fixed (so that this step is run $R$ times in total) and we optimize $S_{r}^{\min}$ after optimizing $S_{r}^{\max}$.

Rather than selecting pairs of points in each iteration in the manner of Halpin’s trick (?; ?), our method pre-computes those pairs, leading to gains in computational efficiency. The similar aspect is that both tricks have hyper-parameters to select to threshold the triggering kernel value. We employ the TELBO for hyper-parameter selection while frequentist methods use the cross validation.

6 Experiments

Evaluation. We employ two metrics: the first is the $\mathbf{L_{2}}$ distance (for cases with a known ground truth), which measures the difference between predictive and truth Hawkes kernels, formulated as $\text{L}_{2}(\phi_{\text{pred}},\phi_{\text{true}})=(\int_{\mathcal{T}}(\phi_{\text{pred}}(\bm{x})-\phi_{\text{true}}(\bm{x}))^{2}\text{ d}\bm{x})^{0.5}$ and $\text{L}_{2}(\mu_{\text{pred}},\mu_{\text{true}})=|\mu_{\text{pred}}-\mu_{\text{true}}|$; the second is the hold-out log likelihood (HLL), which describes how well the predictive model fits the test data, formulated as $\log p(\mathcal{D}_{\text{Test}}=\{\bm{x}_{i}\}_{i=1}^{N}|\mu,f)=\sum_{i=1}^{N}\log\lambda(\bm{x}_{i})-\int_{\mathcal{T}}\lambda$. To calculate the HLL for each process, we generate a number of test sequences by every time randomly assigning each point of the original process to either a training or testing sequence with equal probability; HLLs of test sequences are normalized (by dividing test sequence length) and averaged.

Prediction. We use the pointwise mode of the approximate posterior triggering kernel as the prediction because it is computationally intractable to find the posterior mode at multiple point locations (?). Besides, we exploit the mode of the approximate posterior background intensity as the predictive background intensity.

Baselines. We use the following models as baselines. (1) A parametric Hawkes process equipped with the sum of exponential (SumExp) triggering kernel $\phi(\bm{x})=\sum_{i=1}^{K}a_{1}^{i}a_{2}^{i}\exp(-a_{2}^{i}\bm{x})$ and the constant background intensity. (2) The ordinary differential equation (ODE) based non-parametric non-Bayesian Hawkes process (?). The code is publicly available (?). (3) Wiener-Hopf (WH) equation based non-parametric non-Bayesian Hawkes process (?). The code is publicly available (?). (4) The Gibbs sampling based Bayesian non-parametric Hawkes process (Gibbs Hawkes) (?). For fairness, the ARD kernel is used by Gibbs Hawkes and corresponding eigenfunctions are approximated by Nyström method (?), where regular grid points are used as VBHP. Different from batch training in (?), all experiments are conduct on single sequences.

6.1 Synthetic Experiments

Synthetic Data. Our synthetic data are generated from three Hawkes processes over $\mathcal{T}=[0,\pi]$, whose triggering kernels are sin, cos and exp functions respectively, shown as below, and whose background intensities are the same $\mu=10$:

[TABLE]

As a result, for any generated sequence, say $\{x_{i}\}_{i=1}^{N}$, $\mathcal{T}_{i}=[0,\pi-x_{i}]$ is used in the CELBO and the TELBO.

Model Selection. As the marginal likelihood $p(\mathcal{D}|\bm{\theta})$ is a key advantage of our method over non-Bayesian approaches (?; ?), we investigate its efficacy for model selection. Fig.2b shows the contour plot of the approximate log marginal likelihood (the TELBO) of a sequence. It is observed that the contour plot of the TELBO has agreement to the contour plots of $\text{L}_{2}(\phi)$ (Fig.2c) — GP hyper-parameters with relatively high marginal likelihoods have relatively low $\text{L}_{2}$ errors. Fig.2a plots the posterior triggering kernel corresponding to the maximal approximate marginal likelihood. Similar agreement is also observed between the TELBO and the HLL (Fig.3a, 3b). This demonstrates the practical utility of both the marginal likelihood itself and our approximation of it.

Evaluation. To evaluate VBHP on synthetic data, 20 sequences are drawn from each model and 100 pairs of train and test sequences drawn from each sample to compute the HLL. We select GP hyper-parameters of Gibbs Hawkes and of VBHP by maximizing approximate marginal likelihoods. Table 2 shows evaluations for baselines and VBHP (using 10 inducing points for trade-off between accuracy and time, so does Gibbs Hawkes) in both $\text{L}_{2}$ and HLL. VBHP achieves the best performance but is two orders of magnitudes slower than Gibbs Hawkes per iteration (shown as Fig.3c and 3d). The TELBO performs closely to the CELBO in the $\text{L}_{2}$ error and this is also reflected in Fig.2c where the maximum points of the TELBO and the CELBO overlap. In contrast, the TELBO consistently improves the performance of VBHP in the HLL, which is also reflected in Fig.3b where hyper-parameters selected by the TELBO tend to have a higher HLL. Interestingly, when the parametric model SumExp uses the same triggering kernel (a single exponential function) as the ground truth $\phi_{\exp}$, SumExp fits $\phi_{\exp}$ best in $\text{L}_{2}$ distance while due to learning on single sequences, the background intensity has relatively high errors. Although our method is not aware of the parametric family of the ground truth, it performs well. Compared with non-parametric frequentist methods which have strong fitting capacity but suffer from noisy data and have difficulties with hyper-parameter selection, our Bayesian solution overcomes these disadvantages and achieves better performance.

6.2 Real World Experiments

Real World Data. We conclude our experiments with two large scale tweet datasets. ACTIVE (?) is a tweet dataset which was collected in 2014 and contains $\sim$41k (re)tweet temporal point processes with links to Youtube videos. Each sequence contains at least 20 (re)tweets. SEISMIC (?) is a large scale tweet dataset which was collected from October 7 to November 7, 2011, and contains $\sim$166k (re)tweet temporal point processes. Each sequence contains at least 50 (re)tweets.

Evaluation. Similarly to synthetic experiments, we evaluate the fitting performance by averaging HLL of 20 test sequences randomly drawn from each original datum. We scale all original data to $\mathcal{T}=[0,\pi]$ (leading to $\mathcal{T}_{i}=[0,\pi-x_{i}]$ used in the CELBO and the TELBO for a sequence $\{x_{i}\}_{i=1}^{N}$) and employ 10 inducing points to balance time and accuracy. The model selection is performed by maximizing the approximate marginal likelihood. The obtained results are shown in Table 2. Again, we observe similar predictive performance of VBHP: the TELBO performs better the CELBO; VBHP achieves best scores. This demonstrates our Bayesian model and novel VI schema are useful for flexible real life data.

Fitting Time. We further evaluate the fitting speed222The CPU we use is Intel(R) Core(TM) i7-7700 CPU @ 3.60GHz and the language is Python 3.6.5. of VBHP and Gibbs Hawkes on synthetic and real world point processes, which is summarized in Fig.3c and 3d. The fitting time is averaged over iterations and we observe that the increasing trends with the number of inducing points and with the data size are similar between Gibbs Hawkes and VBHP. Although VBHP is significantly slower than Gibbs Hawkes per iteration, VBHP converges faster, in 10$\sim$20 iterations (Fig.5 (?)), giving an average convergence time of 549 seconds for a sequence of 1000 events, compared to 699 seconds for Gibbs Hawkes. The slope of VBHP in Fig.3d is 1.04 (log-scale) and the correlation coefficient is 0.96, so we conclude that the fitting time is linear to the data size.

7 Conclusions

In this paper, we presented a new Bayesian non-parametric Hawkes process whose triggering kernel is modulated by a sparse Gaussian process and background intensity is Gamma distributed. We provided a novel VI schema for such a model: we employed the branching structure so that the common ELBO is maximized by the expectation-maximization algorithm; we contributed a tighter ELBO which performs better in model selection than the common one. To address the difficulty with scaling with the data size, we utilize the stationarity of the triggering kernel to reduce the number of possible parents for each point. Different from prior acceleration methods, ours enjoys higher efficiency. On synthetic data and two large Twitter diffusion datasets, VBHP enjoys linear fitting time with the data size and fast convergence rate, and provides more accurate predictions than those of state-of-the-art approaches. The novel ELBO is also demonstrated to exceed the common one in model selection.

8 Acknowledgements

We would like to thank Dongwoo Kim for helpful discussions. Rui is supported by the Australian National University and Data61 CSIRO PhD scholarships.

Appendix A Omitted Derivations

A.1 Deriving Eqn.(19)

As per Eqn.(8), there is

[TABLE]

where the KL term can be simplified as

[TABLE]

We utilise the likelihood $p(\mathcal{D},B|\mu,f)$ by the reconstruction term and the KL term w.r.t. $B$

[TABLE]

where $H_{B}=-\sum_{B}q(B)\log q(B)$ is the entropy of $B$ and further computed as follows. We adopt $q(B)$ from Eqn.(16) and the close form expression of $H_{B}$ is derived as

[TABLE]

A.2 Closed Forms of KL Terms in the ELBO

[TABLE]

A.3 Extra Closed Form Expressions for Eqn.(23)

[TABLE]

[TABLE]

[TABLE]

where we explicitly express $\mathcal{T}_{i}$ as a Cartesian product $\mathcal{T}_{i}\equiv\bigtimes_{r=1}^{R}[\mathcal{T}_{i,r}^{\min},\mathcal{T}_{i,r}^{\max}]$ for $R$-dimensional data.

Appendix B Additional Experiment Results