Multi-Task Multi-Dimensional Hawkes Processes for Modeling Event Sequences

1 Introduction

In practical information systems (e.g., social networks, financial systems, IPTV systems), a temporal event sequence of a certain user can be modeled as a point process [Rodriguez et al., 2011; Yang and Zha, 2013; Linderman and Adams, 2014; Li and Zha, 2014]. Typically, an event in the sequence may trigger a series of later events, forming a triggering pattern of the event sequence. On the other hand, multiple event sequences may exhibit similar triggering patterns and therefore can be characterized as a cluster. For example, in an IPTV system, users buy and watch various TV programs, and their watching behaviors (when and what they watch) form a large number of temporal event sequences. For each user, watching an early episode of a drama may trigger the events of watching its following episodes (i.e., self-triggering) and other related news (i.e., mutual-triggering). The triggering pattern among watching behaviors reflects the watching preferences of the user. Furthermore, the sequences of the users having similar preferences can be clustered according to the similarity of their triggering patterns. Simultaneously learning the triggering pattern of each individual event sequence and the clustering structure across all the sequences has great practical significance, because it enables both local and global depictions of the entire dynamic system. However, such a learning task is very challenging as it demands concurrent modeling point processes individually and globally.

∗These two authors contribute equally.

Mathematically, suppose that each user's behaviors can be represented as a event sequence {(t1,E1), ..., (tN,EN)} (0 < ti ≤ T). The event Ei occurs at time ti, which is an element of a event set E = {1, ...,C}. Given a set of event sequences from different users, we aim to 1) find the triggering pattern among events for each user, including the influence of the event c on the event c′ , c, c′ ∈ E, denoted as Au = [acc′ ] ∈ RC×C + and the temporal dynamic of the influence; 2) explore the clustering structure of Au 's.

In this paper, we propose a novel point process model, namely Multi-task Multi-dimensional Hawkes Process (MMHP), to learn triggering patterns and clustering structures from a number of event sequences. Specifically, given a set of temporal event sequences, MMHP models the corresponding dynamic system using an intrinsic intensity matrix, a structured infectivity tensor and a triggering kernel. The intrinsic intensity matrix captures the basic instantaneous happening rate of various events; the tensor represents the infectivity among events, which reveals the triggering pattern among events; the triggering kernel measures the time decay effect of the infectivity. As the learning of event sequences from the same cluster should share implicit relatedness, we learn our model under the multi-task learning framework [Evgeniou and Pontil, 2007; Liu et al., 2009]. We impose sparse and low rank constraints on the infectivity tensor so as to induce the clustering structure among sequences, and propose an effective algorithm that takes advantage of alternating direction method of multipliers (ADMM), majorization minimization and EulerLagrange equations to solve the optimization problem.

Different from traditional multi-dimensional Hawkes process models, which only consider the triggering pattern within sequences, MMHP considers both the triggering pattern within sequences and the clustering structure across sequences. In MMHP, the triggering patterns of sequences and their similarities are captured by the structured infectivity tensor. An advantage of MMHP is that they can avoid overfitting for the sequences where only a few events are observed. For those cases, the structure of the infectivity tensor facilitates the knowledge transfer from sequences with many observed events to those short of observations. We evaluate MMHP and compare it with other state-of-the-art methods. Our experimental results on both synthetic and real datasets demonstrate the superior performance and robustness of MMHP.

2 RelatedWork

Point processes. Point processes [Daley and Vere-Jones, 2007] are popular models for sequential and temporal data, which have been successfully applied to model the occurrences of earthquakes [Ogata, 1988; 1999], the transactions in the stock market [Chavez-Demoulin and McGill, 2012; Bacry et al., 2013], asset management [Yan et al., 2013], meme tracking [Yang and Zha, 2013] and user interactions [Zhou et al., 2013a] on social networks. A point process is typically represented as an event sequence {(t1,E1), ..., (tN,EN)} (0 < ti ≤ T), where the event Ei occurs at time ti. Denote N(t) as the number of points (i.e., events) happening in the time interval (-∞, t] and Ht = {Ei|ti < t} as the set of events happening before t, a point process is characterized by its conditional intensity function

where E(dN(t)|Ht) is the expectation of the number of events happening in the interval (t, t+dt] given historical observations Ht. The conditional intensity function represents the expected instantaneous rate of events at time t.

Hawkes processes. A Hawkes process [Hawkes, 1971; Hawkes and Oakes, 1974] is a point process having a selftriggering property, that is, the occurrences of previous events trigger the occurrences of future events. Its conditional intensity function is defined as follows,

