Improving Maximum Likelihood Estimation of Temporal Point Process via Discriminative and Adversarial Learning

1 Introduction and RelatedWork

A major line of research has been devoted to modeling event sequences, especially exploring the continuous timestamp information to learn the underlying dynamics, whereby point process has been a powerful and elegant framework. There is rich literature in point process learning under the maximum likelihood estimation (MLE) framework, is aimed to model the joint distribution of events in the sequence. The learned point processes, with their parameters carrying certain implications, can be either used for relational mining [Zhou et al., 2013], or for event prediction by generating the future events [Du et al., 2016]. This paper is aimed to further improve the prediction capability and stability by introducing new learning techniques to the point processes. The approach is based on the general idea of discriminative learning and adversarial learning, which are relatively ignored in the point process learning literature. We briefly introduce the background.

∗Changsheng Li is correspondence author. The work is supported by NSFC 61602176, 61672231 and NSFC-Zhejiang Joint Fund for the Integration of Industrialization and Informatization U1609220.

Maximum likelihood model While notable progress has been made in generative modelling, in many applications we are still far from producing realistic samples. One of the key open questions is what objective functions one should use to train and evaluate generative models [Theis et al., 2016]. The model likelihood is often considered the most principled training objective and most research in the past decades have been focused on maximum likelihood estimation (MLE) and approximations thereof. Recently we have witnessed several new learning techniques such as those based on adversarial networks [Goodfellow et al., 2014] and kernel moment matching [Dziugaite et al., 2015], and they are not (at least on the surface) related to maximum likelihood. Such a deviation fromMLE is often technically motivated by the fact that the exact likelihood is often intractable. Moreover there exist recent methods which further deliberately differ from MLE to avoid undesired behaviour even when MLE is tractable. One such example is the scheduled sampling strategy [Bengio et al., 2015] which obtains empirically improved performance.

Adversarial learning Adversarial learning techniques e.g. generative adversarial networks (GANs) [Goodfellow et al., 2014; Arjovsky et al., 2017] have been recently applied for data generation in different domains. In theory, GANs can be used to model an arbitrarily complex probability distribution, leading to state-of-the-art results on challenging tasks like image super-resolution [Ledig et al., 2017], and video prediction [Mathieu et al., 2015]. In this paper, we are interested in using adversarial learning for point process based event prediction, which is rarely addressed in existing literature [Xiao et al., 2017a], especially for parametric point process.

Temporal point processes Temporal point processes are fundamental mathematical tool to model event sequences in continuous time space. One key component is its conditional intensity function, defined as the probability of observing an event in an infinitesimal window given the history [Aalen et al., 2008]. Over the decades, different parametric forms of the conditional intensity function are specified, including the Hawkes processes [Hawkes, 1971], self-correcting process [Isham and Westcott, 1979] etc. Based on the given parametric form, existing learning based methods are often devised to optimize the objective of MLE with particular algorithms e.g. sampling technique [Ertekin et al., 2015b], majorization-minimization solver [Lewis and Mohler, 2011] and the ADMM solver [Zhou et al., 2013] etc. Recently there are also emerging works using recurrent neural networks [Du et al., 2016; Xiao et al., 2017b], with the advantage of dismissing the prespecified parametric assumption which is nontrivial and often requires strong domain expert knowledge. As a matter of fact, parametric point processes still have important advantages as prior knowledge can be incorporated (especially in face of insufficient training data), and there is often a clear physical meaning for the parameters [Zhou et al., 2013].

