Row-clustering of a Point Process-valued Matrix
Authors: Lihao Yin
Ganggang Xu
Huiyan Sang
Yongtao Guan
Venue: NeurIPS 2021 Poster (2021)
Keywords: Expectation-Solution Algorithm
Functional Principal Component Analysis
Marked Point Processes
Model-based Clustering
Semiparametric Model

Abstract: Structured point process data harvested from various platforms poses new challenges to the machine learning community. To cluster repeatedly observed marked point processes, we propose a novel mixture model of multi-level marked point processes for identifying potential heterogeneity in the observed data. Specifically, we study a matrix whose entries are marked log-Gaussian Cox processes and cluster rows of such a matrix.  An efficient semi-parametric Expectation-Solution (ES) algorithm combined with functional principal component analysis (FPCA) of point processes is proposed for model estimation. The effectiveness of the proposed framework is demonstrated through simulation studies and real data analyses.


Row-clustering of a Point Process-valued Matrix
Lihao Yin 1.2
Ganggang Xu
Huiyan Sang²
Yongtao Guan
Institute of Statisties and Big Data, Renmin University, Beijing, China, DDpartment of Statisties, Texas A&M University, College Station, Texas, UUniversity of Miami, Coral Gables, Florida, (lihao, buiyan)@stat tamu. edu, (gangxu, yguan) @bus mi miami... edu
Abstract
Structured point process data harvested from various platforms poses new chal- lenges to the machine learning community. To cluster repeatedly observed marked point processes, we propose novel mixture model of multi-level marked point processes for identifying potential heterogeneity the observed data. Specifically, we study matrix whose entries marked log- Gaussiann Cox processes and cluster rows of such a matrix. An efficient semi-parametric Expectation-Solutio (ES) algorithm combined with functional principal component analysis (FPCA) of point processes is proposed for model estimation. The effectiveness of the proposed framework is demonstrated through simulation studies and real data analyses.
Introduction
Large-scale, high-resolution, and irregularly scattered event time data has attracted enormous research interest recently in many applications, including medical visiting records (13), financial transaction ledgers (28] and server logs (II). Given collection of event time sequences, one research thread is to identify groups displaying similar patterns. In practice, the significance of this task emerges in multifarious scenarios. For example, matching users with similar activity patterns on social media platforms is beneficial to ads recommendations; clustering patients by their visiting records may help predict the course of the disease progression.
Our study is motivated by a dataset we collected from Twitter, which consists of posting times of 500 university official accounts from April 15, May 14th, 2021. Figure|] posting time stamps of seven selected accounts in five consecutive days. While the daily posting patterns vary across different accounts, the posting date seems also play an important role. Specifically, all accounts cascade barrage of postings on April 16th while few postings appear on April 18th. Lastly, each posting associated with specific type of activity namely, tweet, retweet, or reply. Our main interest is to cluster these multi multi-category, dynamic posting patterns into subgroups.
To characterize the highly complex posting patterns, we propose mixture model of Multi-level Marked Point Processes (MM-MPP). We assume that the event sequences from each cluster are realizations of multi-level log-Gaussian Cox process (LGCP) [16), which has been demonstrated useful for modeling repeatedly observed event sequences 28 We here extend their work to the of mixture models and propose semiparametrie Expeetation-SSoutiin algorithm to learn the underlying cluster structure. The proposed learning algorithm avoids iterative numerical optimizations within each ES step and hence computationally efficient. In particular, the expectation step is carried out using MCMC samples based on the FPCA of the latent Gaussian processes, which imposes minimal parametric assumptions on the proposed model. Finally, we design an algorithm that can take advantage of array programming and GPU acceleration to further speed up computation.
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

