Neural Temporal Point Processes: A Review

First, we need to represent each event $(t_{j},m_{j})$ as a feature vector ${\bm{y}}_{j}$ that can then be fed into an encoder neural network (Section 3.2). We consider the features ${\bm{y}}^{\textrm{time}}_{j}$ based on the arrival time $t_{j}$ (or inter-event time $\tau_{j}$) and ${\bm{y}}^{\textrm{mark}}_{j}$ based on the mark $m_{j}$. The vector ${\bm{y}}_{j}$ is obtained by combining ${\bm{y}}^{\textrm{time}}_{j}$ and ${\bm{y}}^{\textrm{mark}}_{j}$, e.g., via concatenation.

The core idea of autoregressive neural TPP models is that event history ${\mathcal{H}}_{t_{i}}$ (a variable-sized set) can be represented as a fixed-dimensional vector ${\bm{h}}_{i}$ Du et al. (2016). We review the two main families of approaches for encoding the past events $({\bm{y}}_{1},\dots,{\bm{y}}_{i-1})$ into a history embedding ${\bm{h}}_{i}$ next.

The main advantage of recurrent models is that they allow us to compute the history embedding ${\bm{h}}_{i}$ for all $N$ events in the sequence in $O(N)$ time. This compares favorably even to classical non-neural TPPs, such as the Hawkes process, where the likelihood computation in general scales as $O(N^{2})$. One downside of recurrent models is their inherently sequential nature. Because of this, such models are usually trained via truncated backpropagation through time, which only provides an approximation to the true gradients Sutskever (2013).

For simplicity, we start by considering the unmarked case and postpone the discussion of marked TPPs until the next section. An autoregressive TPP models the distribution of the next arrival time $t_{i}$ given the history ${\mathcal{H}}_{t_{i}}$. This is equivalent to considering the distribution of the next inter-event time $\tau_{i}$ given ${\mathcal{H}}_{t_{i}}$, which we denote as $P^{*}_{i}(\tau_{i})$.¹¹1We again use $*$ to denote conditioning on the history ${\mathcal{H}}_{t_{i}}$. The distribution $P^{*}_{i}(\tau_{i})$ can be represented by any of the following functions:

1. 1.

   probability density function $f_{i}^{*}(\tau_{i})$

2. 2.

   cumulative distribution function $F_{i}^{*}(\tau_{i})=\int_{0}^{\tau_{i}}f_{i}^{*}(u)du$

3. 3.

   survival function $S_{i}^{*}(\tau_{i})=1-F_{i}^{*}(\tau_{i})$

4. 4.

   hazard function $\phi_{i}^{*}(\tau_{i})=f_{i}^{*}(\tau_{i})/S_{i}^{*}(\tau_{i})$

5. 5.

   cumulative hazard function $\Phi_{i}^{*}(\tau_{i})=\int_{0}^{\tau_{i}}\phi_{i}^{*}(u)du\,.$

In an autoregressive neural TPP, we pick a parametric form for one of the above functions and compute its parameters using the history embedding ${\bm{h}}_{i}$. For example, the conditional PDF $f_{i}^{*}$ might be obtained as

$\displaystyle f^{*}_{i}(\tau_{i})=f(\tau_{i}|{\bm{\theta}}_{i}),$

where

$\displaystyle{\bm{\theta}}_{i}=\sigma({\bm{W}}{\bm{h}}_{i}+{\bm{b}}).$

(2)

Here $f(\cdot|{\bm{\theta}})$ is some parametric density function over $[0,\infty)$ (e.g., exponential density) and ${\bm{W}}$, ${\bm{b}}$ are learnable parameters. A nonlinear function $\sigma(\cdot)$ can be used to enforce necessary constraints on the parameters, such as non-negativity.