Motivation This paper focuses on enabling gradient descent learning for differentiable parametric point process model. One natural idea is adopting the mean squared error (MSE) as the discriminative loss (though there is technical gap for its use in parametric point process as will be addressed in the paper). Moreover since the ℓ2 loss inherently assumes the data is drawn from a Gaussian distribution, it tends to lead to smoothed prediction curves that fit poorly on multimodal distributions [Mathieu et al., 2015]. As a complementary recipe, GAN technique (and its derived new loss e.g. Wasserstein distance) is adopted to measure the authenticity of predicted events. The hope is to push the predicted curves from temporally evenly distributed sequences to those more realistic-looking ones. However, the above gradient based discriminative and adversarial learning paradigm cannot be directly applied for many existing models. Popular models like the Hawkes process [Hawkes, 1971], Reinforced Poisson Process [Shen et al., 2014] predict the future sequence via sampling techniques e.g. the thinning algorithm [Ogata, 1981], where the loss cannot be computed in a closed form (at least in the general case) for gradient backpropogation.

Novelty Beyond the MLE generative learning paradigm, we propose a novel approach by discriminative and adversarial learning of differentiable point processes. Our approach is fundamentally based on the observation that by discretizing the counting process, the prediction can be approximated by recursively computing the intensity integral for prediction a closed-form function. Note in the recent conditional GAN based work [Xiao et al., 2018], a network based generator is used for event prediction. While in this paper, the generator refers to an explicit parametric point process model leading to a different learning mechanism.

2 Preliminaries

2.1 Intensity Function of Temporal Point Processes

In temporal point processes, each observed point (i.e. event) t (a non-negative real-valued timestamp) is an outcome of the process, forming a sequence S = {ti|i ∈ Z+}. Such sequences are inherently different from time series due to their synchronized nature i.e. the timestamp falls in the continuous domain while time series is formed with equal time interval and the fine-grained time information is lost.

A temporal point process can be characterized by the intensity function, and a cascade (i.e. the observed sequence) can be called a realization of the underlying process. Concretely, denote N(t, t′ ) as the number of points during interval (t, t′ ), Λ(t) = E[N(0, t)] as the expected number of points, the intensity function is defined as:

The key to characterize a point process is to find appropriate (parametric) form of intensity functions. Often, two components are considered. One is the intrinsic component whose value reflects the inherent property of the sequence. For instance, consider a paper citation prediction task, i.e. predicting the future citation events of a paper, the properties may refer to the publication venue, topic, author affiliation, etc. The other component is the external effect generating from the previous events, which can be called predecessorsuccessor dependency, thus the intensity functions become history-dependent. For example, it is commonly believed that the citation exhibits a Mathew effect [Wang et al., 2013] such that a paper with more citations tend to have even more in (short) future. This effect is found ubiquitous for instance for equipment failures [Ertekin et al., 2015a], crime [Mohler et al., 2011], merger and acquisition [Yan et al., 2016] etc.

As a popular embodiment, the intensity function of a multidimensional Hawkes process [Eichler et al., 2017] and the variants [Liu et al., 2017] can be seen as a superposition of background and history effect:

where d refers to event taker, t is timestamp, µ is background, Γd(t) represents the strength that the preceding event affects the successors and g is the affecting kernel function, m demotes the event dimension i.e. event type. It is affected by the preceding events with dimension mj . A specification of these terms is given in the experiment.

2.2 MLE Learning and Sampling based Prediction

The traditional way of training a temporal point process model is usually based on maximum likelihood estimation (MLE) [Ogata, 1988; Lewis and Mohler, 2011]. Concretely, for D observed sequences S = {S1, S2, ..., SD}, and each sequence is represented as Sd = {tdj }Nd j=1 , where tdi ≤ tdj if i < j. Let t0 = 0, tNd = Td as the observation window, the probability density for point t is expressed as [Rubin, 1972]:

Thus the log-likelihood Ld on the whole sequence is:

As each timestamp tj is associated with an event type mj i.e. the so-called dimension, each sequence is denoted as Sd = {(tdj ,mdj )}Nd j=1 . Hence, we have M interdependent intensity functions {λmd (t)}Mm=1 for each type. The resulting overall log-likelihood L can be written as:

Once the parameters are estimated, we can perform in-sample future event prediction via simulation by adopting Ogata's thinning algorithm [Ogata, 1981], which is a counting process based method that generates a sequence point by point. In essence, it is a non-deterministic and sampling process which prevents the use of gradient based analytical learning. Readers are referred to [Ogata, 1981] for more details.