Evert Types
reply
reweet
UF
uodevis
Stanford
RutgersU
NCState
MIT
04-18
04-17
01-18
04-19
04-20
04-21
Date
Figure 1: The activities of selected accounts on Twitter.
2
Related Work
Modelling of Event Sequences Point processes have been widely used to model temporal events [5], although rarely does existing work focus repeatedly observed event sequences. One prominent example is the Hawkes process B 32 331, which accounts for temporal dependence among events by self-triggering mechanism. However, the existing Hawkes process may fail to describe our cases for two reasons. First, many human activities naturally have discontinuity by day. So it is unclear how define the triggering mechanism across days with Hawkes processes. Second, the multivariate Hawkes process characters an overall rate of events for different days. The clustering methods based on the Hawkes process [15
14 30 are more likely to distinguish individuals by overall event frequency, other than their intra-day behavior patterns.
In our motivating example, there exist multiple variations for the event sequences, both from individual and day levels. One way to account for variations from multiple sources is to exploit Cox process models, whose intensities are modeled by latent random functions. One popular class of Cox processes is the log-Gaussian Cox process (LGCP) [16], whose latent intensity functions are assumed to be transformed Gaussian processes. Recently, 28 proposed multi-level LGCP model to account for different sources of variations for repeatedly observed event data. However, clustering of repeatedly observed marked event time data was considered in their work.
Clustering of Event Sequences. Extensive research has been done on this topic. To our knowledge, clustering models for point processes can be summarized into two major categories: distance-bart based clustering S [8] and distribution- based clustering 330 13 The former measures the closeness between event sequences based clustering algorithms such as k-means [ 19 or EM algorithms [27]. The second approach, also referred to model-based clustering, assumes that event sequences are derived from a parametric mixture model of point processes. One notable thread is the mixture model of the Hawkes point processes. For example, (30] proposed Dirichlet mixture of Hawkes processes (DMHP) under the Exppecation-Maxinile
extracted features and then uses classical distanced- based
(EM) framework to identify clusters. However, existing EM algorithms for event sequence clustering have common issue that they typically require iterative numerical optimizations within each M-step, which would drastically overburden the computation. This computational issue will be accentuated when event data are repeatedly observed and have marks.
3 Model-based Row-clustering for a Matrix of Marked Point Processes
Notation. Suppose that we observe daily event sequences from accounts during m days. For account i on day j. let Nr,, denote the total number of events, tz.j.1 E (0, T denote the I-th event time stamp, and ri.j. E {I, R} denote the corresponding event types (marks). The activities of account on day j can summarized by a set S.j = {(L.,,.t ri..,)) No. recording the time stamps and types for all Ni, events. This general notation can also describe other marked event sequences which are repeatedly observed on m non-overlapping time slots. We represent the collection of all marked daily event sequences as x matrix S, whose (i, jhth entry is marked event sequence

---

St.. We aim to cluster the rows of S to identify potential heterogeneity in account activity patterns, while taking into account the dependence across rows and columns to characterize the complex event patterns and interactions among accounts, days, and event types.
3.1 A Mixture of Multi-level Marked LGCP Model
Given matrix of daily event sequences S, we can separate each matrix entry S according their marks (event types). Let S = {ta.aa/aa.a = r be the collection of time stamps of event type r E {1,.... R}. We model each S by an inhomogeneous Poisson point process conditional on a latent intensity function \1 (t)A = exp(Ar,,(t) where (t) : [o, T] + IR is the random log intensity function on [0, 77. Following 228 we assume a multi-level model for (t):
N (1) (7() + + (T) 2(, (t),),
t [0, T)
for i = 1,, nn 1, m and 1.... R. In model ①. (T) (T()) and Zi (t) random functions on [T, characterizing the variations of account-level. day-level and the residual deviations, respectively. In addition, we also take into account the dependence across event types when modelling& XT (4), Yÿ (t) and (t). while assuming independence across accounts, that is, for any (r, r'), (I) and X (1) are independent when i + i', (1) and Yj (1) are independent when j J and (1)(1) and y (1) are independent if (i.j) + (f)
We assume that (t()) = (x) is a mixture of multivariate Gaussian processes with C' compo- nents in order to detect heterogeneous clusters. We introduce binary vector w, {wi,i, wc, Y to encode the cluster membership for account i, where We,s = if account i belongs to the c-th cluster and 0 otherwise. In analogy to other model-based clustering approaches, the unobserved cluster membership w; are treated as missing data and assumed to follow a categorical distribution (TI,.. TT), where Te indicates the probability that an account belongs with parameter Conditional on ग, we assume that ([)')'s in different clusters have heteroge- to the e-th cluster.
neous behavioral patterns, characterized by their corresponding cluster-specifie multivariate Gaus- sian processes with mean functions _ (t) = E[ (t) = 1] and cross covariance functions rer (s, Cov(XT X (t)lwe, 1], for s, E [0, and r, characterizes the cluster-specific first-order intensity function, and r5: (s, () describes the temporal dependence patterns within and across event types in the same cluster c, C = 1,... ,C
1,..., ,R Here, (2,a(t)
Similarly, we assume that Y (t) [t) and (t) (t) are both mean-zen multivariate Gaussian processes account for dependence of day-level and residual random effects within and across event types, respectively. The covariance functions take the forms: r7" (t)
Cov(Y3 (Yo (1), YJ' (1)). and r (t) Cov Cov(Zr, (1), ZT (t). As the heterogeneity patterns are assumed to be mainly explained by the account-level effect X, both rr.e' (t) and I' (t) are assumed to be homogencous across all clusters.
Single-level Special Case. When m = 1, our data matrix S only has one column of event sequences. The multi-level model in ① reduces to single-level model:
(M,, Mr = XT(),, EE [0,
where X,(1) = (X] has the same model specification in the multi-level case described earlier. We remark that it is still of importance to consider this special case that has also been studied in the literature (29], as even in this simpler case limited work has been done for the clustering of repeatedly observed marked point processes.
3.2 Likelihood Function
We denote the parameters concerning X,(t) in cluster as Oz,c and the parameters concerning
Y, and Zı.5(() (t) as G, and O., respectively. Therefore, the parameters in model I consist of
Ω = TT, ,y O,, Or.., C l.... C). When m 1, Ω = { Or,e; 1.... C} representing
the parameters involved in model (2). The complete data D consists of the observed data S and the
unobserved latent variables {[w; p C). where L { { (t)), (Y,(), (t)} for model (1
3

---

and L = ({){ for model (2) Let S; be the i-th row of S representing activities of the i-th account. our mixture model, the probability of the observed data S can be written as
((S: E. PP(((. (wi)-=11
(3)
where the expectations taken with respect to the conditional distribution of latent variables L and w/ ss and CC( is the conditional probability of a Poisson point process,
So No Aig)ds)
PP(S: (4)
II 1
X,((1A)))) exp [-
(4)
j=1rm1
where, conditional on L, (1) has the form as
for m > 1 and as
(2 for m 1.
Row-clustering Algorithms
Existing mixture model-based clustering methods typically rely on likelihood- Expectation- Maximization algorithms [1] by treating unobserved latent variables, {[w;)i=1: L} in our case, as missing data. However, standard EM algorithms are computationally intractable for the models we consider here. One computation bottleneck is the numerical optimizations involved in M-steps, which require many iterations due the lack of closed-form solutions when updating parameters. Moreover, the computation burden is severely aggravated by the fact that the expectations in E-step (see (3) for an example) involve an intractable multivariate integration.
In Section 4.
we describe a novel efficient semi- caniiparametric Expectation-Soluie algorithm for linelelevel model in (2) to bypass the computation challenges deseribed above. We then show in Section 4.2 that the learning task of multi-level models in (I) can be transformed and solved by utilizing algorithm similar to that of single-level models.
the single- leel
4.1 Learning of Single-tevel Models
The ES algorithm 6 is general iterative approach to solving estimating equations involving missing data or latent variables. The algorithm proceeds by first constructing estimating equations based on a complete-data summary statistic, which may arise from a likelihood, a guasi-likelihood or other generalized estimating equations. Similar to the EM algorithm, the ES algorithm then iterates between an expectation (E)-step and a solution (S)-step until convergence to obtain parameter estimates. The detailed steps of general ES algorithm framework framework a special case of ES when estimating equations are constructed from full likelihoods and using complete data the summary statistic.
included in Supplementary S.2. The EM
Due to the lack of closed-fe form for the likelihood function (3). we opt to design our algorithm under the more flexible and general ES framework for parameter estimations of the single-level models in (2) i.e., m The algorithm is summarized in Algorithm and detailed below.
As preliminary, we give the form of the expectation of the conditional intensity function given cluster memberships as follows:
+ +T,, p. = (t),,(t)) | We,s = 1] = (t (t),
(5)
The form of the second conditional intensity function is
pax (s) (1) We,t 1| [[exp[X7 XI (1)we, 1}}
p= (t) exp[I x.c t)],
(6)
for =--==,, 1,... ".....,,,,rr1.. =1, R, where the last equality is derived following the moment generating function of Gaussian random variable.
carefully construct estimating equations of unknown parameters with three considerations in mind: (1) the expectation of the estimating equations over the complete data should be zero; (2) the conditional expectation of the estimating equation can be solved efficiently in the S-step; (3) the estimating equations should be fast calculate.
Estimating Equations.

---

Let (-)) be a kernel function and (t) = h h (t/hh) with bandwidth h > 0. We define
usto
Ka(s w) v))
A (s, tt;h)) = Σ Word (s,t;h), where of (s,t;h)=
rg((s; (g(s;))g(t;1h)
WEST
Ka - u)
(t; = Σ Wa by ht;
where (1(t;h)
Σ WEST
ng(t; h)
for c L.... C. and r, r,r' I,... R. where g(z;h) Kn(2 - t)dt. Using the Campbell's
Theorem [17] and the moment generating function of the normal distribution, it is straightforward to
show that E A (s, t; k)]lw s (=pp( (t) expli x,e (s, t)] and that E (B) [)]w] (t)
provided that h sufficiently small. This motivates to consider the following estimating equations:
A5 (s, (2 (s) (t) exp (s, =0
(t: h) (tP
(8)
We,i O.
Expectation (E-step). Given an observed data S and current parameter estimate Q we calculate the conditional expectation of the estimation equations in (8) Note that, conditioned on S, the three estimating equations are all linear with respect fwe.a, 1, C, 1,... n). Therefore, the conditional expectations of the estimating equations are obtained by replacing Wey with its conditional expectation EEw (8: which has the following form:
E ( | = xC [(Silwo. = 1;52")
(9)
_ क.. (S.lwa,, = 1;5²)*
where [(S, 1:22) (PP(S [)
Wed 1:52" Here PP(.) is the conditional distribution
function of 5; given We,s and f2" as defined in its Monte Carlo counterpart,
We propose to approximate f(S (we, = =1;0" by
f(S. We, = 1:5²) R0 Q
Σ
PP(S (t)).
(10)
where Q is the Monte Carlo sample size, x (4) (t)'s are independent samples from the multivariate
Gaussian process with parameters O (see details below), and PP(() numerical quadrature approximation PP(() following 2
PP(S, X()) exp Σ uE
Σ
Vu(yaX"(1)) (u) exp (v)).
(11)
In the above, Sh11 is the union of Sh and of regular grid points, v. is the quadrature weight corresponding to each u and yv = w. A.,. where A, an indicator of whether u an observation (A,, = 1) or grid point (A,, 0).
Solution (S-step). In this step. we update the parameters by finding the solutions the expected estimating equations from the E-step. For = 1,... C, rrr = =,, R and r * r', the solutions take the following closed forms:
72
E [W
(12)
FE [
(8, l; S
(s, log 2,0
(13)
Eu BB (s: h S: E. B (t: (t:
(C) log -:=- "E [ (t: [S: /2 exp v (t, t)/2
(14)
5

---

Sampling Strategies. The multi-dimensional functional form of X! renders the sampling proce- dures in (TO) intractable. Given 2,,, one solution is to find low-rank representation of X, with the functional principal components analysis (FPCA) [20]. Specifically,
approximate the latent Gaussian process X, in cluster e nonparametrically using the Karhunen-LLowe expansion [25] as: (t) M. E." @(t), for 1, is a vector of orthogonal eigenfunctions. Using FPCA, we can obtain the samples of X; by sampling c, indirectly. More detailed sampling procedure via FPCA can be seen in Supplementary S.2.
1,... .....
R, where ç, is vector of normal random variables, and (tt)
Remarks. The most significant advantage of our method that it avoids expensive iterations inside each E-step and S-step, unlike other EM algorithms for mixture point process models [30] The elements o (s, h))ss and by (t; h)'s in 7 can be pre-caleulated before E-S iterations to save computations. Moreover, the S-step is fast to execute thanks to the closed form solutions. We will analyze the overall computation complexity of the learning algorithm in Section 4.3
Algorithm Learning of the Single-level model in
Input: S [S,)/= the number of clusters C, the bandwidth h;
Output: Estimates of model parameters, #, (t). fng'' (s, (), for f 1, C,,,,'
,,,,,,,,,r, "r,'',".,,,,.,..................
11... ,,...,R. ,
Calculate the components a²" (s, 1; hi'ss and ht; given in Initialize Q = (T", O% R.c. c= 1,.., C} randomly;
Repeat:
E-Step:
Calculate Ew [ as
Calculate Ew [A (s, (5,t))::27 12" and E. (t) as linear combinations of E [We, S2''''s.
S-Step:
Update π."
(t) -
and Irr' (s, according
and (44
End;
Until:
Reach the convergence criteria;
á n²", Me. (1) and inr'' (s, () ro.' T
Model Selection. Algorithm requires choosing a proper number of clusters C and bandwidth h. In model-based clustering, one popular method for choosing the number of clusters is based on the Bayes information criterion (BIC) (2]. which can be readily computed for our model since the probability f (SSw) is already calculated in each iteration. The choice of bandwidth h also plays an important role in model estimation. A small h may produce unstable clustering results while
large h would dampen the characteristics of each cluster. With a." (s, t; /s'ss and 69 (s; fi'ss given in pre-caleulated for different candidates of h, we can adaptively choose the h that maximizes the likelihood in each iteration in a computationally efficient manner.
4.2 Learning of Multi-level Models
We now consider developing the learning algorithm for the multi-level model []. assuming we repeatedly observe R types of events from accounts on m days with m > Below, we propose method to transform the learning task of a multi-level model into a problem that can be solved by two-step procedure, where the second step is mathematically equivalent to a single-level model and hence can be conveniently solved by a similar algorithm as in Algorithm 1.
For given account consider the aggregated event sequence Sr. 1 for each row of S and event type r. If we assume a multi-level model for each SY as in ①.. conditional on latent variables C, St. is a superposition of m independent Poisson processes and hence can be viewed as a new Poisson process with intensity functional X (1 = Ey exp (L). We approximate the distribution of ST. by a Poisson process with a marginal intensity function,
(t(() = Ey 2 = exx)(((;())
(15)
where X; (X),.... X new multivariate mixture Gaussian process with mean function (1) 16 + TT (1, + (T. l)/2 and covariance function fr-y' 2,0 (s, = re.e' (s, 1),
account i belongs to cluster When m is large, expect the above approximation to be accurate.

---

Note that the model in (1)) for the aggregated event sequence S. is inherently reduced to single- level model. It allows us to separate the inference of the multi-level model in ( into two steps: (Step 1) learn the parameters in 6, and e₂ and denote the estimated parameters as f (s, t) and inr' (s, t); (Step II) learn the clusters of the single-level model in 15| and estimate the parameters n F and Prr''. s.c Afterwards, the parameters involved in 2z can be obtained by
(z, = (t) " y (t.8)/2- 1 tt,,
Fur'' a.e (s, = F (s.
For the learning task in Step I, (28 developed a observed event sequences. estimate G, and G in our mixture multi-level model (I) and provide the details in Supplementary S.1. For step II, we resort to the single-level model algorithm described Section
algorithm to learn the repeatedly analogy to their work, we propose a similar inference framework to
4.3 Computational Complexity and Acceleration
Assume that the training event sequences belong to n accounts and C' clusters and are repeatedly observed on m time slots. We also assume that the data contains types of events and each sequence consists of I time stamps on average. Let Q be the sampling size used in the Monte Carlo integration in 10 In numerical implementation, we divide the interval [T,T] T into D equally spaced grid points D (0- U1 <
up T). In Step I, it requires ( R D2) computation complexity to
estimate Oy and ..,, according to [28]. Computation complexity to pre-calculate a." (s. t; h)'s and
ht; in 7 for all s,t D is of the order (rmmi D²) if we decomposition af"" (s, (s,t; h) as:
of of (s,t;h)
Σ WEST
Kn(s- u) g(s:h)
Kn(t - v) g(t:h)
Σ
Kn u) Kn v) 9(s;h)g(t::n h)
PEST
NESTOS?"
In Step II, for each E-S iteration, need O(CQR*) for sampling and O(nCIQR²) for other calculations. Therefore, the overall computational complexity is O (nm.2 + +QR nCIQ)). To further reduce computation. we use array programming and GPU acceleration calculate the high-dimensional integration in the Monte Carlo EM framework (26) to reduce the runtime of (9) The details are included Supplementary S.2, and numerical demonstration is given in Section
Numerical Examples
We examine the performance of our MM-MPP framework for clustering event sequences via synthetic data examples and real-world applications and compare the performances between the proposed method and two other state-of-the-ar methods. One competing method is a discrete Fréchet distance- based method (DF) by [18]. Unlike other distance-based clustering methods, the DF cluster can characterize interactions among events. Another is model-based clustering method based on the Dirichlet mixture of Hawkes processes (DMHP) by [30]. DMHP is chosen as a competitor due to its capability of accounting for complex point patterns while performing clustering and making efficient variational Bayesian inference algorithms under nested EM framework.
5.1 Synthetic Data
Setting. We generate the synthetic data from the proposed mixture model of log-Gaussian Cox
processes in ( and (2). in which there are R = 5 event types and daily time stamps reside in [0, 2]. We set the number of clusters C from 2 to 5 and set the number of accounts in each cluster to 500. We experiment with an increasing number of replicates (m 1, 20 or 100), to check the convergence of our method. When m 1, we generate event sequences from the single-level model in (2) without day-level variations. In this case, we compare the clustering results of DF, DMHP with those of the single-level model (MS-MPP) When m or 100, we generate data from the multi-level model in (1) and use the MM-MPP method to model the scenario where event sequences are repeatedly observed. However, the two competing methods, DF and DMHP, are not directly applicable for repeated event sequences. Therefore, in this case, we concatenate {ST J sequentially into new event sequence S. on [0, mI] and then apply DF and DMGP to this new sequence. The detailed
7

---

settings of (t)''ss, s. YJ (t)'s and (t)'s is and other details of synthetic data examples are elaborated in Supplementary S.3.
Results. We evaluate the clustering performance of each method over 100 repeated experiments under each setting, using clustering purity (21) as evaluation metric. Table I reports the averaged clustering purity of each method on the synthetic data. When m = 1, MS-MPP obtains the best clustering result in terms of purity consistently across different numbers of clusters. Especially when C' increases, in which case there are more overlaps among clusters, the advantage of MS-MPP becomes more prominent. When m 20, 100, MM-MPP still significantly outperforms the other two competitors. It is also noticeable that the performance of DF and DMPH, in general, deteriorates as m increases, although more repeated event sequences offer more information for clustering. One explanation is that both DF and DMHP may incur bias due to ignoring different sources of variations for repeatedly observed event times. Another reason may be that many existing Hawkes process models, such as DMHP, assume a constant triggering function over time, which may not be flexible enough to characterize the data generated from models and
Table Clustering Purity on Synthetic Data.
20
100
C
DF
DMPH
MS-MPP
DF
DMPH
MM-MPP
DF
DMPH
MM-MPP
2
0.597
0.537
0,831
0.536
0,513
0.947
0.532
0.522
0,988
3
0,514
0.466
0.767
0.465
0.423
0.902
0.477
0.394
0.967
4
0.443
0.421
0.714
0.422
0.356
0.874
0.436
0.285
0,944
5
0.379
0.354
0.675
0.351
0.298
0.835
0.333
0.276
0,919
Our code can be accessed hiahttps:///ithhhh cow com/LihhoVin//MMMP To show the computational advantage of the proposed ES algorithm over the EM algorithm, 2] the computation times of CPU-based EM, CPU-based ES, and GPU-based ES algorithms for 20 iterations in the estimation of model [ with n = 500, 100, m = 20, R = and C = 3. For each iteration, 10, 000 MCMC samples are drawn to approximate (TO). Table 2 demonstrates that with the GPU acceleration, the computation time of the proposed ES can be reduced by more than 20 folds in this case scenario compared to the EM algorithm, which
suitable for array programming
Table 2: Running Time (in seconds) on Synthetic Data
Methods and devices
n = 500
n 1000
GPU-ES (RTX 8000 48G GPU)
30.09
51.42
CPU-ES (i7-7700H CPU) CPU-EM (17-7700HQ CPU)
275.87
505.07
568.36
1105.46
5.2 Real-world Data
In this section, we apply our method to the following real-world datasets.
Twitter Dataset. The Twitter dataset consists of the postings of the official accounts of America's top 500 universities from April 15, 2021, to May 14, 2021. The data set was scraped from Twitter with the API rtweet 02). The dataset involves three categories of postings (tweet, retweet, and reply). indicating At = 3 in this study. As a result, the dataset contains 500 Twitter accounts for m 30 consecutive days with total of 233, 465 time stamps.
Chicago City Taxi Dataset The City of Chicago collected the information of all taxi rides in Chicago since 2013 Each trip record in the dataset consists of drivers' encrypted IDs, pick-up/droppoff time stamps, and locations (in the form of latitude/longitude coordinates). We gathered the trips of 9,000 randomly selected taxi drivers from Jan to Dec 31, 2016, and more than 19 million trip records were picked. We mapped the pick-up coordinates to their corresponding zoning types according to Chicago Zoning Map = which divides the city into nine basic zoning districts including Residence (R), Business (B), Commercial (C), Manufacturing (M), etc. For this data set, we have
n = 9000, m = 366, and R 9.
org/
http:://sconndctyyoonnggg zones/
8

---

Credit Card Transaction Dataset. The dataset contains 641, 914 transaction records of 5, 000 European credit card customers (n 5000) during the period covering January December 31, 2016 (m = 366). We applied the univariate model (R = 1) without event marks to the dataset.
We evaluate and compare clustering stability based
measure called clustering consistency via K- trial cross-validations [23] [24], as there are no ground truth clustering labels. The detailed definition of clustering consistency and other real data example details are included in Supplementary S.4.
Results. We compare the performance of DE, DMHP, and MM-MPP in terms of clustering consistency for three data sets with K 100 trials. The results in Table suggest that MM-MPP outperforms its competitors notably, demonstrating that our model can better characterize the postings patterns and offer more stable and consistent clustering than other methods. Figure 2 shows the histograms of the number of learned clusters for each method. For the Twitter dataset, the median numbers of learned clusters are 3, and 8 for MM-MPP, DMHP, and DF respectively. Besides, the distribution of the number of clusters from our method seems to be the least variable, indicating robustness in clustering. The robustness of our method may be partly attributed to the flexibility of the latent conditional intensity functions that account for multi-level deviations within each account. In contrast, other methods that fail account for different sources of deviations may treat them sources of heterogeneity and consequently result in more clusters.
Table 3: Clustering Consistency on
Datasets.
Method
DF
DMHP
MM-MPP
Twitter
0.096
0.275
0.394
Credit Card
0.102
0.331
0.378
Chicago Taxi
0.045
0.142
0.153
Wetood DMP VWARP
II IF B DUHP MAMP
Counts
Counts
10
Number Clusters
Number of Clusters
Figure 2: Histogram of the number of clusters. Left: Twitter dataset; Right: Credit Card dataset;
More stories can be told by the estimated posting patterns. Given predicted membership of account i by C arg max, En (S: 52), Figure 3 displays the estimated curves of f for tweet events Recall AE., is
(r 1), retweet events (r 2) and reply events interpreted as the baseline of intensity functions. This figure shows three different activity modes for the selected Twitter accounts. The universities cluster marked by curves in Figure in general have lower frequency of posting retweets and replies, especially during the daytime. This group includes the most top university in America, such as MIT, Harvard, and Stanford. In contrast, the accounts in cluster are relatively more active in all three types of postings. We further find that many accounts in this cluster belong to the universities with middle ranks.
3) respectively for
We further applied the proposed MM-MPP to the Chicago Taxi dataset. As suggested by BIC, the 9000 taxi drivers are clustered into groups, whose averaged daily pick-up log intensity functions are illustrated in Figure 4 We can see that the taxi drivers are clustered not only according to their pick-up frequency but also by their working schedules. For example, the black curves on Figure 4a) corresponds to the most dominating group, which occupies 23.2% of the sample. Figure displays the curves of average log intensity (the black line) and log intensity for each driver (gray lines) the selected cluster. Figure 4 show the estimated fex (t) for pick-up in commercial, residence and

---

Groupe
Cla11
CLa22
Chel3
Group:
Cust
Cler2
Diss
Group: Clat1
Clar2
Char3
Time
Time
Time
Figure 3: Curves of f2 (t). Left: tweet events; Mid: retweet events; Right: reply events;
manufacturing districts, respectively. While the pick-up events more likely to occur in commercial districts for this group during the daytime, they also tend to pick up passengers at the residential district in the morning and to appear the manufacturing district in the afternoon. These patterns consistent with the schedules of passengers who commute between homes and workplaces.
More results and discussions chase credit dataset are included in our Supplementary file.
Overall
Commercial
(x)
(
(c)
Log Intensity
0.0
1.0
0.0
1.0
Residence
Manufacturing
(d)
0.0
0.2
0.4
0.6
0.8
1.0
??
0.0
1.0
0.0
1.0
Time
Figure 4: Left: Overall log-intensities for all clusters; Right: Log-intensity one selected cluster;
6 Concluding Remarks
We propose a mixture of multi-level marked point processes to cluster repeatedly observed marked event sequences. A novel and efficient learning algorithm is developed based ES algorithm. The proposed method is demonstrated to significantly outperform other competing methods in simulation experiments and real data analyses.
The current model only focuses on events over temporal domains. However, clustering of spatial patterns on 2- 3-dimensional domains has also attracted much research interest [10 31 9 will
be an interesting research topic to extend the current model to such settings.
This work has foreseeable negative societal impacts, but users should be cautious when giving interpretation on clustering results avoid any misleading conclusions.
Acknowledgement
We thank the anonymous reviewers for their constructive comments that help improve the manuscript significantly. Xu's research supported by NSF grant SES-1902199 Guan's research supported by NSF grant SES-1755575, and Sang's research is supported by NSF grant DMS-185415
10

---

References
[1] Murray Aitkin and Donald B Rubin. Estimation and hypothesis testing in finite mixture models. Journal of the Royal Statistical Society: Series B (Methodological), 47(1):67-75, 1985.
[2] Mark Berman and Rolf Turner Approximating point process likelihoods with glim. Journal of the Royal Statistical Society: Series C (Applied Statistics), 41(1):31-38, 1992.
[3] Donald J Berndt and James Clifford. Using dynamic time warping to find patterns in time series. In Proceedings of the 3rd International Conference on Knowledge Discovery and Data Mining, pages 359-370, 1994.
[4] Paul S Bradley and Usama M Fayyad. Refining initial points for k-means clustering. In ICML, volume 98, pages 91-99, Citeseer, 1998.
[5] Daryl J Daley and David Vere-Jones, An introduction to the theory of point processes: volume I: elementary theory and methods. Springer, 2003.
[6] Michael Elashoff and Louise Ryan. An EM algorithm for estimating equations. Journal of Computational and Graphical Statistics, 13(1):48-65, 2004.
[7] Charles Harris, K Jarrod Millman, Stéfan van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J Smith, et al. Array programming with numpy. Nature, 385(7825):35--36,, 2020.
[8] Alan G Hawkes and David Oakes. A cluster process representation of a self-exciting process. Journal of Applied Probability, pages 493-503, 1974.
[9] Kristian Bjørn Hessellund, Ganggang Xu, Yongtao Guan, and Rasmus Waagepetersen. Semi- parametric multinomial logistic regression for multivariate point pattern data. Journal of the American Statistical Association, pages 1-16, 2021.
[10] Anders Hildeman, David Bolin, Jonas Wallin, and Janine B Illian. Level set Cox processes. Spatial statistics, 28:169 12211699999 193, 2018.
[11] Husna Sarirah Husin, Lishan Cui, Herny Ramadhani Husny Hamid, and Norhaiza Ya Abdullah. Time series analysis of web server logs for an online newspaper. In Proceedings of the 7th International Conference on Ubiquitous Information Management and Communication, pages 1-4, 2013.
[12] Michael W Kearney. rtweet: Collecting and analyzing twitter data. Journal of Open Source Software, 4(42):1829, 2019.
[13] Thomas A Lasko. Efficient inference of gaussian gaussian-proces: prouess-modulated renewal processes with application medical event data. In Uncertainty artificial intelligence: proceedings of the... conference. Conference on Uncertainty in Artificial Intelligence, volume 2014, page 469. NIH Public Access, 2014.
[14] Liangda Li and Hongyuan Zha. Dyadic event attribution social networks with mixtures of hawkes processes. In Proceedings of the 22nd ACM international conference Information & Knowledge Management, pages 1667-1672, 2013.
[15] Dixin Luo, Hongteng Xu, Yi Zhen, Xia Ning, Hongyuan Zha, Xiaokang Yang, and Wenjun Zhang. Multi-task milei-ddimeiionall hawkes processes for modeling event sequences. In Tounty-Fourth International Joint Conference on Artificial Intelligence, 2015.
[16] Jesper Moller, Anne Randi Syversveen, and Rasmus Plenge Waagepetersen. Log gaussian cox processes. Scandinavian jonddinavianoouraa journal of statistics, 25(3):451-42,, 1998.
[17] Jesper Moller and Rasmus Plenge Waagepetersen. Statistical inference and simulation for spatial point processes. CRC Press, 2003.
[18] Tao Pei, Xi Gong, Shih-Lung Shaw, Ting Ma, and Chenghu Zhou. Clustering of temporal event processes, International Journal of Geographical Information Science, 574351484-510, 2013.

---

[19] Jie Peng, Hans-Georg Müller, et al. Distance- clustering of sparsely observed stochastic processes, with applications to online auctions. Annals of Applied Statisties, 1056-10761 2008.
[20] JO Ramsay and BW Silverman. Principal components analysis for functional data. Functional data analysis, pages 147-172, 2005.
[21] Hinrich Schütze, Christopher D Manning, and Prabhakar Raghavan. Introduction to information retrieval, volume 39. Cambridge University Press Cambridge, 2008
[22] Gideon Schwarz al. Estimating the dimension of model. Annals of statistics, 662):461-464, 1978.
[23] Robert Tibshirani and Guenther Walther. Cluster validation by prediction strength. Journal of Computational and Graphical Statistics, 14(3):511-522 2005.
[24] Ulrike von Luxburg. Clustering stability: An overview. Machine Learning, 2(3):23--244, 2009.
[25] Satosi Watanabe. Karhunen-loeve expansion and factor analysis: theoretical remarks and application. In Trans. 4th Prague Conf. Information Theory, Statistic Decision Functions, and Random Processes Prague, pages
1965.
[26] Hong-Zhong Wu, Jun-Jie Zhang, Long-Gard Gang Pang, and Qun Wang. Zmcintegral: A package for multi-dimensional monte carlo integration on multi-gpus. Computer Physics Communications, 248: 148:106960 106962, 2020.
[27] Weichang Wu, Junchi Yan, Xiaokang Yang, and Hongyuan Zha. Discovering temporal patterns for event sequence clustering via policy mixture model. IEEE Transactions on Knowledge and Data Engineering. 2020.
[28] Ganggang Xu, Ming Wang, Jiangze Bian, Hui Huang, Timothy Burch, Sandro Andrade, Jingfei Zhang, and Yongtao Guan. Semi cemi-parametric learning structured temporal point processes. Journal of machine learning research, 2020.
[29] Ganggang Xu, Chong Zhao, Abdollah Jalitian, Rasmus Waagepetersen, Jingfei Zhang, and Yongtao Guan. Nonparametric estimation of the pair correlation function of replicated inhomo- geneous point processes. Electronic Journal of Statistics, 14(2)33730-3765, 2020.
[30] Hongteng Xu and Hongyuan Zha. A dirichlet mixture model of hawkes processes for event sequence clustering, arXiv preprint arXiv:1701.09117, 2017.
[31] Fan Yin, Guanyu Hu, and Weining Shen. Analysis of professional basketball field goal attempts via a bayesian matrix clustering approach. arXiv preprint arXiv:2010.08-955 2020.
[32] Zhang, C Walder, MA Rizoiu, and L Xic. Efficient non non-parametric bayesian hawkes processes. In IJCAI International Joint Conference Artificial Intelligence, 2019.
[33] Feng Zhou, Zhidong Li, Xuhui Fan, Yang Wang, Arcot Sowmya, and Fang Chen. Efficient inference for nonparametric hawkes processes using auxiliary latent variables. Journal of Machine Learning Research, 21(241):1-31, 2020.
12


Reviews:
[{"review_id": "fX7hma8qpyL", "reviewer": "Reviewer_DpLt", "summary": "The problem is equivalent to a marked point process, i.e. a set of times with associated categorical marks. In the main application of Twitter data modelling, the times represent when a tweet was made, and the marks are the cartesian product of twitter handle and tweet type (tweet, retweet, reply).\n\nThe data is split by day and by twitter handle, yielding a 2D matrix of marked point processes with the time (without date) and tweet type for each element of the matrix (so each element of this matrix is a sequence of random length). A heirarchical model is fit to this data, where the intensity for any element of this matrix is a linear combination of that for the day, for the twitter handle, and a third term for the interaction between both, with latent GP for each such term. So far this model comes mainly from [26]. Then, a cluster assignment variable omega is introduced which selects from multiple GP models for the twitter handle component of the intensity function, thereby yielding a mixture model with mixture parameters being GP distributed and capturing the intensity of tweets throughout the day, for different clusters of users.\n\nInference is performed with a form of EM, namely ES or Expectation Solution. This is rather non-trivial and involves significant analytical and numerical adeptness to execute. The paper is therefore technically impressive, and demonstrates excellent technical prowess in probabilistic inference.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is interesting and the model seems interesting for an ML audience, given the rise of social media and associated data.\n\nThe main contribution is in the form of a non-trivial inference scheme. Given the depth of numerical challenges here, it is hard to evaluate the paper. It would have been nice to have received code. At the very least, failing providing the code, it is essential that the experimental setup be reproducible so that subsequent authors can compare. Also, it is a slight disappointment that the paper includes just one significant real-world  numerical evaluation, on a dataset which is somewhat arbitrary in construction and evaluation. However, it remains impressive that the authors could put together this non-trivial scheme and ostensibly make it work.\n\nSome literature is missing. On line 40: \"most existing work rely on strong parametric assumptions\" is questionable. Much recent work uses latent GP intensities for Poisson models and also for latent GP Hawkes' process triggerring kernels. Please cite these works and remove this misleading statement.\n\nThe model is fun but has some strong assumptions, chiefly it is US centric given the segmentation by \"day\". \n\nThe title is a little obfuscated and might make the paper easy to miss.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "MW66ooUkNt6", "reviewer": "Reviewer_J9LW", "summary": "This paper studies clustering based on multi-dimensional data array that can be represented by a point process, which is motivated by Twitter activity data. It proposes a mixture model, where each mixture is a multivariate Gaussian process plus noise and addition structure. In order to identify the mixtures and the components, this paper uses the E-S algorithm, together with several approximation and sampling tricks. Experiments on the Twitter activity data show that the universities under study can be grouped into three clusters, which has meaningful interpretations.", "questions": "", "limitations": "", "rating": 7, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This paper seems to provide a useful statistical methodology for a complex heterogeneous dataset. While I do not claim expertise in this field, this paper does seem to present a practical E-S algorithm for statistical analysis. One drawback of this paper is that the motivating data example looks quite limiting, and that the data example does not reveal interesting discoveries beyond obvious facts (universities can be grouped according to their tiers, etc). I am not criticizing but I feel that a better data example would lead to a much wider appeal to the general audience.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Not applicable.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "8maz7Ln_moD", "reviewer": "Reviewer_qLjo", "summary": "This paper proposes a mixture model of multi-level marked point processes (MM-MPP), which is inspired by a real-world Twitter dataset. A computationally efficient expectation-solution (ES) algorithm is proposed and shown to outperform two existing methods in numerical experiments.", "questions": "", "limitations": "", "rating": 6, "confidence": 1, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This work extends the existing MM-MPP to a mixture model setting, accounting for heterogeneity in data. This is largely inspired by a Twitter dataset, though it is potentially suitable for many other applications as well. \n\nPro: the overall structure of this paper is clean and easy to follow. The models (Sec 3) and algorithms (Sec 4), as well as settings and results of numerical experiments (Sec 5), are presented and explained in detail. In particular, the proposed ES algorithm admits close-form solution in its intermediate steps, which allows efficient computation.\n\nCon: the proposed ES algorithm doesn't come with any theoretical guarantees, and one of my concern is about its convergence. It is well known that typically EM-type algorithm can get stuck at local minimums; the current paper has no discussion on this issue. I would suggest checking (at least) empirically whether the proposed algorithm is indeed robust/insensitive to random initialization, e.g. with 100 random initializations, how many times does it converge to a good solution?\n\nA minor suggestion: in Tables 1 and 2, it would be better to include not just a single number, but also a confidence interval, to account for uncertainty/randomness in data and/or algorithms. It is also helpful (and standard) to use bold font to highlight the best performance under each setting.\n\nTypo: Line 7, \"Exponential-Solution\" -> \"Expectation-Solution\".", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "N/A", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "0gTfCPO1GrG", "reviewer": "Reviewer_3ktE", "summary": "This paper studies an interesting problem, which aims to identify potential heterogeneity in the observed data where structured event data are harvested from various platforms. The authors imposes a matrix structure to repeatedly observed marked point processes, and propose a mixture model of multi-level marked point processes. In particular, the data is structured into a matrix, where each entry contains a sequence of events; these events are modeled by marked log-Gaussian Cox processes; the ultimate goal is to cluster rows of such a matrix. They also propose an efficient semi-parametric Exponential-Solution (ES) algorithm combined with functional principal component analysis of point processes for model estimation.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is clearly written and easy to follow. The main contribution of this paper is introducing the ES algorithm, which is similar to EM algorithm, aiming to facilitate the computation in estimating the latent variables by adopting the kernel method. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "- I don't quite understand how the event data can be divided by days? Would this separation cause discontinuity to the original sequence? And how to evaluate the influence of such data preprocessing to the final results?\n- Since the entire paper is based on a key assumption that event sequences from each cluster are realizations of a multi-level log-Gaussian Cox process, can the authors add more explanation why this assumption applies to the most of the applications?\n- The choices of the bandwidth is significant to the performance of the proposed model. The author mentioned the bandwidth is determined adaptively at each iteration via MLE. It seems to be expensive to carry out in practice. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}]

Author rebuttal:

Meta-review: 
The paper studies structured point process data. The model fitting is achieved using the ES algorithm (an extension of the EM algorithm). The problem considered is interesting and relevant to the machine learning community. The benefit of leveraging the structure for point process models is demonstrated using real data. However, it may be challenging to develop a performance guarantee for such a type of algorithm. So the contribution is maybe in developing a structure-aware model and algorithms for model fitting. The paper addressed the computational efficiency of the algorithm.
Decision: Accept (Poster)