The chosen parametric function should satisfy certain constraints. First, the inter-event times $\tau_{i}$ are by definition strictly positive, so it should hold that $\smash{\int_{-\infty}^{0}f_{i}^{*}(u)du=0}$ or, equivalently, $\Phi_{i}^{*}(0)=0$. Second, a reasonable default for many applications is to assume that the TPP is non-terminating, which means that the next event $t_{i}$ will arrive eventually with probability one Rasmussen (2011). This corresponds to the constraint $\int_{0}^{\infty}f_{i}^{*}(u)du=1$ or $\lim_{\tau\to\infty}\Phi_{i}^{*}(\tau)=\infty$. The alternative case—a terminating point process—means that with a probability $\pi>0$ no new events will arrive after $t_{i-1}$. This corresponds to the statement $\int_{0}^{\infty}f_{i}^{*}(u)du=1-\pi$ or, equivalently, $\lim_{\tau\to\infty}\Phi_{i}^{*}(\tau)=-\log\pi$.

Specifying one of the functions (1) – (5) listed above uniquely identifies the conditional distribution $P^{*}_{i}(\tau_{i})$, and thus the other four functions in the list. This, however, does not mean that choosing which function to parametrize is unimportant. In particular, some choices of the hazard function $\phi_{i}^{*}$ cannot be integrated analytically, which becomes a problem when computing the log-likelihood (as we will see in Section 5). In contrast, it is usually trivial to obtain $\phi_{i}^{*}$ from any parametrization of $\Phi_{i}^{*}$, since differentiation is easier than integration Omi et al. (2019). More generally, there are three important aspects that one has to keep in mind when specifying the distribution $P^{*}_{i}(\tau_{i})$:

• •

  Flexibility: Does the given parametrization of $P_{i}^{*}(\tau_{i})$ allow us to approximate any distribution, e.g., a multimodal one?

• •

  Closed-form likelihood: Can we compute either the CDF $F_{i}^{*}$, SF $S_{i}^{*}$ or CHF $\Phi_{i}^{*}$ analytically? These functions are involved in the log-likelihood computation (Section 5), and therefore should be computed in closed form for efficient model training. Approximating these functions with Monte Carlo or numerical quadrature is slower and less accurate.

• •

  Closed-form sampling: Can we draw samples from $P_{i}^{*}(\tau_{i})$ analytically? In the best case, this should be done with inversion sampling Rasmussen (2011), which requires analytically inverting either $F_{i}^{*}$, $S_{i}^{*}$ or $\Phi_{i}^{*}$. Inverting these functions via numerical root-finding is again slower and less accurate. Approaches based on thinning Ogata (1981) are also not ideal, since they do not benefit from parallel hardware like GPUs. Moreover, closed-form inversion sampling enables the reparametrization trick Mohamed et al. (2020), which allows us to train TPPs using sampling-based losses (Section 5.2).

Existing approaches offer different trade-offs between the above criteria. For example, a simple unimodal distribution offers closed-form sampling and likelihood computation, but lacks flexibility Du et al. (2016). One can construct more expressive distributions by parametrizing the cumulative hazard $\Phi_{i}^{*}$ either with a mixture of kernels Okawa et al. (2019); Zhang et al. (2020b); Soen et al. (2021) or a neural network Omi et al. (2019), but this will prevent closed-form sampling. Specifying the PDF $f_{i}^{*}$ with a mixture distribution Shchur et al. (2020a) or $\Phi_{i}^{*}$ with invertible splines Shchur et al. (2020b) allows to define a flexible model where both sampling and likelihood computation can be done analytically. Finally, parametrizations that require approximating $\Phi_{i}^{*}$ via Monte Carlo integration are less efficient and accurate than all of the above-mentioned approaches Omi et al. (2019).

Recently, Soen et al. (2021) have proved that several autoregressive neural TPP models are capable of approximating any other TPP model arbirarily well. However, in practice such flexible parametrizations listed above might be more difficult to train or more prone to overfitting. Therefore, the choice of the parametrization remains an important modeling decision that depends on the application.