3 Discriminative and Adversarial Learning

Despite the popularity of the MLE models, in many applications one is more interested in predicting future events rather than modeling the joint distribution of past events. This suggests the potential value for a discriminative paradigm to directly boost the prediction performance. Specifically, the discriminative loss used in the paper is the widely used mean square error (MSE) between the predicted event distribution and ground truth in a specified time window (via bining e.g. by year). However, it is nontrivial to adapt existing prediction methods e.g. Ogata's thinning algorithm [Ogata, 1981]. This paper aims to mitigate this gap.

3.1 Enabling Prediction Error Backpropogation

Maximum likelihood estimation (MLE) learning of point process has achieved empirical success in a number of real-world applications, such as social media popularity dynamics prediction [Shen et al., 2014], citation forecasting [Liu et al., 2017] etc. MLE can obtain a learned generative point process model to fit the observed event sequence. As discussed above, however, in prediction tasks the model is used for future event prediction rather than explaining history. We are aimed to improve the prediction capability via discriminative learning.

In fact, adopting discriminative loss for learning is nontrivial to point process models. The main obstacle lies in the prediction (used to compute the loss against ground truth) is non-deterministic in a forward time window because the prediction in the beginning can affect the prediction later in the window e.g. the Hawkes process like models in Eq. 2. Thus existing methods [Ogata, 1981; Møller and Rasmussen, 2006; 2005; Dassios et al., 2013] dominantly use the traditional counting process based simulation algorithms to estimate the predictions over time, which disallows the error signal of predicted sequence from gradient backpropagating.

For prediction, denote ̂ S = {̂ti}Ii=1 as the prediction generated via counting process simulation, whose point number in (t1, t2) is denoted by ̂ N(t1, t2). It is desirable to update point process parameters via error propagation by computing the deviation from the actual count N(t1, t2). However, ̂ N(t1, t2) is in general computed by random sampling from the intensity function, which disallows the analytical gradient computing - more details are given as follows.

We propose an approximation method to enable gradient backpropagation. Recall the intensity function in Eq. 2, as well as many other popular forms as listed in Table 1, the main difficulty is that λmd (t) is nondeterministic over the prediction window (t1, t2). For instance, in the self-exciting Hawkes process, sampling method can generate a few events by some chance during (t1, t1 +∆t), which can in turn increase the intensity during (t1 + ∆t, t2) (refer to the second term on the right of Eq. 2). Hence, the prediction cannot be exactly computed by the integral ∫ t2 t=t1 λmd (t) in a closed form. Because this will ignore the events that may occur during (t1, t2). Similar case also happens for self-correcting point process [Isham and Westcott, 1979] where the recent occurred events decrease the successive event occurrence chance.

Based on the above analysis, our key idea is that for a future time window, one can split it into multiple time units, and the integration is recursively performed on a rolling basis from earlier time units to the later - to get the whole picture, see Eq. 7 where ̂cmd (i) is a function of {̂cmd (j)}j<i. As such, the events occurring in earlier units can be approximately accounted for the prediction in later ones in a more tight fashion.

Formally assume there are D sequences S = {S1, ..., SD}, with each sequence Sd = {(tdj ,mdj )}Nd j=1 where mdj is the event type, we discretize the continuous timestamp tdj by:

where τ is the predefined length of the (short) time interval for discretization. The rounding down operation ⌊·⌋ transforms a sequence Sd with continuous timestamp within [0, Td] into a one indexed with discrete time unit (i.e. bining):

where cmd (i) = N(iτ, iτ + τ) is the number of events of type m occurring during the (i + 1)th time interval unit [iτ , iτ + τ) when the time interval is in the past. For predicted future points, ̂cmd (i) = E[N(iτ, iτ + τ)] is the expected event number which in general is a real number.