where µ ∈ R+ is an intrinsic intensity for the occurrences of events, g : R+ → R+ is a triggering kernel function quantifying the triggering effects from previous events.

Multi-dimensional Hawkes processes. In most situations, there can be multiple types of events happening sequentially, e.g., multiple users post information and interact with each other on a social network, multiple items are sold on an online shop. In those cases, each event not only triggers the events of its type, but also triggers the events of other different types. Multi-dimensional Hawkes processes are used to model such processes. Specifically, given C types of events, the conditional intensity λ(t) = [λ1(t), ..., λC(t)]⊤ is a sizeC vector, where λc(t) is the conditional intensity for type-c events defined as follows,

where µc is the intrinsic intensity of the event of type c. Compared with (1), an infectivity matrix A ∈ RC×C + is introduced to measure the influence across events of different types, that is, acc′ represents the infectivity of type-c′ events to type-c events. ∑j:tj<t accj g(t - tj) represents the influence of historical events on the instantaneous rate of event at time t.

Multi-dimensional Hawkes processes have been proposed and applied to analyze the topic diffusion [Rodriguez et al., 2011; Du et al., 2013; Yang and Zha, 2013] and the user interactions [Zhou et al., 2013a; Blundell et al., 2012] on social networks, the transactions in the stock market [ChavezDemoulin and McGill, 2012; Bacry et al., 2013], etc. In these works, the event types can correspond to topics, users, transaction types and any other objects. Among these works, Blundell et al. proposed a multi-dimensional Hawkes processes with the Infinite Relational Model (IRM) [Blundell et al., 2012] to simulate and predict the social interactions among users on social networks. Zhou et al. proposed a multi-dimensional Hawkes process model, which learns an infectivity matrix of users explicitly with sparse and lowrank constraints [Zhou et al., 2013a]. Furthermore, the triggering kernel of the model is learned by nonparametric estimation in [Zhou et al., 2013b]. Yang et al. proposed to use multi-dimensional Hawkes process to perform diffusion network inference and meme tracking jointly [Yang and Zha, 2013]. Recently, an online learning algorithm for multi-dimensional Hawkes processes is proposed in [Hall and Willett, 2014], which approximates continuous Hawkes processes in a discrete manner. Multi-dimensional Hawkes processes have achieved promising results in many challenging tasks. However, most of the existing works focus on learning triggering patterns of sequences while few of them consider the clustering structure across sequences. The recent works in [Du et al., 2013; Li and Zha, 2014; Linderman and Adams, 2014] start to explore the relationship among sequences by learning parametric models. Different from these works, we take advantage of the multi-task learning strategy [Evgeniou and Pontil, 2007; Liu et al., 2009; Jacob et al., 2009] and add structural regularization directly to the proposed model.

3 Multi-taskMulti-dimensional Hawkes Processes

We propose the following Multi-task Multi-dimensional Hawkes Process (MMHP) models to learn triggering patterns of each event sequence and structures across various sequences jointly, where each event sequence includes multiple types of events. Given U sequences with C event types, we represent the uth event sequence as Su = {(tui , cui )}nu i=1 , where u = 1, 2, · · · ,U. cui ∈ {1, ...,C} and tui ∈ (0, Tu] are the event type and time stamp of the i-th event in Su, respectively. Tu is the time span of Su, and nu is the number of events in Su. Each sequence Su can be modeled by a multidimensional Hawkes process. Based on (2), the conditional intensity function for Su on type-c event at time t is as follows,

Different from the conventional multi-dimensional Hawkes processes as in (2), the model parameters in (3) fall in two categories: 1) a global triggering kernel g(t), which reflects the attenuating influential effects from historical events and is shared by all the sequences; 2) a local infectivity matrix Au = [aucc′ ] and a local natural intensity vector µu = [µu1 , ..., µuC ]⊤, which are specific to sequence Su. We represent such Au 's and µu 's into an infectivity tensor A = [aucc′ ] ∈ RC×C×U + and an intrinsic intensity matrix µ = [µuc ] ∈ RC×U + , respectively.

3.1 Structural Constraints onMMHPs

We impose structural constraints on the model parameters in (3) so as to enable pattern learning across sequences. Specifically, we impose sparse and low-rank constraints on the flattened matrix A = [vec(A1 ), ..., vec(AU )] ∈ RC2×U from tensor A. The sparsity constraint is based on the observations that typically within each sequence, only a subset of event types happens and triggers others. The intuition behind the low-rank constraint is to uncover clustering structures when the overall triggering patterns of multiple sequences in terms of their infectivity effects are similar and thus introduce similar Au 's. Here we assume that the similarity is described by the self-representation property of infectivity matrices, which implies the low-rank structure of A. These two structural constraints are reasonable for practical systems. For example, in an IPTV system, each user typically has a preference over a small number of program categories, and such a preference can be common among a large number of users.