Lastly, we would like to point out that the view of a TPP as an autoregressive model naturally connects to the traditional conditional intensity characterization (Equation 1). The conditional intensity $\lambda^{*}(t)$ can be defined by stitching together the hazard functions $\phi_{i}^{*}$

$\displaystyle\lambda^{*}(t)=\begin{cases}\phi_{1}^{*}(t)&\text{if }0\leq t\leq t_{1}\\ \phi_{2}^{*}(t-t_{1})&\text{if }t_{1}<t\leq t_{2}\\ \vdots\\ \phi_{N+1}^{*}(t-t_{N})&\text{if }t_{N}<t\leq T\end{cases}$

(3)

In the TPP literature, the hazard function $\phi_{i}^{*}$ is often called “intensity,” even though the two are, technically, different mathematical objects.

In a marked autoregressive TPP, one has to parametrize the conditional distribution $P^{*}_{i}(\tau_{i},m_{i})$ using the history embedding ${\bm{h}}_{i}$. We first consider categorical marks, as they are most often used in practice.

While conditionally independent models require fewer parameters to specify $P^{*}_{i}(\tau_{i},m_{i})$, recent works suggest that this simplifying assumption may hurt predictive performance Enguehard et al. (2020). There are two ways to model dependencies between $\tau_{i}$ and $m_{i}$ that we consider below.

Negative log-likelihood (NLL) is the default training objective for both neural and classical TPP models. NLL for a single sequence $X$ with categorical marks is computed as

$\displaystyle\begin{split}-\log p(X)=&-\sum_{i=1}^{N}\sum_{k=1}^{K}\mathbb{1}(m_{i}=k)\log\lambda^{*}_{k}(t_{i})\\ &+\sum_{k=1}^{K}\left(\int_{0}^{T}\lambda^{*}_{k}(u)du\right).\end{split}$

(9)

The log-likelihood can be understood using the following two facts. First, the quantity $\lambda^{*}_{k}(t_{i})dt$ corresponds to the probability of observing an event of type $k$ in the infinitesimal interval $[t_{i},t_{i}+dt)$ conditioned on the past events ${\mathcal{H}}_{t_{i}}$. Second, we can compute the probability of not observing any events of type $k$ in the rest of the interval $[0,T]$ as $\exp\left(-\int_{0}^{T}\lambda^{*}_{k}(u)du\right)$. By taking the logarithm, summing these expressions for all events $(t_{i},m_{i})$ and event types $k$, and finally negating, we obtain Equation 9.

Computing the NLL for TPPs can sometimes be challenging due to the presence of the integral in the second line of Equation 9. One possible solution is to approximate this integral using Monte Carlo integration Mei and Eisner (2017) or numerical quadrature Rubanova et al. (2019); Zuo et al. (2020). However, some autoregressive neural TPP models allow us to compute the NLL analytically, which is more accurate and computationally efficient. We demonstrate this using the following example.

Suppose we model $P_{i}^{*}(\tau_{i},m_{i})$ using the “time conditioned on marks” approach from Section 3.4. That is, we specify the distribution $P_{i}^{*}(\tau_{i}|m_{i}=k)$ for each mark $k$ with a cumulative hazard function $\Phi_{ik}^{*}(\tau)$. By combining Equation 5 and Equation 9, we can rewrite the expression for the NLL as

$\displaystyle\begin{split}-\log p(X)=&-\sum_{i=1}^{N}\sum_{k=1}^{K}\mathbb{1}(m_{i}=k)\log\phi_{ik}^{*}(\tau_{i})\\ &+\sum_{i=1}^{N+1}\sum_{k=1}^{K}\Phi_{ik}^{*}(\tau_{i}).\end{split}$

(10)