This transformation enables to generate future points in a differentiable way w.r.t. the point process parameters. More concretely, without loss of generality given the intensity function in Eq. 2, we can write out the analytical approximation formula for each interval (iτ, iτ + τ):

where P(·) is a distribution determined by the underlying temporal point process, g (and its integration G) is the kernel function, µ (and its integration U) represents intrinsic properties, and Γ represents the predecessor-successor dependency among events. Note the approximate equation in the above formula is a lower-bound (for the Hawkes like self-exciting process) to the true expectation because we only consider the triggering effect from the history citations before iτ ignoring the possible points occurring during (iτ , iτ + τ ). As we assume τ is short such an approximation is reasonable and in practice often becomes natural since the timestamps have limited resolution e.g. year, month in raw data collection.

The above approximation treatment leads to a deterministic and closed-form expectation. We recursively predict future points by Eq. 7 from i-th unit to (i +1)th. In each step, ̂cmd is a real value regarding with the point process parameters (see Eq. 7). The total event count until time t is written as:

3.2 Discriminative and Adversarial Learning

Figure 1 gives the main flow of our approach. We first apply MLE on the whole observed sequence as pre-training to obtain an initial model. This is because we empirically find compared with directly discriminative or adversarial learning, MLE is more insensitive to initial point and can often obtain a good model for further refinement. In fact, there have been welldeveloped MLE learning techniques [Lewis and Mohler, 2011;

Liu et al., 2017] that one can leverage. In contrast, it has been well known that GAN models are practically difficult to train [Arjovsky et al., 2017]. We leave on-the-fly learning of parametric point process model via GAN for future work.

As shown in Fig. 1, we split the whole observed and transformed (see Eq. 6) sequence Sd into two segments (recall that we perform in-sample prediction): i) training sequence Stn d = {cmd (i)}⌊T tn d /τ⌋ i=0 ; ii) validation sequence Svd d = {cmd (i)}⌊T vd d /τ⌋ i=⌊T tn d /τ⌋+1 , i.e. the rest events to validate

the correctness of model's prediction YΘ(Stn d ) conditioned on the input Stn d . By this protocol, we can compute discriminative loss e.g. mean squared error (MSE) between these two event count vectors with size ⌊T vd d /τ⌋ - ⌊T tn d /τ⌋:

where P(Stn d ) is the distribution of training sequence Stn d and Nacc(·) denotes the accumulated total event count in a sequence or period. Given Stn d , the model YΘ(·) recursively predicts {̂cmd (i)}⌊T vd d /τ⌋ i=⌊T tn d /τ⌋+1 based on Eq. 7, which establishes the analytical relation between Θ and the MSE loss.

Due to the limitation that the MSE loss inherently assumes the data is drawn from a Gaussian distribution, and works poorly with multimodal distributions [Mathieu et al., 2015], we apply adversarial training technique to provide the multimodality for generated events. The objective function of our GAN model is written as:

where Pr(·) is the distribution of observed real sequence Sd. The critic FW (·) is a network capable of measuring the Wasserstein distance [Arjovsky et al., 2017] between distribution of generated sequence Pg and that of real one Pr. FW(·) is trained to distinguish from real sequences from model generation. The resulting weighted loss is:

In summary, given a parameterized model YΘ(·) and observed sequences S, first we estimate parameter Θ via MLE on S to obtain a good baseline model. Then we refine YΘ(·) via the proposed adversarial and discriminative learning. The overall learning method is sketched in Algorithm 1.

Note we use the Wasserstein GAN technique [Arjovsky et al., 2017] (specifically the gradient clipping see line 15 in Algorithm 1) instead of the improved WGAN [Gulrajani et al., 2017]. Because we encounter the persistent gradient vanishing issue regarding with the discriminator's penalty terms accounting for the Lipschitz condition as required in improved WGAN. We get more stable convergence by WGAN.

Algorithm 1 Improving MLE based temporal point process modeling via adversarial and discriminative learning.