3.2 Learning Algorithms for MMHPs

We learn our MMHP model via Maximum Likelihood Estimation (MLE). Specifically, we learn the parameters µ, A and g(t) by solving the following optimization problem.

where L(µ,A, g) is the negative log-likelihood of MMHP that can be written as:

The second termR(g) in (4), defined as

regularizes the triggering kernel [Zhou et al., 2013b] to ensure the triggering kernel is smooth and differential energy limited (i.e., ∫ ∞ 0 [g′ (t)]2dt < ∞). The nuclear norm and ℓ1 norm on A in (4) impose low rank and sparsity on A.

To solve the optimization problem (4), we apply the scheme of ADMM [Boyd et al., 2011; Ouyang et al., 2013; Zhou et al., 2013a] and introduce two auxiliary variables Z1

and Z2, two dual variables U1 and U2, and thus convert the problem in (4) as follows.

We solve the problem in (6) by an iterative algorithm, which updates A and µ, g(t), Z1 and Z2, U1 and U2 iteratively.

Step 1: Update A and µ. We first update A and µ by a majorization-minimization algorithm. Given the parameters of k-th iteration as Θ(k) , we find a surrogate function of objective function by Jensen's inequality as follows.

where

Considering the terms related to {µ,A} and setting ∂Q ∂µuc = 0 and ∂Q ∂aucc′ = 0, we obtain close-form solutions for µ = [µuc ] and A = [aucc′ ]:

where

zu(k) cc′ and uu(k) cc′ are the elements in Zu(k) and U u(k) w.r.t. aucc′ , respectively.

Step 2: Update g(t). The triggering kernel g can be updated in an infinite dimensional space [Zhou et al., 2013b]. Considering the terms of (7) related to g(t), the solution of g(t) satisfies the Euler-Lagrange equation:

where g′′ (t) is the second-order derivative of g(t) and

I(·) is the indicator function which returns 1 if the input predicate is true and 0 otherwise. We solve (10) numerically using the following efficient Seidel-type iterations. Specifically, setting the sampling interval as ∆t, we discretize the differential equation over small intervals m∆t, for m = 1, ...,M, as follows:

where gm = g(m∆t), Fm = F (m∆t) and Dm = D(m∆t). M is the number of samples of g and M∆t is the length of g. Therefore, we can solve for gm by fixing all other gm′ , m′ ̸= m and solving the above quadratic equation.

Step 3: Update Z1 and Z2. Given the updated A(k+1) from Step 1, we solve for Z1 by solving the following optimization problem.

The solution is obtained by shrinking the singular values of A(k+1) + U (k) 1 by soft-thresholding [Zhou et al., 2013a] as

where Sλ1/ρ(·) is the soft-thresholding function with the threshold λ1/ρ. Similarly, we solve for Z2 by solving the following optimization problem,

and the solution is

where Eλ2/ρ(·) is the soft-thresholding function of matrix's elements. The threshold is λ2/ρ.

Step 4: Update U1 andU2. The dual variables U1 and U2 are solved as follows based on ADMM.

The whole learning algorithm is summarized in Algorithm 1.

Algorithm 1 MMHP Learning Algorithm

Input: event sequences {Su}, parameters {λ1, λ2,M, α} Output: A, µ and g(t)

repeat

k = k +1

Update {µ(k) ,A(k)} by (8) and (9)

repeat

for m = 1 : M do

m by (11)

Update g(k)

end for

until convergence

Update U (k)

1 , U (k)

until convergence

4 Experimental Results

We evaluate the performance of MMHP on both synthetic and real-world data. Specifically, we compare the learning algorithm MMHP in Algorithm 1 with the following alternatives:

• Full: the infectivity tensor A has no structures (i.e., λ1 = λ2 = 0 in (4)). This method is similar to that in [Zhou et al., 2013b], where they learn triggering patterns for each sequence independently.
• Sparse: only the sparsity constraint is imposed on A (i.e., λ1 = 0 in (4)).
• LowRank: only the low-rank constraint is imposed on A (i.e., λ2 = 0 in (4)).

We use the following metrics to evaluate the performance of various methods:

• LogLik: the log-likelihood of testing data using the trained model.
• EstErr: the averaged estimation error of instantaneous infectivity aucc′g(t), defined as