Assuming that our parametrization allows us to compute $\Phi_{ik}^{*}(\tau)$ analytically, we are now able to compute the NLL in closed form (without numerical integration). Remember that the hazard function $\phi_{ik}^{*}$ can be easily obtained by differentiation as $\phi_{ik}^{*}(\tau)=\frac{\partial}{\partial\tau}\Phi_{ik}^{*}(\tau)$. Finally, note that the NLL can also be expressed in terms of, e.g., the conditional PDFs $f_{ik}^{*}$ or survival functions $S_{ik}^{*}$ (Section 3.3).

Evaluating the NLL in Equation 10 can be still computationally expensive when $K$, the number of marks, is extremely large. Several works propose approximations based on noise-contrastive-estimation that can be used in this regime Guo et al. (2018); Mei et al. (2020).

TPPs can be trained using objective functions other than the NLL. Often, these objectives can be expressed as

$\displaystyle\mathbb{E}_{X\sim p(X)}\left[f(X)\right].$

(11)

Such sampling-based losses have been used by several approaches for learning generative models from the training sequences. These approaches aim to maximize the similarity between the training sequences in ${\mathcal{D}}_{\textrm{train}}$ and sequences $X$ generated by the TPP model $p(X)$ using a scoring function $f(X)$. Examples include procedures based on Wasserstein distance Xiao et al. (2017a), adversarial losses Yan et al. (2018); Wu et al. (2018) and inverse reinforcement learning Li et al. (2018).

Sometimes, the objective function of the form (11) arises naturally based on the application. For instance, in reinforcement learning, a TPP $p(X)$ defines a stochastic policy and $f(X)$ is the reward function Upadhyay et al. (2018). When learning with missing data, the missing events $X$ are sampled from the TPP $p(X)$, and $f(X)$ corresponds to the NLL of the observed events Gupta et al. (2021). Finally, in variational inference, the TPP $p(X)$ defines an approximate posterior and $f(X)$ is the evidence lower bound (ELBO) Shchur et al. (2020b); Chen et al. (2021a).

In practice, the gradients of the loss (Equation 11) w.r.t. the model parameters usually cannot be computed analytically and therefore are estimated with Monte Carlo. Earlier works used the score function estimator Upadhyay et al. (2018), but modern approaches rely on the more accurate pathwise gradient estimator (also known as the “reparametrization trick”) Mohamed et al. (2020). The latter relies on our ability to sample with reparametrization from $P_{i}^{*}(\tau_{i},m_{i})$, which again highlights the importance of the chosen parametrization for the conditional distribution, as described in Section 3.3.

On a related note, sampling-based losses for TPPs (Equation 11) can be non-differentiable, since $N$, the number of events in a sequence, is a discrete random variable. This problem can be solved by deriving a differentiable relaxation to the loss Shchur et al. (2020b).

Prediction is among the key tasks associated with temporal models. In case of a TPP, the goal is usually to predict the times and marks of future events given the history of past events. Such queries can be answered using the conditional distribution $P_{i}^{*}(\tau_{i},m_{i})$ defined by the neural TPP model. Nearly all papers mentioned in previous sections feature numerical experiments on such prediction tasks. Some works combine elements of neural TPP models (Section 3) with other neural network architectures to solve specific real-world prediction tasks related to event data. We give some examples below.

Recommender systems are a recent application area for TPPs. Here, the goal is to predict the next most likely purchase event, in terms of both the time of purchase and the type (i.e., item), given a sequence of customer interactions or purchases in the past. Neural TPPs are especially well-suited for this task, as they can learn embeddings for large item sets using neural networks, similar to other neural recommendation models. Moreover, representing the temporal dimension of purchase behavior enables time-sensitive recommendations, e.g., used to time promotions. For example, Kumar et al. (2019) address the next item prediction problem by embedding the event history to a vector, but are concerned only with predicting the next item type (mark). This approach is extended to a full model of times and events, with a hierarchical RNN model for intra- and inter-session activity represented on different levels, in Vassøy et al. (2019).