Require: Observation {Sob d = [Stn d , Svd d ]}Dd=1 (see Fig. 1); its discretization Sob d ; event taker d's profile features xd;

Require: Notations: YΘ(·): generator; FW(·): adversarial critic network; ξ: the MLE learning early stop threshed; Nacc(·): accumulated number of events in a period.

• 1: Initialize point process model parameter Θ randomly;
• 2: for t in 1 to tmle do
• 3: Update Θ by EM based estimation for L(Θ) by Eq. 4 on the observed sequence Sob d ;
• 4: Compute prediction error on validation sequence:

rt = ∑Dd=1 |Nacc

(YΘ(Stn

d )) -Nacc(Svd

• 5: If rt-1 - rt < ξ, break (early stopping);
• 6: end for
• 7: Initialize critic network parameter W randomly;
• 8: while Θ not converge do
• 9: Sample {Stn d }D0 d=1 ∼ P(Stn d ) a batch of D0 sequences

• 11: for t in 1 to tpost do
• 12: Sample {Svd d }D0 d=1 ∼ Pr for validation sequences;
• 13: Sample {Stn d }D0 d=1 ∼ P for training sequences;

• 15: W ← clip(W,-c, c)
• 16: end for
• 17: end while

3.3 Comparison to Network basedMethods

Now we discuss the specific difference of our method to recent work [Xiao et al., 2017a] and its more recent conditional GAN version [Xiao et al., 2018] using network based Wasserstein learning for point process: our model predicts the individual's future events' timestamp and marker (superscript m in Eq. 2) by given its observed preceding sequences. It also means our method allows for personalized model learning for each individual event taker (subscript d in Eq. 2). While in [Xiao et al., 2017a], the model is devised inherently shared for all event takers. Moreover it cannot perform in-sample prediction i.e. for event taker d, predicting its future events based on the history. Instead, the whole sequence can only be generated from scratch in [Xiao et al., 2017a] which hinders the applicability in real-world prediction problems.

More importantly, models in [Xiao et al., 2017a; 2018] are all neural networks making gradient descent trivial. In fact how to adapt to parametric point process calls for more intellectual challenge, which is addressed in this paper. Since parametric point process can have its advantage in terms of better leveraging prior knowledge in face of small data, and higher interpretability compared with black box neural networks, we believe our contribution is orthogonal to [Xiao et al., 2017a; 2018]. In fact, to our best knowledge, we identified no relevant work on adversarial learning of parametric point process.

d )|/Nacc(Svd

d )

4 Embodiment and Experiments

4.1 Prediction Performance EvaluationMetrics

Denote ̂ Cd(t) = ̂ Nd(0, t) as the predicted event number for sequence d before t and Cd(t) as its actual number. Two popular metrics are used to measure the long-term prediction capability in literature [Shen et al., 2014; Liu et al., 2017]:

Mean Absolute Percentage Error (MAPE) It measures the mean relative deviation between predicted and true points count. MPAE is given by (the lower the better):

Accuracy It measures the fraction of sequences correctly predicted for a predefined error tolerance ϵ. The accuracy of popularity prediction is defined by (the higher the better):

Note the above definitions are regarding with a time period [0, t) whereby the total event number is accounted for.

Intensity function modeling We use the intensity function model i.e. Eq. 2 as developed in [Liu et al., 2017] for sequence prediction. Here we briefly introduce the components and readers are referred to the papers for details:

(14)

where model parameter Θ includes α, β, θ and ω, except for xd ∈ RP denoting properties associated with event performer e.g. a paper with its publication venue, institution and its citation events. Note we omit the index p in summation of all profile features for βm, xd, θm d . By Eq. 4 and Eq. 2, the overall MLE log-likelihood function is:

where r(Θ) is an optional sparsity regularizer. An EM algorithm can be used for MLE based parameter estimation. Details can be found in [Liu et al., 2017].

4.2 Dataset and Settings