where { ̂ A, ˆ g(t)} represents real parameters and {A, g(t)} represents the corresponding estimates.

• RankCorr: the averaged Kendall's rank correlation coefficient between each row of the real ̂ A and that of the estimated A.
• ClusAcc: the clustering accuracy, defined as the percentage of sequences clustered correctly based on the learned infectivity tensor. It is only used for synthetic data.
• ClusDiff : a metric of clustering accuracy for real-world data, where the ground truth of parameters and the clustering indices are unavailable. In specific, we first cluster sequences by applying k-means [Zha et al., 2001; Ng et al., 2002] on the estimated A. After computing centers of clusters from estimated A, we then construct a clustered infectivity tensor ˜ A, where the u-th slide ˜ Au is the center of the

2 by (14) and (15)

cluster the Au belongs to. Denote the set of possible methods as I (Here I = {Full, Sparse, LowRank, MMHP}). For a certain method, ClusDiff measures the difference of the log-likelihood calculated from the ever best performing method and the log-likelihood calculated using the constructed ˜ A from that method, that is, ClusDiff of the i-th method is defined as

If a clustering result is good, then each cluster center should be representative for capturing the dynamics of the sequences belonging to the cluster, and the difference of log-likelihood caused by replacing the specific infectivity matrices with the cluster centers should be small.

4.1 Experimental Results on Synthetic Data

We generate a synthetic dataset in which there are U = 40 sequences and C = 5 event types. The sequences are generated so as to fall into two clusters of equal size. We generate the flatten version of the infectivity tensor as A = [thres(u1v⊤ 1 ), thres(u2v⊤ 2 )], where u1, u2 ∈ RC2 + and v1, v2 ∈ RU/2 + are four random vectors with values uniformly distributed over [0, 1]. The function thres(X) randomly sets half of the rows in X as zero. In this way, the generated A is low-rank and sparse, and inherently represents two clusters. We generate an intensity matrix µ from a uniform distribution over [0, 0.001]. We use an exponential kernel g(t) = exp(-t) with t ∈ (0, 20] as the triggering kernel for the synthetic data. Given the above parameters, we simulate 100 training sequences and 100 testing sequences respectively. Each training sequence contains 2500 events, and each testing sequence contains 500 events.

Parameter Studies. We first conduct a parameter study on the four parameters {λ1, λ2, M, α} involved in our MMHP algorithm. The parameters λ1 and λ2 control the rank and sparsity level of the infectivity tensor, and M and α control the sampling rate and the smoothness of the triggering kernel. We first identify the optimal parameter configuration by doing a grid search in the parameter space of λ1 ∈ [0, 0.8], λ2 ∈ [0, 0.008], M ∈ [50, 1500] and α ∈ [1, 10000] on a training set where each sequence has 500 events. With the identified optimal configuration as λ1 = 0.2, λ2 = 0.004, M = 500, α = 100, we fix 3 parameters as their optimal values each time and alter the fourth parameter to train a different MMHP model. Fig. 1 represents the LogLik of such MMHPs as the different parameters vary. As the value of λ1 or λ2 grows larger, the LogLik of MMHPs first increases and then decreases, and this demonstrates the effectiveness of λ1 and λ2 in controlling the rank and sparsity of the infectivity tensor, respectively, and thus the model quality. Similar trends apply for M corresponding to the fact that extremely small M leads to a coarse estimation of the triggering kernel while extremely large M leads to over-fitting. The performance is relatively stable for α ∈ [1, 1000] as for larger α values, the strong smoothness regularization leads to an over-smoothed triggering kernel. Overall, the relatively small performance changes around the optimal parameter configuration demonstrate the robustness of our algorithm with respect to its parameters.

Performance Comparisons. We compare MMHP with the Full, Sparse and LowRank methods on LogLik, EstErr and RankCorr, respectively. For each method, we evaluate its performance when the number of events in each training sequence varies. We run the experiment 10 times with various parameter configurations. For each configuration, except the λ1, λ2 in Full, the λ1 in Sparse and the λ2 in LowRank, which are fixed to be 0's, the rest parameters are sampled from a small neighborhood of the optimal configuration shown in the parameter studies. Fig. 2 presents the averaged results of the 10 runs, which shows that MMHP consistently achieves significantly better performance (i.e., higher LogLik, lower EstErr and higher RankCorr) than other methods over different training sets. This is particularly true when the training sequences have fewer events. The experimental results demonstrate that by learning from multiple sequences concurrently, MMHP is able to leverage information from other sequences for the sequences with fewer events, and thus effectively prevent over-fitting.