Another common application of event sequence prediction is within the human mobility domain. Here, events are spatio-temporal, featuring coordinates in both time and space, potentially along with other marks. Examples include predicting taxi trajectories in a city (e.g., the ECML/PKDD 2015 challenge²²2http://www.geolink.pt/ecmlpkdd2015-challenge/), or user check-ins on location-based social media. Yang et al. (2018) address this task with a full neural TPP, on four different mobility data sets. This extends the approach in DeepMove Feng et al. (2018), where the authors similarly use an RNNs to compute embeddings of timestamped sequences, but limit the predictions to the next location alone.

Other applications include clinical event prediction Enguehard et al. (2020), predicting timestamped sequences of interactions of patients with the health system; human activity prediction for assisted living Shen et al. (2018), and demand forecasting in sparse time series Türkmen et al. (2019a).

In prediction tasks we are interested in the conditional distributions learned by the model. In contrast, in structure discovery tasks the parameters learned by the model are of interest.

For example, in latent network discovery applications we observe event activity generated by $K$ users, each represented by a categorical mark. The goal is to infer an influence matrix ${\bm{A}}\in\mathbb{R}^{K\times K}$ that encodes the dependencies between different marks Linderman and Adams (2014). Here, the entries of ${\bm{A}}$ can be interpreted as edge weights in the network. Historically, this task has been addressed using non-neural models such as the Hawkes process Hawkes (1971). The main advantage of network discovery approaches based on neural TPPs Zhang et al. (2021) is their ability to handle more general interaction types.

Learning Granger causality is another task that is closely related to the network discovery problem. Eichler et al. (2017) and Achab et al. (2017) have shown that the influence matrix ${\bm{A}}$ of a Hawkes process completely captures the notion of Granger causality in multivariate event sequences. Recently, Zhang et al. (2020b) generalized this understanding to neural TPP models, where they used the method of integrated gradients to estimate dependencies between event types, with applications to information diffusion on the web.

Neural TPPs have also been used to model information diffusion and network evolution in social networks Trivedi et al. (2019). The DyRep approach by Trivedi et al. generalizes an earlier non-neural framework, COEVOLVE, by Farajtabar et al. (2017). Similarly, Know-Evolve Trivedi et al. (2017) models dynamically evolving knowledge graphs with a neural TPP. A related method, DeepHawkes, was developed specifically for modeling item popularity in information cascades Cao et al. (2017).

Lack of standardized experimental setups and high-quality benchmark datasets makes a fair comparison of different neural TPP architectures problematic.

Each neural TPP model consists of multiple components, such as the history encoder and the parametrization of the conditional distributions (Sections 3 and 4). New architectures often change all these components at once, which makes it hard to pinpoint the source of empirical gains. Carefully-designed ablation studies are necessary to identify the important design choices and guide the search for better models.

On a related note, the choice of baselines varies greatly across papers. For example, papers proposing autoregressive models (Section 3) rarely compare to continuous-time state models (Section 4), and vice versa. Considering a wider range of baselines is necessary to fairly assess the strengths and weaknesses of different families of approaches.

Finally, it is not clear whether the datasets commonly used in TPP literature actually allow us to find models that will perform better on real-world prediction tasks. In particular, Enguehard et al. point out that two popular datasets (MIMIC-II and StackOverflow) can be “solved” by a simple history-independent baseline. Also, common implicit assumptions, such as treating the training sequences as i.i.d. (Section 5), might not be appropriate for existing datasets.

To conclude, open-sourcing libraries with reference implementations of various baseline methods and collecting large high-quality benchmark datasets are both critical next steps for neural TPP researchers. A recently released library with implementations of some autoregressive models by Enguehard et al. (2020) takes a step in this direction, but, for instance, doesn’t include continuous-time state models.

Many of the metrics commonly used to evaluate TPP models are not well suited for quantifying their predictive performance and have subtle failure modes.

For instance, consider the NLL score (Section 5), one of the most popular metrics for quantifying predictive performance in TPPs. As mentioned in Section 3.4, the NLL consists of a continuous component for the time $t_{i}$ and a discrete component for the mark $m_{i}$. Therefore, reporting a single NLL score obscures information regarding the model’s performance on predicting marks and times separately. For example, the NLL score is affected disproportionately by errors in marks as the number of marks increase. Moreover, the NLL can be “fooled” on datasets where the arrival times $t_{i}$ are measured in discrete units, such as seconds—a flexible model can produce an arbitrarily low NLL by placing narrow density “spikes” on these discrete values Uria et al. (2013).

More importantly, the NLL is mostly irrelevant as a measure of error in real-world applications—it yields very little insight into model performance from a domain expert’s viewpoint. However, other metrics, such as accuracy for the mark prediction, and mean absolute error (MAE) or mean squared error (MSE) for inter-event time prediction are even less suited for evaluating neural TPPs. These metrics measure the quality of a single future event prediction, and MAE/MSE only take a point estimate of $\tau_{i}$ into account. One doesn’t need to model the entire distribution over $\tau_{i}$ (as done by a TPP) to perform well w.r.t. such metrics. If only single-event predictions are of interest, one could instead use a simple baseline that only models $P_{i}^{*}(m_{i})$ or produces a point estimate $\tau_{i}^{\textrm{pred}}$. This baseline can be trained by directly minimizing absolute or squared error for inter-event times or cross-entropy for marks. TPPs are, in contrast, probabilistic models trained with the NLL loss; so comparing them to point-estimate baselines using above metrics is unfair.

The main advantage of neural TPPs compared to simple “point-estimate” methods is their ability to sample entire trajectories of future events. Such probabilistic forecasts capture uncertainty in predictions and are able to answer more complex prediction queries (e.g., “How many events of type $k_{1}$ will happen immediately after an event of type $k_{2}$?”). Probabilistic forecasts are universally preferred to point estimates in the neighboring field of time series modeling Gneiting and Katzfuss (2014); Alexandrov et al. (2020). A variety of metrics for evaluating the quality of such probabilistic forecasts for (regularly-spaced) time series have been proposed Gneiting et al. (2008), but they haven’t been generalized to marked continuous-time event sequences. Developing metrics that are based on entire sampled event sequences can help us unlock the full potential of neural TPPs and allow us to better compare different models. More generally, we should take advantage of the fact that TPP models learn an entire distribution over trajectories and rethink how these models are applied to prediction tasks.

While most recent works have focused on developing new architectures for better prediction and structure discovery, other applications of neural TPPs remain largely understudied.

Applications of classical TPPs go beyond the above two tasks. For example, latent-variable TPP models have been used for event clustering Mavroforakis et al. (2017); Xu and Zha (2017) and change point detection Altieri et al. (2015). Other applications include anomaly detection Li et al. (2017), optimal control Zarezade et al. (2017); Tabibian et al. (2019) and fighting the spread of misinformation online Kim et al. (2018).

Moreover, most papers on neural TPPs consider datasets originating from the web and related domains, such as recommender systems and knowledge graphs. Meanwhile, traditional application areas for TPPs, like neuroscience, seismology and finance, have not received as much attention. Adapting neural TPPs to these traditional domains requires answering exciting research questions. To name a few, spike trains in neuroscience Aljadeff et al. (2016) are characterized by both high numbers of marks (i.e., neurons that are being modeled) and high rates of event occurrence (i.e., firing rates), which requires efficient and scalable models. Applications in seismology usually require interpretable models Bray and Schoenberg (2013), and financial datasets often contain dependencies between various types of assets that are far more complex than self-excitation, commonly encountered in social networks Bacry et al. (2015).

To summarize, considering new tasks and new application domains for neural TPP models is an important and fruitful direction for future work.