Synthetic data We simulate D = 1000 synthetic sequences by using the above intensity function over [0, T) for T = 25.

The observation window spans [0, 15], consisting of [0, 10] as training window and [10, 15] as validation. The testing period is [15, 25]. For model's parameters in Eq. 2 and Eq. 14, event take d's profile feature is of size P = 15 i.e. βm, xd ∈ R15 . For each event taker, its profile value xd is sampled by a uniform distribution U(0, 1). The profile encoding parameter βm is sampled from U(0, 3), and θm d is shared among all takers d and is sampled from U(0, 0.5). The triggering parameter αd is set varying over takers and has the form of αd = c ∗

exp (-0.01√d) where c is a sampled from U(0, 6) and d is the d-th sequence from 1 to D. Given the above parameters, the event sequence can be readily simulated using off-the-shelf algorithm e.g. the thinning method [Ogata, 1981].

To study the robustness of our model, we further add noise to the timestamp of each event. Not for learning, we use the noise-add data, while for performance evaluation of the test period, we compare with the ground truth without noise.

Real-world dataset Microsoft Academic Graph (MAG) [Sinha et al., 2015] consists scientific publications with their individual profiles e.g. authors, institutions, venues etc. and citation records over time. We have collected 10,000 publications satisfying: i) published in computer science venues during 1980-1990; ii) each paper has at least 5 citations. For each of these papers, we collect a span of 25 years of citation records since its publishing year. Following the protocol in Fig. 1, the 25-year is split into the observed period (from year 1 to 15) and the test period (from year 16 to 25). Moreover, the observed period is divided into the training period (from year 1 to 10) and the validation period (from year 11 to 15).

We also use NYSE transaction [Du et al., 2016] containing 0.7 million high-frequency transaction records at a stock one day, forming 3200 sequences. We pick the top 1000 sequences with most events for our experiments.

Protocols We use convolutional neural network (CNN) [LeCun et al., 1998] as the critic for classifying genuine sequence Svd d from the generated prediction YΘ(Stn d ). The CNN has a conv-pooling-conv-FC structure with 128 filters. For WGAN [Arjovsky et al., 2017], we set parameter clip c = 0.1, learning stepsize η = 0.01, mixed loss weight δ = 0.5; For detecting overfitting and early stopping, we set ξ = 0.001. RMSProp [Tieleman and Hinton, 2012] is used for backpropagation. For evaluation metrics, we use MAPE (Eq. 12) and ACC (Eq. 13) with different tolerance ϵ in line with [Shen et al., 2014; Liu et al., 2017]. We test MSE (mean square error) as discriminative loss. Three settings are evaluated: i) MLE alone, ii) MLEMSE, iii) MLEMSE+WGAN where the subscripts denote the posterior learning after MLE.

4.3 Results and Discussion

Synthetic data results From the results as shown in Table 2 one can observe our technique can improve MLE.

Real-world data results Table 3 reports both the MAPE and ACC over the first 5 years and on the whole 10 years

(including the first 5-year) respectively. The results show the posterior learning, especially with the mixed loss involving both discriminative and adversarial learning can often improve the MLE baseline. Fig. 2 discloses the yearly performance over the 10-year testing period. One can see combing MSE with WGAN outpeforms MSE only, which suggests of adopting adversarial loss in addition with MSE. Similar results can be found in Table 3 for the NYSE transaction dataset. Moreover, we study the behavior of the weight δ between MSE and WGAN loss in Eq. 11 as shown in Table 4.

5 Conclusion

This paper presents a novel technique for improving maximum likelihood based estimation of temporal point processes. Its utility is verified by testing multi-dimensional Hawkes-like point process on synthetic and real-world sequences. Our model is parametric in contrast to those full RNNs models.

For synergizing the MSE and WGAN losses, in this paper we only tried a naive weighted linear sum with moderate improvement to each other. Seeing the practical complexity of training GAN, better tricks e.g. gradually adding GAN loss may be developed which we leave for future work.