Clustering Effects Comparisons. We cluster the training sequences by applying k-means clustering algorithm on the learned infectivity tensors and evaluate the clusters using ClusAcc. The results in Table. 1 demonstrate that MMHP outperforms others in uncovering structures across sequences and the performance difference is more significant when training sequences have fewer events, indicating the effectiveness of MMHP in preventing over-fitting.

4.2 Experimental Results on Real IPTV Systems

We apply MMHP to model the data from a real IPTV system. The dataset is collected from China Telecom, in Shanghai,

Bold numbers correspond to the best performance.

China [Luo et al., 2014], which consists logs of TV program watching events frommultiple users, time stamps for the beginning and endings of each watching session, and the names and the categories (labeled manually) of the TV programs. The dataset contains 2967 users (i.e., U = 2967) and 9,000 TV programs belonging to 25 categories (i.e., C = 25) that these users ever watched during 11 months in 2012.

Modeling Watching Behaviors. We model the watching behaviors of users using MMHP. Specifically, for each user u, the event sequence corresponding to her watching behaviors is Su = {(tui , cui )}nu i=1 , where the event (tui , cui ) represents that the user u watches a program of cui -th class at time tui . Given the sequences of all the users during a period, we learn models using MMHP, Full, Sparse and LowRank, respectively. For all the methods, we set the length of the triggering kernel as 11520 minutes (8 days) and the sampling interval ∆t as 20 minutes(M = 576). Such a configuration ensures that 1) for daily and weekly TV programs, the triggering kernel will capture their periodic influence on its own, and 2) for most TV programs, which are 20-40 minutes in length, the triggering kernel has a good resolution to capture the influence from a previous watching event.

For all the methods, we learn the models from the first N months and test them on the data from the (N + 1)-th month. Running the experiments in the same way as for the synthetic data in Section. 4.1, we obtain the averaged LogLik with respect to N = 2, ..., 10 shown in Fig. 3(a). Similar to the results on the synthetic data, MMHP produces better LogLik than other methods. Additionally, we visualize the triggering kernel learned by MMHP in Fig. 3(b), which clearly shows that the triggering kernel captures the temporal influence decay of a program on its following programs, which corresponds to the expected nature of user watching behaviors well. There are 8 spikes in the kernel with 1 spike per day periodically, which corresponds to the self-triggering pattern of daily programs. The first spike is the highest and corresponds to the influence of previous watching events on that day. It indicates that the mutual-triggering patterns among various program categories mainly exist in the watching behaviors happening in the same day. With time elapsing, the intensity of the spikes is reduced gradually, which corresponds to the decay of influences over time. However, the spike of the 7th day is a little higher than those of its adjacent days, indicating the existence of the self-triggering pattern of weekly programs.

Learning User Clusters. We evaluate the performance of various methods on clustering users using ClusDiff . We cluster the users into 5 clusters by considering the infectivity matrix {Au}Uu=1 as the feature of a certain user and applying k-means. Table. 2 shows that compared with the other methods, MMHP not only achieves the best learning result (maxLogLik(A)) given the data from 11 months, but also obtains the smallest ClusDiff . It means that the clustering results based on the proposed MMHP learning method are reasonable, where the clustering centers are representative for most of users' watching behaviors.

Bold numbers correspond to the best performance.

5 Conclusion

In this paper, we propose a multi-task multi-dimensional Hawkes process model and the corresponding learning algorithm. By considering the sparse and low-rank structure of the infectivity tensor, the proposed model captures the triggering patterns within sequences and the clustering structure across sequences jointly. The proposed learning method has superior performance compared to the other methods on both synthetic data and real-world data. In the future, we plan to propose practical applications based on the model, e.g., personalization of IPTV service.

Acknowledgement. This work is supported in part by NIH R01 GM108341,NSFC-61221001 and NSFC-61420106008.

Appendix

Feasibility of Surrogate Function. The surrogate function Q(Θ|Θ(k) ) in Eq. (7) is induced as follows. Based on Jensen's inequality, we have

The equality holds if and only if µ = µ(k) and A = A(k) . Similarly, we also have

The equality holds if and only if g(t) = g(k) (t) and A = A(k) . Denote the objective function as L(Θ). As a result, the surrogate function satisfies

Therefore,

Euler-Lagrange Equations. Similar to that in [Zhou et al., 2013b], the optimization of Q(Θ|Θ(k) ) w.r.t. to g(t) is equivalent to minimize ∫ ∞ 0 f(g, g′ )dt, where

By Euler-Lagrange equation, the solution satisfies

which corresponds to Eq. (10).