In the derivation of the standard Dirichlet-Multinomial posterior predictive,
we consider a single draw from the Multinomial distribution (i.e., a
categorical distribution) with a Dirichlet prior on the parameter $\vec{p}$.
Usually, this prior is linearly dependent on previous draws from the
distribution. We propose to modify this assumption by using a Dirichlet prior
that depends non-linearly on the history of draws as:

$Dir_{r}(\vec{p}|\vec{\alpha},\vec{N})=\frac{1}{B(\vec{\alpha}+\vec{N}^{r})}\prod_{k}p_{k}^{\alpha_{k}+N_{k}^{r}-1}$
(50)

In Eq.50, the vector $\vec{N}^{r}$ shifts the parameter $\vec{\alpha}$
according to the count of draws allocated to each cluster $k$ up to the nth
draw. The parameter $r\in\mathbb{R}$ controls the intensity of this shift for
each entry of $\vec{N}$.

We demonstrate that the Powered Dirichlet distribution is a conjugate prior of
the Multinomial distribution, by writing Eq.50 as:

$\begin{split}Dir_{r}(&\vec{p}|\vec{\alpha},\vec{N})=\frac{1}{B(\vec{\alpha}+\vec{N}^{r})}\prod_{k}p_{k}^{\alpha_{k}-1}\prod_{k}p_{k}^{N_{k}^{r}}\\\
&\stackrel{{\scriptstyle\text{Eqs.\ref{eq-
Mult}}}}{{=}}\frac{B(\vec{\alpha})\prod_{k}N_{k}^{r}!}{B(\vec{\alpha}+\vec{N}^{r})(\sum_{k}N_{k}^{r})!}Dir(\vec{p}|\vec{\alpha})Mult(\vec{N}^{r}|\sum_{k}N_{k}^{r},\vec{p})\\\
&\propto
Dir(\vec{p}|\vec{\alpha})Mult(\vec{N}^{r}|\sum_{k}N_{k}^{r},\vec{p})\end{split}$
(51)

where the prior on vector $\vec{N}$ is a regular Multinomial distribution of
parameter $N=\sum_{k}N_{k}^{r}$. Note that for certain values of $r$, the
vector $\vec{N}^{r}$ might not be made of integer values; the resulting
Multinomial prior on $\vec{N}^{r}$ must then be expressed in terms of $\Gamma$
functions (see Eq.48) to be valid for $\vec{N}^{r}\in\mathbb{R}^{|\vec{N}|}$.
Distributions of non-integer counts are not new in the literature [70, 84, 39]
and are essentially allowed by the generalized definition of the factorial
function in terms of the gamma function. When $r=1$, we recover the standard
Dirichlet-Multinomial prior on $\vec{p}$ for the $n^{th}$ draw; the history of
draws $\vec{N}$ can be expressed as the result of $N$ independent draws of
equal probability $\vec{p}$. When $r\neq 1$, the prior on $\vec{N}$ is sampled
from a Multinomial distribution in which the number of samples drawn depends
on $r$ as $\sum_{k}N_{k}^{r}$. For instance, let $\vec{N}=(1,2)$ and $r=2$:
the resulting powered conditional Dirichlet prior would then be sampled from a
Multinomial distribution $Mult(\vec{N}=(1,4)|N=5,\vec{p}=(p,1-p))$.

#####  .17.4.c Posterior predictive

We now derive the posterior distribution for the $n^{th}$ draw to belong to a
cluster $c$ given all previous draws. We assume that $\vec{C_{-}}$ represents
all previous realizations up to $n-1$, that is, the cluster to which each
previous draw has been associated. For simplicity of notation, we define the
population of a cluster $k$ at time $n-1$ as
$N_{k}=|\\{C_{i}|i=k\\}_{i=1,2,...,n-1}|$. We are now looking at the
probability distribution of its $n^{th}$ draw to belong to $c$. It is
expressed as the probability of a draw from the categorical distribution given
all previous observations (because there is only one new draw, it is the same
as a Multinomial distribution with parameter $N=1$) combined with the powered
Dirichlet prior defined Eq.50. Then:

$\begin{split}DirCat_{r}(C_{n}=c|\vec{\alpha},\vec{C_{-}})&=\int_{\vec{p}}Cat(C_{n}=c|\vec{p})\underbrace{Dir_{r}(\vec{p}|\vec{\alpha},\vec{N})}_{\textbf{Eq.\ref{eq-
DirMultPrior}}}\\\
=&\int_{\vec{p}}\frac{1}{B(\vec{\alpha}+\vec{N}^{r})}\prod_{k}p_{k}^{c_{k}+\alpha_{k}+N_{k}^{r}-1}\\\
&=\frac{B(\vec{c}+\vec{\alpha}+\vec{N}^{r})}{B(\vec{\alpha}+\vec{N}^{r})}\end{split}$
(52)

where $\vec{c}$ is a vector of the same length as $\vec{\alpha}$ whose
$c^{th}$ entry equals to 1, and 0 anywhere else. Alternative demonstrations of
this result are possible [130, 111].

#####  .17.4.d Powered Chinese Restaurant process

We finally derive an expression for the Powered Chinese Restaurant process
from Eq.52. We recall that $N_{k}=|\\{C_{-i}|i=k\\}_{i=1,2,...,n-1}|$. Taking
back the conditional probability for the $n^{th}$ observation to belong to
cluster $c$ (Eq.52), we have:

$\begin{split}p(C_{n}=&c|\vec{C_{-}},\vec{\alpha})=DirCat_{r}(C_{n}=c|\vec{C_{-}},\vec{\alpha})\\\
=&B(\vec{c}+\vec{N^{r}}+\vec{\alpha})/B(\vec{N^{r}}+\vec{\alpha})\\\
=&\Gamma(N_{c}^{r}+\alpha_{c}+1)\frac{\prod_{k\neq
c}\Gamma(N_{k}^{r}+\alpha_{k})}{\Gamma(1+\sum_{k}N_{k}^{r}+\alpha_{k})}\frac{\Gamma(\sum_{k}N_{k}^{r}+\alpha_{k})}{\prod_{k}\Gamma(N_{k}^{r}+\alpha_{k})}\\\
=&\frac{(N_{c}^{r}+\alpha_{c})}{\sum_{k}N_{k}^{r}+\alpha_{k}}\frac{\prod_{k}\Gamma(N_{k}^{r}+\alpha_{k})}{\Gamma(\sum_{k}N_{k}^{r}+\alpha_{k})}\frac{\Gamma(\sum_{k}N_{k}^{r}+\alpha_{k})}{\prod_{k}\Gamma(N_{k}^{r}+\alpha_{k})}\\\
=&\frac{N_{c}^{r}+\alpha_{c}}{\sum_{k}N_{k}^{r}+\alpha_{k}}\end{split}$ (53)

Every cluster with $N_{c}=0$ (empty clusters) has an identical probability of
getting chosen. Besides, the result is identical if either of them gets
chosen. Therefore, we can express the probability of choosing any empty
clusters as a function of $\alpha=\sum_{k}\alpha_{k}$.

Finally, taking the limit $K\rightarrow\infty$ and defining the limit value
$\lim_{K\rightarrow\infty}\sum_{k}^{K}\alpha_{k}=\alpha$, we find the Powered
Chinese Restaurant Process ($PCRP$):

$PCRP(C_{i}=c|\alpha,C_{1},C_{2},...,C_{i-1})=\begin{cases}\frac{N_{c}^{r}}{\alpha+\sum_{k}^{K}N_{k}^{r}}\text{
if c = 1, 2, ..., K}\\\ \frac{\alpha}{\alpha+\sum_{k}^{K}N_{k}^{r}}\text{ if c
= K+1}\end{cases}$ (54)

The formal derivation of the Powered Chinese Restaurant process in Eq.54 and
the demonstration of its link to the conditional Dirichlet prior on $\vec{p}$
are the first main contribution of this section. Besides, this demonstration
uncovers the link between the prior in Eq.50 and an exotic formulation of the
Multinomial distribution, which has never been considered before. Most
importantly, it highlights that it originates from sampling every new
observation from independent weighted Multinomial distributions of different
parameters $N_{r}=\sum_{k}x_{k}^{r}$. As stated in the introduction, special
cases of the process have already been used in some occasions [104, 65, 124]
but never demonstrated. Furthermore, this formulation generalizes the Uniform
process when $r\rightarrow 0$ [124], the Dirichlet process when $r\rightarrow
1$ and the Pitman-Yor process when
$r_{k}(N_{k})=(log(1-\beta/N_{k})+log(N_{k}))/log(N_{k})$ and
$\alpha(K)=\alpha+\beta K$ (see Eq.47, we recall that $e^{\log x}=x$). The
present expression explicitly allows for controlling the importance of the
“rich-get-richer” property as well as recovering state-of-the-art processes.

Figure 25: Effect of $r$ on the Powered Chinese Restaurant process prior
probability. The grey bars represent the population in each cluster at a given
time. The solid lines represent the prior probability of the next observation
to belong to each of these. For $r=1$, we recover the Dirichlet Process prior.

We illustrate the change in prior probability for an existing cluster to get
chosen induced by the Powered Chinese Restaurant process in Fig.25 – we do not
plot the prior probability for a new cluster to be created. This figure plots
the population of clusters (grey bars) and their associated prior probability
of getting chosen. When $r>1$, the most populated clusters are associated with
a more significant prior probability than in the standard CRP, whereas the
less populated ones have even fewer chances to get chosen; rich-get-more-
richer, the prior on population is peakier on large clusters. On the other
hand, when $r<1$, most populated clusters have fewer chances to get chosen
than in CRP, whereas less populated ones have an increased chance of getting
chosen; rich-get-less-richer, the prior on population is flatter across
clusters of different sizes. In the limit case $r=0$, the clusters’ population
does not play any role anymore; rich get-no-richer, the prior is flat over all
clusters. Note that if we wanted to represent the Pitman-Yor process prior in
this figure, it would correspond to the plot for $r=1$ vertically shifted of
$-\beta$ (such as defined Eq.47) leading to an increased probability of
creating a new cluster of $\beta K$ (not represented in the plot) [93].
Varying the parameter $\alpha$ of $\Delta\alpha$ plays a similar role as
$\beta$ in this situation. It would uniformly shift the prior probability for
each existing cluster to get chosen by $\frac{\Delta\alpha}{K}$ and increase
the probability of creating a new one by $\Delta\alpha$. For Both Pitman-Yor
and Dirichlet processes, the linear dependence of each cluster’s population
does not change.

In Fig. 25, we understand that the Powered Chinese Restaurant process allows
for defining priors from clusters population that are not possible when tuning
the Chinese Restaurant or Pitman-Yor processes. Introducing non-linearity in
the dependence on previous observations allows giving more or less importance
to the “rich-get-richer” property.

####  .17.5 Properties of the Powered Chinese Restaurant process

We will now investigate some key properties of the Powered Chinese Restaurant
process. We recall that $N_{k}$ is the population of the cluster $k$, and
$N=\sum_{k}N_{k}$.

#####  .17.5.a Convergence

######  Proposition 2.

For $N\rightarrow\infty$, the Powered Chinese Restaurant process converges
towards a stationary distribution. When $r<1$, it converges towards a uniform
distribution over all the possible clusters, and when $r>1$, it converges
towards a Dirac distribution on a single cluster.

###### Proof.

We consider a simple situation where only 2 clusters are involved. The
generalization to the case where $K$ clusters are involved is straightforward.
When the clusters’ population is large enough, we make the following Taylor
approximation:

$\begin{split}(N_{i}+1)^{r}&=N_{i}^{r}(1+\frac{1}{N_{i}})^{r}=N_{i}^{r}+rN_{i}^{r-1}+\mathcal{O}(N^{r-2})\end{split}$
(55)

Since the population of a cluster $N_{i}$ is a non-decreasing function of $N$,
we assume that first order Taylor approximation holds when
$N\rightarrow\infty$. Given clusters population at the $N^{th}$ observation,
we perform a stability analysis of the gap between probabilities $\Delta
p(N)=p_{1}(N)-p_{2}(N)$. We recall that the probability for cluster $i$ to get
chosen is $p_{i}(N)=N_{i}^{r}/(\sum_{k}N_{k}^{r})$ and that either of the
clusters is chosen with this probability at the next step (at step $N+1$,
$\Delta p(N+1)=p_{1}(N+1)-p_{2}(N)$ with probability $p_{1}(N)$ and $\Delta
p(N+1)=p_{1}(N)-p_{2}(N+1)$ with probability $p_{2}(N)$). Explicitly the
variation of the gap between probabilities when $N$ grows is written as:

$\begin{split}&\frac{p_{1}(N)(p_{1}(N+1)-p_{2}(N))+p_{2}(N)(p_{1}(N)-p_{2}(N+1))-\Delta
p(N)}{\Delta p(N)}\\\ &\stackrel{{\scriptstyle\text{Eq.\ref{eq-
taylorapprox}}}}{{\approx}}\frac{1}{p_{1}(N)-p_{2}(N)}\times\left(p_{1}(N)\frac{N_{1}^{r}-N_{2}^{r}+rN_{1}^{r-1}}{N_{1}^{r}+N_{2}^{r}+rN_{1}^{r-1}}+p_{2}(N)\frac{N_{1}^{r}-N_{2}^{r}-rN_{2}^{r-1}}{N_{1}^{r}+N_{2}^{r}+rN_{2}^{r-1}}\right)\\\
&\ \
=\frac{2rN_{1}^{r}N_{2}^{r}}{(N_{1}^{r}+N_{2}^{r}+rN_{1}^{r-1})(N_{1}^{r}+N_{2}^{r}+rN_{2}^{r-1})}\left(\frac{N_{1}^{r-1}-N_{2}^{r-1}}{N_{1}^{r}-N_{2}^{r}}\right)\end{split}$
(56)

We see in Eq.56 that the sign of the variation of the gap between
probabilities depend only on the term
$\frac{N_{1}^{r-1}-N_{2}^{r-1}}{N_{1}^{r}-N_{2}^{r}}$. We can therefore
perform a stability analysis of the Powered Chinese Restaurant process using
only this expression.

When $0<r<1$, the sign becomes negative because the following relation holds:
${N_{1}^{r-1}-N_{2}^{r-1}<0\Leftrightarrow N_{1}^{r}-N_{2}^{r}>0\ \forall
N_{1},N_{2}}$; that makes right hand side of Eq.56 negative. Therefore, adding
a new observation statistically reduces the gap between the probabilities of
the two clusters. We could forecast this prediction from Eq.55 by seeing that
adding a new observation to a large cluster increases its probability to get
chosen lesser than for a small cluster– rich-get-less-richer. Moreover, we see
from Eq.55 that a crowded cluster (such as $N_{1}^{r}\gg N_{2}^{r}$) see its
probability evolve as $N^{r-1}$. Asymptotically, the only fixed point of Eq.56
when $N\rightarrow\infty$ is $N_{1}\rightarrow N_{2}$, which implies a uniform
distribution.

On the contrary, when $r>1$ we have the following relation:
$N_{1}^{r-1}-N_{2}^{r-1}>0\Leftrightarrow N_{1}^{r}-N_{2}^{r}>0\ \forall
N_{1},N_{2}$; ; that makes right hand side of Eq.56 positive. Adding a new
observation statistically increases the gap between probabilities. From Eq.55,
we see that adding an observation to a large cluster increases its probability
with its population – rich-get-more-richer. In this case, Eq.56 has $K+1$
fixed points, with $K$ the number of clusters. The uniform distribution is an
unstable fixed point, while $K$ Dirac distributions (each on one cluster) are
stable fixed points of the system. It means the gap converges to $1$, that is
a probability of 1 for one cluster and a probability of 0 for the others.

When $r=1$, the right-hand side of Eq.56 is null. It means the gap remains
statistically constant $\forall N_{i}$, which is a classical result for the
regular Dirichlet process. This convergence has already been studied on many
occasions [35, 6].

We note that as $r\rightarrow 0$, Eq.56 is not defined anymore. That is
because the probability for a cluster to be chosen does not depend on its
population anymore. In this case, $p_{1}(N)-p_{2}(N)\propto
N_{1}^{0}-N_{2}^{0}=0$: the probability for any cluster to be chosen is equal,
hence the Uniform process – “rich-get-no-richer”.

∎

#####  .17.5.b Expected number of tables

######  Proposition 3.

When $N$ is large, $\sum_{k}N_{k}^{r}$ varies with $N$ as
$N^{\frac{r^{2}+1}{2}}$ when $r<1$, and with $N^{r}$ when $r\geq 1$.

###### Proof.

Taking back Eq.54, we are interested in the variation of
$p_{i}=\frac{N_{i}^{r}}{\sum_{k}N_{k}^{r}}$ according to $N$ when $N_{i}^{r}$
is large:

$\begin{split}p_{i}(N+1)-p_{i}(N)\approx\begin{cases}\frac{rN_{i}^{r-1}+\mathcal{O}(N^{r-2})}{\sum_{k}N_{k}^{r}}&\text{
if $N_{i}$ grows}\\\ 0&\text{ otherwise}\end{cases}\end{split}$ (57)

We see in Eq.57 that for $r<1$, the larger $N_{i}$ the slower the variation of
$p_{i}$. It means that for large $N_{i}^{r}$, we can write $N_{i}\propto
Np_{i}$, with $p_{i}$ a constant of $N$. Since $N_{i}$ is either way a non-
decreasing function of $N$, we reformulate the constraint $N_{i}^{r}$ large in
$N^{r}$ large.

For $r>1$, the probability $p_{i}$ varies greatly with $N$ and quickly
converges to 1 for large $N$ (see Proposition 2), and so $N_{i}\approx N$ for
cluster $i$ and $N_{j\neq i}\ll N_{i}\ \forall j$.

Since the sum $\sum_{k}N_{k}^{r}$ essentially varies according to large
$N_{k}$, we can approximate $\sum_{k}N_{k}^{r}\approx N^{r}\sum_{k}p_{k}^{r}$
for large $N^{r}$.

Besides, we showed in Proposition 2 that for large $N$ the process converges
towards a uniform distribution for $r<1$ and towards a Dirac distribution when
$r>1$. Therefore, we can express $\sum_{k}^{K}p_{k}^{r}$ as:

$\sum_{k}^{K}p_{k}^{r}\stackrel{{\scriptstyle N\gg
1}}{{\approx}}\begin{cases}K^{1-r}&\text{for $r<1$}\\\ 1&\text{for $r\geq
1$}\end{cases}$ (58)

Based on the demonstration of Eq.4 in [124], we suppose that $K$ evolves with
$N$ as $N^{\frac{1-r}{2}}$ when $r<1$. We verify that this assumption holds in
the Experiment section.

Therefore, we can write:

$\sum_{k}N_{k}^{r}\approx
N^{r}\sum_{k}^{K}p_{k}^{r}\approx\begin{cases}N^{r}\left(N^{\frac{1-r}{2}}\right)^{1-r}=N^{\frac{1+r^{2}}{2}}&\text{for
$r<1$}\\\ N^{r}&\text{for $r\geq 1$}\\\ \end{cases}$ (59)

∎

######  Proposition 4.

The expected number of tables of the Powered Chinese Restaurant process
evolves with $N\gg 1$ as $H_{\frac{r^{2}+1}{2}}(N)$ for $r<1$ and as
$H_{r}(N)$ when $r\geq 1$, where $H_{m}(n)$ is the generalized harmonic
number.

###### Proof.

In general, the expected number of clusters at the $N^{th}$ step can be
written as:

$\mathbb{E}(K|N,r)=\sum_{1}^{N}\frac{\alpha}{\sum_{k}N_{k}^{r}+\alpha}\stackrel{{\scriptstyle
N^{r}\gg 1}}{{\propto}}\sum_{1}^{N}\frac{1}{\sum_{k}N_{k}^{r}}$ (60)

We showed in Proposition 3 that we can rewrite $\sum_{k}N_{k}^{r}\propto
N^{\frac{r^{2}+1}{2}}$ when $r<1$ and $\sum_{k}N_{k}^{r}\propto N^{r}$ when
$r\geq 1$. Injecting this result in Eq.60 for $r$, we get:

$\mathbb{E}(K|N,r)\stackrel{{\scriptstyle N^{r}\gg
1}}{{\propto}}\begin{cases}\sum_{1}^{N}\frac{1}{N^{\frac{r^{2}+1}{2}}}=H_{\frac{r^{2}+1}{2}}(N)\\\
\sum_{1}^{N}\frac{1}{N^{r}}=H_{r}(N)\end{cases}$ (61)

∎

For $r=1$, $\mathbb{E}(K|N,r=1)\propto H_{1}(N)\approx\gamma+\log(N)$ where
$\gamma$ is the Euler–Mascheroni constant, which is a classical result for the
regular Dirichlet process.

When $r>1$ and $N\rightarrow\infty$, the term $H_{\frac{r^{2}+1}{2}}(N)$
converges towards a finite value and the sum $\sum_{k}p_{k}^{r}$ goes to 1
(see Proposition 2). By definition $\mathbb{E}(K|N,r>1)\stackrel{{\scriptstyle
N\rightarrow\infty}}{{\propto}}\zeta(\frac{r^{2}+1}{2})$, where $\zeta$ is the
Riemann Zeta function.

When $r<1$, we can approximate the harmonic number in a continuous setting. We
rewrite Eq.61 as:

$\begin{split}\mathbb{E}(K|N,r)\stackrel{{\scriptstyle N^{r}\gg
1}}{{\propto}}\sum_{n=1}^{N}\frac{1}{n^{\frac{r^{2}+1}{2}}}&\stackrel{{\scriptstyle
N^{r}\gg 1}}{{\approx}}\int_{1}^{N}n^{-\frac{r^{2}+1}{2}}dn\\\ &\ \ \ \ \ \
=\frac{2}{1-r^{2}}(N^{\frac{1-r^{2}}{2}}-1)\end{split}$ (62)

It can be shown that
$\frac{N^{1-x}-1}{1-x}=H_{x}(N)+\mathcal{O}(\frac{1}{N^{x}})$. Therefore, the
Powered Chinese Restaurant process exhibits a power-law behaviour similar to
the Pitman-Yor process Eq.47 for $r=\sqrt{1-2\beta}$ for $0<r<1$. For values
of $r>1\Leftrightarrow\beta<0$, the equivalent Pitman-Yor process is not
defined unlike the Powered Chinese Restaurant process. Note that there is a
priori no reason for $r$ to be constrained in the domain of real number.
Complex analysis of the process might be an interesting lead for future works.

Figure 26: Numerical validation of Propositions 2 (left), 3 (middle), 4
(right). In the first plot, $K$ is the number of non-empty clusters. In the
second and third plots, the theoretical results are the solid lines, and the
associated numerical results are the transparent lines of same colour. Except
for small $N$, the difference between theory and experiments is almost
indistinguishable.

####  .17.6 Experiments

#####  .17.6.a Numerical validation of propositions

First of all, we present numerical confirmations of propositions stated above
(Propositions 2, 3, 4) by simulating 100 independent Powered Chinese
Restaurant processes with parameter $\alpha=1$ for various values of $r$. Note
that in all the theorems we verify here, $\alpha$ acts as a scaling factor
without modifying the shape of the results discussed here. We present the
results of numerical simulations in Fig.26.

In the left part, we plot the evolution of the probability for each cluster to
be chosen as $N$ grows for $r=0.2$ for one run. We see that the probabilities
do not remain constant but instead diminish as the number of clusters grows.
The figure suggests they all converge to a common value (a uniform
probability) as shown in Proposition 2. The black line shows the probability
of a uniform distribution. We chose not to show the results for $r>1$; in this
case, one probability goes to 1 as the other fades to 0 as $N$ grows, as
expected.

In the middle part of the figure, we plot the expression for
$\sum_{k}N_{k}^{r}$ derived in Proposition 3 (solid lines) versus the value of
the sum from experimental results (transparent lines), averaged over 100 runs.
Note that plots are in a log-log scale and that curves have been shifted
vertically for visualization purposes. As assumed in Proposition 3, the
approximation holds for all values of $r$.

Finally, in the right picture, we plot the evolution of the number of clusters
$K$ versus $N$ according to Proposition 4 (solid lines) and experiments
(transparent lines). The error bars correspond to the standard deviation over
the 100 runs. We see that the expression derived in Proposition 4 accounts
well for the evolution of the number of clusters. Note that plots are in a
log-log scale and that curves have been shifted vertically for visualization
purposes. We must point out that there is a constant shift from experiments to
the theory that does not appear on the plot (because of the rescaling). This
shift comes from the approximation of large $N^{r}$ which is not valid at the
beginning of the process. However, it does not play any role in the evolution
of $K$ as $N$ grows large enough.

#####  .17.6.b Use case: infinite Gaussian mixture model

We now illustrate the usefulness of a prior that alleviates the “rich-get-
richer” property on several synthetic datasets and on a real-world
application. We choose to consider as an illustration its use as a prior in
the infinite Gaussian mixture model. We choose this application to ease visual
understanding of the implications of the PCRP, but the argument holds for
other models using DP priors as well (text modelling, gene expression
clustering, etc.).

Figure 27: Application on synthetic data. (Top) Original datasets used for the
experiments (Density, Diamond, and Grid). (Bottom) Results for various values
of $r$; the x and y axes are all the same. The dashed line indicates the
regular DP prior as $r=1$. The error corresponds to the standard error of the
mean over all runs.

We consider a classical infinite Gaussian mixture model coupled with a Powered
Dirichlet process prior. We fit the data using a standard collapsed Gibbs
sampling algorithm for IGMM [106, 124, 140], with a Normal Inverse Wishart
prior on the Gaussians’ parameters. The input data is shuffled at each
iteration to reduce the ordering bias from the dataset. Note that we cannot
completely get rid of the bias because the Powered Dirichlet Process is not
exchangeable for all $r$. The problem has been addressed on numerous occasions
(Uniform process [124], distance-dependent CRP [17, 40], spectral CRP [113])
and shown to induce negligible variations of results in the case of Gibbs
sampling. We stop the sampler once the likelihood of the model reaches
stability ; we repeat this procedure 100 times for each value of $r$. Finally,
the parameter $\alpha$ is set to 1 in all experiments (see Section .17.2).

Note that we choose not to compare to other types of clustering algorithms. In
this section, we demonstrate the power of alternative forms of the Dirichlet
process. The argument on a simple model (here a regular DP combined with IGMM)
extends to other priors built on Dirichlet processes (Hierarchical and Nested
Dirichlet processes). Besides, comparison of DP-based priors to other
clustering methods (KNN, DBScan, Spectral clustering, etc.) has already been
done numerous times and is out of this section’s scope.

######  Synthetic data

Synthetic datasets are represented in Fig.27-top and comprise $N$=1000
observations each. They have been generated by sampling from 2D Gaussian
distributions whose relative parameters are explicit from Fig.27.

Table 10: Numerical results of the Powered Dirichlet process, Uniform process and Dirichlet process priors coupled to a standard Infinite Gaussian Mixture Model, for the 3 synthetic datasets plotted Fig. 27 and 2 real-world datasets (100 runs each). We see that using PDP as prior makes the model outperform the baselines consistently on every metric. The standard error is given in shorthand form – $0.123(12)\Leftrightarrow 0.123\pm 0.012$.  |  |  Adj.MI |  Adj.RI |  Norm.VI |  $\frac{K_{inf}-K_{true}}{K_{true}}$ |
---|---|---|---|---|---|---
Density |  PDP (r=0.60) |  0.992(1) |  0.980(2) |  0.006(1) |  0.045(5) |
DP |  0.951(4) |  0.797(17) |  0.037(3) |  0.128(10) |
UP |  0.939(2) |  0.854(4) |  0.050(1) |  0.548(1) |
Diamond |  PDP (r=0.50) |  0.982(2) |  0.956(5) |  0.011(1) |  0.063(7) |
DP |  0.909(7) |  0.731(19) |  0.053(4) |  0.202(12) |
UP |  0.927(2) |  0.844(6) |  0.051(2) |  0.544(2) |
Grid |  PDP (r=0.85) |  0.997(1) |  0.990(2) |  0.003(1) |  0.014(2) |
DP |  0.995(1) |  0.977(4) |  0.004(1) |  0.018(3) |
UP |  0.811(1) |  0.517(3) |  0.154(1) |  2.120(1) |
Iris |  PDP (r=0.90) |  0.868(4) |  0.866(7) |  0.057(2) |  0.000(0) |
DP |  0.843(6) |  0.820(12) |  0.065(2) |  0.030(10) |
UP |  0.544(2) |  0.295(3) |  0.303(2) |  2.777(32) |
Wines |  PDP (r=0.10) |  0.712(15) |  0.637(20) |  0.102(5) |  0.157(17) |
DP |  0.589(19) |  0.461(16) |  0.128(4) |  0.327(13) |
UP |  0.713(17) |  0.657(21) |  0.103(5) |  0.147(17) |
Cancer |  PDP (r=0.10) |  0.254(17) |  0.278(21) |  0.118(1) |  0.000(0) |
DP |  0.085(16) |  0.094(19) |  0.108(2) |  0.000(0) |
UP |  0.271(17) |  0.300(21) |  0.118(1) |  0.000(0) |
20-NG |  PDP (r=0.80) |  0.421(4) |  0.119(3) |  0.477(3) |  - |
DP |  0.404(4) |  0.105(4) |  0.491(3) |  - |
UP |  0.000(4) |  0.000(0) |  0.830(3) |  - |

We present the results on synthetic data in Fig.27-bottom and in Table10. We
consider standard metrics in clustering evaluation with a non-fixed number of
clusters: mutual information score and rand index both adjusted for chance
(Adj.MI and Adj.RI), normalized variation of information (Norm.VI, lower is
better), Fowlkes-Mallow score, marginal likelihood (normalized for
visualization) and absolute relative variation of the inferred number of
clusters according to the number used to generate the data
($\mathbf{\frac{K_{inf}-K_{true}}{K_{true}}}$, lower is better). Note that we
purposely chose stereotypical cases to illustrate the argument better. The
Density dataset on the left of Fig.27 is informative about the change induced
by $r$. Here, clusters are distributed at various scales in the dataset; we
see that the lower the value of $r$, the better the results. Indeed, when $r$
is small, the model can distinguish clusters in the dense area better, whereas
when $r$ is closer to 1, the clusters in the dense area are put together in a
larger cluster. The same happens with the Diamond dataset in the middle of
Fig.27, where clusters are distributed according to two different scales.
Finally, on the Grid dataset on the right part of Fig.27, we see an optimum
$r$ exists to distinguish the clusters distributed on a grid; it makes sense
since only one scale in clusters distribution is involved in this dataset. In
Table. 10, we explicitly report the values of the PDP optimal $r$ and compare
them to the values yielded by the same model using either a Dirichlet Process
(DP) prior or a Uniform Process (UP) prior.

Figure 28: Application to spatial clustering on geolocated data for $r=0.8$
(left), $r=1$ (right) and $r=0.9$ (middle). We see that the model using a
Powered Dirichlet process prior for $r=0.9$ and $r=0.8$ describes the data
better than the same model using a Dirichlet process prior ($r=1$).

######  Real-world data

In Table 10, we report the results for the 20Newsgroup (20-NG) real-world
dataset, which is a collection of 18,000 user posts published on Usenet,
organized in 20 Newsgroup (which are our target thematic clusters). As a
model, we consider a modified version of LDA [19] that uses a PDP prior
instead of DP in the words sampling step. Note that because the number of
clusters must be provided to LDA, we do not compute
$\mathbf{\frac{K_{inf}-K_{true}}{K_{true}}}$. We also run an additional
experiment on three additional well-known datasets: Iris (4 attributes, 3
classes), Wines (13 attributes, 3 classes) and Cancer (30 attributes, 2
classes). We see PDP yields improved performances on every dataset.

We now illustrate the interest of using an alternate form of prior for the
Infinite Gaussian Mixture model on real-world data. We consider a dataset of
4.300 roman sepulchral inscriptions comprising the substring “Antoni” that
have been dated between 150AC and 200AC and assigned with map coordinates. The
dates correspond to the reign of Antoninus Pius over the Roman empire. The
dataset is available on Clauss-Slaby repository [26]. It was common to give
children or slaves the name of the emperor; the dataset gives a global idea of
the main areas of the roman empire at that time [53]. The task is to discover
spatial clusters of individuals named after the emperor. We expect to find
geographical clusters around: Italy, Egypt, Gauls, Judea, and all along the
limes (borders of the roman empire, which concentrate lots of sepulchral
inscriptions for war-related reasons) [52]. We present the results in Fig.28.

We see that when $r=1$, the classical DP prior is not fit for describing this
dataset, as it misses most of the clusters. On the other hand, when $r=0.9$,
the infinite Gaussian mixture model retrieves the expected clusters. It also
makes some clusters that were not expected, such as the north Italian cluster
or the long cluster going through Spain and France that corresponds to the
layout of the Roman roads (via Augusta and via Agrippa; it was common to bury
the dead on roads edges). Finally, when $r=0.8$, we get even more detail: some
of the main clusters are broken into smaller ones (Italy breaks into Rome,
North Italy, and South Italy; Britain becomes an independent cluster, etc.).
In this case, changing $r$ controls the level of details of the clustering.
Note that we do not compute metrics for this experiment in the absence of
“ground-truth” clustering; there is no such thing as the right clustering in
this case. Applied to spatial data, the PDP prior allows to control the
clustering granularity. We see how different results can be according to the
extent the model relies on the “rich-get-richer” prior and how its control is
needed to make modelling relevant to a given situation.

####  .17.7 Conclusion

In this section, we discussed the necessity of controlling the “rich-get-
richer” property that arises from the classical Chinese Restaurant Process
usual formulation. We showed cases where this modelling hypothesis must be
alleviated or strengthened to describe data more accurately. To this end, we
derive the Powered Chinese Restaurant Process from a powered version of the
Dirichlet-Multinomial distribution. This formulation allows reducing the
expected number of clusters, which is not possible in the standard Pitman-Yor
processes, while generalizing the standard Dirichlet process and the Uniform
process.

The main feature of this formulation is that it allows for direct control of
the “rich-get-richer” priors’ importance. We derive elementary results on
convergence and the expected number of clusters of the new process.

Finally, we show that it yields better results on synthetic data when coupled
to a standard Gaussian mixture model and illustrates a possible use case with
real-world data. For future works, it might be interesting to investigate
cases where $r$ takes non-positive values (which might lead to a “poor-get-
richer” kind of process), and to develop a procedure to infer it automatically
for specific problems (by minimizing a dispersion criterion for instance).

The regular Chinese Restaurant process has been used for decades as a prior in
many real-world applications. However, alternate forms for this prior have
been little explored. We are convinced that controlling the impact of the
“rich-get-richer” hypothesis will bring significant changes in many state-of-
the-art models.

We now propose to illustrate this claim and to further answer our task at
hand. To this end, we use PDP to redefine the Dirichlet-Hawkes Process in
terms of the Powered Dirichlet-Hawkes Process (PDHP) and investigate the new
possibilities brought by this reformulation.

###  .18 Powered Dirichlet-Hawkes Process – Modelling self interacting
clusters

This work has been published, see [99]

####  .18.1 Introduction

#####  .18.1.a PDHP as an answer to DHP’s limits

In this section, we develop the Powered Dirichlet-Hawkes process (PDHP) as an
answer to the limits of DHP raised in Section .16. It makes use of the
redefinition of Dirichlet Processes detailed in Section .17.

We highlighted earlier that DHP has several limits. It fails to model data
correctly when textual information is scarce (e.g., short texts such as
tweets, or when topics’ vocabulary overlap) [141] and when temporal
information is scarce (overlapping dynamics, few observations). In addition,
the usual assumption is that publication dynamics and textual content are
perfectly correlated, which cannot be true in real-world applications. For
instance, two identical textual contents will trigger different publication
dynamics depending on who published them, or when they have been published.

#####  .18.1.b Contributions

We overcome all these limitations by developing the Powered Dirichlet-Hawkes
Process (PDHP), which yields better results than DHP on every dataset
considered (up to +0.3 NMI). In particular, when textual or temporal
information is scarce, PDHP allows to control the importance given to one or
the other and thus retrieves better clusters. PDHP is the DHP equivalent of
controlling the importance of the “rich-get-richer” hypothesis in DP using
PDP, only using temporal information as a prior. PDHP also allows to
distinguish textual clusters from temporal clusters (documents that follow the
same dynamic independently from their content).

Our contributions are listed below:

  * •

We highlight and explain the limitations of the DHP prior: it does not handle
weakly informative temporal and textual information and it is not designed to
consider different dynamics between text and time.

  * •

We derive the Powered Dirichlet-Hawkes process (PDHP) as a new prior in
Bayesian non-parametric for the temporal clustering of a stream of textual
documents, which is a generalization of the Dirichlet-Hawkes process (DHP) and
of the Uniform process (UP).

  * •

We show how the PDHP prior performs better than DHP and UP priors through a
thorough evaluation and comparison on several synthetic datasets and real-
world datasets from Reddit.

  * •

We show that PDHP prior allows choosing whether clusters are based more on
textual or temporal information, or a mixture of both. We can favour their
generation more or less according to documents’ textual content or their
publication dynamics.

  * •

We perform a detailed analysis of PDHP’s results on real-world datasets from
Reddit. We illustrate our approach with examples of real-world topics
uncovered by our method. We systematically analyse the influence of the
hyperparameter $r$ introduced in PDHP. Our method allows us to recover more or
less bursty events from data streams depending on $r$.

####  .18.2 Model and algorithm

#####  .18.2.a Dirichlet prior and alternatives

We briefly recall the definition of a Dirichlet prior (see Section .17 for an
extensive discussion on DP). A Dirichlet prior for clustering implements the
assumption that the more a cluster is populated, the more chances a new
observation belongs to it (“rich-get-richer” property). Besides, there is
still a chance that a new observation gets assigned to a newly created
cluster. It is often expressed using a metaphor, the Chinese Restaurant
Process (CRP), and it goes as follows: if an $i^{th}$ client arrives in a
Chinese restaurant, they will sit at one of the $K$ already occupied tables
with a probability proportional to the number of persons already sat at this
table. They can also sit alone at a new table $K+1$ with a probability
inversely proportional to the total number of clients in the restaurant. When
their choice is made, the next client arrives, and the process is repeated.
Let $c$ be the cluster chosen by the $i^{th}$ customer, $\vec{C^{-}}$ the
table assignment of previous customers up to $i-1$, $N_{c}$ the population of
table $c$, $C$ the number of already occupied tables and
$\alpha_{0}\in\mathbb{R}^{+}$ the concentration parameter. The process can be
written formally as:

$\text{CRP}(C_{i}=c|\vec{C^{-}},\alpha_{0})=\begin{cases}\frac{N_{c}}{\alpha_{0}+N}\text{
if c = 1, 2, ..., C}\\\ \frac{\alpha_{0}}{\alpha_{0}+N}\text{ if c =
C+1}\end{cases}$ (63)

The Uniform process [124] has been proposed as an alternative to the DP prior.
In this context, a new customer entering the restaurant has an identical
chance to sit at either of the occupied tables, and a chance to sit at an
empty table inversely proportional to the number of occupied tables. Formally:

$\text{U-CRP}(C_{i}=c|\vec{C^{-}},\alpha_{0})=\begin{cases}\frac{1}{\alpha_{0}+C}\text{
if c = 1, 2, ..., C}\\\ \frac{\alpha_{0}}{\alpha_{0}+C}\text{ if c =
C+1}\end{cases}$ (64)

Finally, the Powered Dirichlet process that we detailed Section .17
generalizes the two processes above, stating that the probability for a new
client to sit at a new table depends arbitrarily on the number of customers
already sat at this table:

$\text{P-CRP}(C_{i}=c|r,\vec{C^{-}},\alpha_{0})=\begin{cases}\frac{N_{c}^{r}}{\alpha_{0}+\sum_{c^{\prime}}N_{c^{\prime}}^{r}}\text{
if c = 1, 2, ..., C}\\\
\frac{\alpha_{0}}{\alpha_{0}+\sum_{c^{\prime}}N_{c^{\prime}}^{r}}\text{ if c =
C+1}\end{cases}$ (65)

where $r\in\mathbb{R}^{+}$ is a hyper-parameter. Varying $r$ allows to give
more or less importance to the “rich-get-richer” hypothesis of DP. Note that
$P-CRP(r=0,\vec{C^{-}},\alpha_{0})=U-CRP(\vec{C^{-}},\alpha_{0})$ and that
$P-CRP(r=1,\vec{C^{-}},\alpha_{0})=CRP(\vec{C^{-}},\alpha_{0})$. We will use
this more general form in the rest of this section and make $r$ vary to
compare those priors in the experimental section.

#####  .18.2.b Hawkes processes

The Hawkes process has already been introduced and discussed in Section
.16.2.b. We recall that Hawkes processes are defined as self-stimulating
temporal point processes. They are fully characterized by their intensity
function $\lambda(t)$, which is related to the probability
$P(t_{events}\in[t;t+\Delta t])$ of an event happening between $t$ and
$t+\Delta t$ by $\lambda(t)=\lim_{\Delta t\rightarrow
0}\frac{P(t_{events}\in[t;t+\Delta t])}{\Delta t}$.

Same as in [30], the intensity of the Hawkes process associated with cluster
$c$ is defined as:

$\lambda_{c}(t|\mathcal{H}_{<t,c})=\sum_{\mathcal{H}_{<t,c}}\vec{\alpha_{c}}^{T}\cdot\vec{\kappa}(t_{i,c})$
(66)

We recall that $\vec{\alpha_{c}}$ is a vector of weights and $\vec{\kappa}(t)$
is a vector of user-defined kernel functions with the same dimension as
$\vec{\alpha}$. The kernel functions are defined at the beginning of the
algorithm and are not modified afterwards. We will want to infer the weights
vector $\vec{\alpha}$ to determine which entries of the kernel vector are the
most relevant for a given situation.

#####  .18.2.c Powered Dirichlet-Hawkes process

In [30], the authors create DHP by substituting the counts $N_{k}$ in the DP
with the intensities of the Hawkes processes associated with each cluster.
Instead, we simply substitute the counts in PDP with the intensities of the
Hawkes processes associated with each cluster, forming the Powered Dirichlet-
Hawkes Process (PDHP):

$P(C_{i}=c|t_{i},r,\lambda_{0},\mathcal{H}_{<t_{i},c})=\begin{cases}\frac{\lambda_{c}^{r}(t_{i})}{\lambda_{0}+\sum_{c^{\prime}}\lambda_{c^{\prime}}^{r}(t_{i})}\text{
if c$\leq$C}\\\
\frac{\lambda_{0}}{\lambda_{0}+\sum_{c^{\prime}}\lambda_{c^{\prime}}^{r}(t_{i})}\text{
if c=C+1}\end{cases}$ (67)

where $t_{i}$ is the arrival time of document $i$. We reformulate the
Dirichlet-Hawkes process to allow nonlinear dependence ($r$) on the non-
integer counts ($\vec{\lambda}$). We illustrate how $r$ influences the prior
probability of cluster selection in Fig. 29.

Figure 29: Temporal prior on a new observation cluster — (Top) Prior
probability for two clusters to get chosen at all times according to DHP
(dotted lines) (Bottom) Same prior probability for each cluster as a function
of $r$ in the PDHP at a given time.

#####  .18.2.d Textual modelling

To compare to DHP on solid ground, we consider the same textual model the
authors used in the original paper, the Dirichlet-Multinomial model. This
model is detailed in Section .16.2.d and it considers the textual content as a
bag of words.

The likelihood of the $i^{th}$ document belonging to cluster $c$ reads:

$\begin{split}\mathcal{L}(C_{i}=c|N_{<i,c},n_{i},\theta_{0})&=P(n_{i}|C_{i}=c,N_{<i,c},\theta_{0})\\\
&=\frac{\Gamma(N_{c}+\theta_{0})}{\Gamma(N_{c}+n_{i}+\theta_{0})}\prod_{v}\frac{\Gamma(N_{c,v}+n_{i,v}+\theta_{0,v})}{\Gamma(N_{c,v}+\theta_{0})}\end{split}$
(68)

where $N_{c}$ is the total number of words in cluster $c$ from observations
previous to $i$, $n_{i}$ is the total number of words in document $i$,
$N_{c,v}$ the count of word $v$ in cluster $c$, $n_{i,v}$ the count of word
$v$ in document $i$ and $\theta_{0}=\sum_{v}\theta_{0,v}$.

#####  .18.2.e Posterior distribution

The resulting posterior distribution of the $i^{th}$ document over clusters is
calculated using Bayes theorem. It is proportional to the product of the
textual likelihood (Eq.68) and the temporal Powered Dirichlet-Hawkes prior
(Eq.67):

$\begin{split}&P(C_{i}=c|r,n_{i},t_{i},N_{c},\mathcal{H}_{<t,c})\\\
\propto&\underbrace{P(n_{i}|C_{i}=c,N_{<i,c},\theta_{0})}_{\text{Textual
likelihood}}\underbrace{P(C_{i}=c|t_{i},r,\lambda_{0},\mathcal{H}_{<t_{i},c})}_{\text{Temporal
prior}}\\\
=&\frac{\Gamma(N_{c}+\theta_{0})}{\Gamma(N_{c}+n_{i}+\theta_{0})}\prod_{v}\frac{\Gamma(N_{c,v}+n_{i,v}+\theta_{0,v})}{\Gamma(N_{c,v}+\theta_{0})}\\\
&\times\begin{cases}\frac{\lambda_{c}^{r}(t_{i})}{\lambda_{0}+\sum_{c^{\prime}}\lambda_{c^{\prime}}^{r}(t_{i})}\text{
if c = 1, ..., C}\\\
\frac{\lambda_{0}}{\lambda_{0}+\sum_{c^{\prime}}\lambda_{c^{\prime}}^{r}(t_{i})}\text{
if c = C+1}\end{cases}\end{split}$ (69)

We recall that $\lambda_{c}(t)$ is defined Eq. 66. The textual likelihood to
open a new cluster $C+1$ is computed by setting $N_{C+1,v}=0$ – because it is
empty before being opened.

#####  .18.2.f Algorithm for parameters inference

We use a similar algorithm to the one in [30]. Briefly, the algorithm is a
sequential Monte-Carlo (SMC) that takes one document at a time in their order
of arrival. The algorithm starts with a number $N_{part}$ of particles whose
weights are $\omega_{p}=\frac{1}{N_{part}}$, each of which will keep track of
a hypothesis on documents clusters. After a few iterations, particles that
contained unlikely allocation hypotheses are discarded and replaced by more
likely ones. The likeliness of a hypothesis is encoded in the weights of each
particle $\omega_{p}$. Such genetic algorithms are favoured for optimizing
DHP-based models. These models are not convex, and genetic algorithms allow to
tackle the problem sequentially.

Figure 30: Schematic workflow of the SMC algorithm — For each new observation
from a stream of document, we run steps 1 (sample document’s cluster), 2
(update sampled cluster’s internal dynamics) and 3 (update particle
hypothesis’ likeliness) for each particle, and then discard particles
containing the less likely hypothesis on cluster allocation.

For each particle, when a new document arrives, (1) the cluster of the
document is sampled according to a Categorical distribution over all clusters,
whose weights are determined by Eq. 69. After the cluster of the new document
has been sampled, (2) the kernel weights $\vec{\alpha}$ from Eq. 66 are
updated using Eq. 42. For efficiency purpose, we sample $\vec{\alpha}$ from a
set of $N_{s}$ pre-computed $\vec{\alpha}$ vectors. We finally (3) update the
weights $\omega_{p}$ of each particle according to the posterior Eq. 69 such
as $\omega_{p}^{(n+1)}=\omega_{p}^{(n)}\times\text{Eq.\leavevmode\nobreak\
\ref{eq-likModelTot}}$. If the weight of a particle falls below a value
$\omega_{thres}$, the particle is discarded and replaced by another existing
one with sufficient weight. The full process is illustrated in Fig. 30. By
updating incrementally the likelihood associated with each of the pre-computed
$\vec{\alpha}$ sample vectors, the algorithm processes each new observation in
constant time $\mathcal{O}(1)$.

The task of updating kernel coefficients (2) is the same as in any Hawkes
process, and the task of updating particle weights and resampling them (3) is
common to any SMC algorithm. The change induced by the PDHP compared to the
DHP happens in step (1). First of all, we note that for $r=1$ the PDHP prior
is identical to the DHP prior. From Section .17, lowering the value of $r$
reduces the “rich-get-richer” aspect of the PDP (“rich-get-less-richer”),
whereas increasing it leads to a “rich-get-more-richer” effect. These
metaphors can be translated as follows in our temporal context: for lower
values of $r$, the relative difference between cluster’s temporal intensities
plays a less significant role in cluster selection, whereas higher values of
$r$ tend to exacerbate these differences and make the temporal aspect of the
greatest consequence on the choice of a cluster. In other words, tuning the
value of $r$ allows giving more or less importance to the temporal aspect of
the clustering. This is illustrated in Fig. 31.

In Fig. 29, we plot the situation when a new observation gets assigned a
cluster. The associated Hawkes intensities are the base to compute the prior
probability for either cluster. This quantity is then modulated using $r$ to
give more or less importance to intensity differences between clusters. In
Fig. 31, we plot the probability for various clusters to be chosen (which is
directly proportional to the posterior distribution, see Eq. 69) according to
$r$ when their textual likelihood and Hawkes process intensity are known. Note
that for $r=0$, the probability for any cluster to get chosen is linearly
proportional to its textual likelihood (Dirichlet-Uniform process), whereas
when $r$ increases, the probability of getting chosen gets closer to a
selection based on the temporal aspect only.

Figure 31: Effect of $r$ on cluster selection probabilities — The probability
for each cluster to get chosen (solid lines) for several values of $r$ and
fixed individual textual likelihood (blue bars) and Hawkes intensity (orange
bars).

This makes the main interest of the PDHP model. Tuning the parameter $r$
allows choosing whether inferred clusters are based on textual or temporal
considerations. It generalizes several state-of-the-art works, which are
special cases of the PDHP for different values of $r$. The DHP [30] is
equivalent to PDHP for $r=1$; the UP [124] is equivalent to PDHP when $r=0$.
In the following sections, we show how fine-tuning $r$ systematically yields
significantly better results than setting it to $r=0$ or $r=1$ (up to a gain
of 0.3 on our experiments’ normalized mutual information metric). We also show
how varying it allows to recover one kind of clustering or the other (textual
or temporal) with high accuracy and see how it affects clustering results on
several real-world datasets.

####  .18.3 Experiments

#####  .18.3.a Synthetic data

######  Synthetic data generation

We simulate a case where only two clusters are considered. Each cluster has
its own vocabulary distribution over 1,000 words and its own kernel weights
$\vec{\alpha}$, with Gaussian Hawkes kernel functions $\vec{\kappa(t)}$ of
parameters $(\mu,\sigma)$=(3, 0.5), (7, 0.5) and (11, 0.5) (see Eq. 66).
Finally, we set $\lambda_{0}=0.05$. We first simulate one independent Hawkes
process per cluster using the Tick Python library [7]. The processes are
stopped at time $t=1500$, which makes a rough average of 7,000 events per run.
Then we associate each simulated observation with a sample of 20 words drawn
from the corresponding cluster’s word distribution. The inference has been
performed using an 8 core processor (i7-7700HQ) with 8GB of RAM on a laptop,
which underlines how scalable the algorithm is. As stated before, the
algorithm processes each new document in constant time $\mathcal{O}(1)$, which
ranged from 0.05s on synthetic data to maximum 1s on real-world data. Note
that this number is directly proportional to the number of active inferred
clusters, and thus depends strongly on the dataset.

Figure 32: Overlaps — (Left) Temporal overlap is defined as the ratio between
the area common to two Hawkes intensities and the total area under the
intensity functions. (Right) Textual overlap is defined as the proportion of
vocabulary that is common to two clusters, weighted by the probability of
words within their respective cluster.

We generate ten such datasets for every considered value of vocabulary overlap
and Hawkes intensities overlap, which leave us with 200 datasets (5 values of
textual overlap x 4 values of temporal overlap x 10 datasets). Overlap is
defined as the shared area of two distributions, normalized by the total area
under the distributions. For instance, if the vocabulary of one cluster ranges
from words "1" to "100" with uniform distribution, and the vocabulary of
another cluster from words "50" to "150" with uniform distribution, the
overlap equals 50%. We define the overlap of Hawkes process intensity in the
same way. If the triggering Hawkes kernel of one cluster is a Gaussian
function with $(\mu,\sigma)=(3,1)$ and one associated observation at $t=0$,
and the triggering kernel of the other is also a Gaussian function but with
$(\mu,\sigma)=(5,1)$ also with an associated observation at $t=0$, the overlap
equals 32% (see Fig. 32). When computing the Hawkes intensity overlap, every
observation within a cluster and its associated timestamp are considered. The
definition of overlaps is illustrated in Fig. 32. To enforce a given
vocabulary overlap (Fig. 32-right), we shift the word distributions of the
clusters from which events’ vocabulary is sampled. To enforce a given Hawkes
intensities overlap (Fig. 32-left), we shift the event times of every event in
one of the clusters until we get the correct overlap ($\pm 5\%$).

Note that we consider ten different datasets instead of considering ten runs
per dataset for two reasons. Firstly, the generation of Hawkes processes is
highly stochastic, so a model might perform significantly better on a single
dataset only by chance. Secondly, given the way the SMC algorithm works, the
standard deviation between runs is small: at each iteration, $N_{part}$
clustering hypotheses are tested, which is equivalent to running $N_{part}$
times a single clustering algorithm. We heuristically set $N_{part}=8$, as we
observe no significant improvement using more particles.

The other parameters we use for clustering synthetic data are:
$\alpha_{0}=0.1$, ${\theta_{0}=1}$,
$\vec{\kappa(t)}=[\mathcal{G}(t;3,0.5),\mathcal{G}(t;7,0.5),\mathcal{G}(t;11,0.5)]$
with $\mathcal{G}(t;\mu,\sigma)$ the Gaussian function, $N_{samples}=2.000$
and $\omega_{thres}=\frac{1}{2N_{part}}$.

Figure 33: PDHP yields good NMI values — Normalized mutual information (NMI)
for various values of $r$, intensities overlap and vocabulary overlap, for one
dataset per combination. The results for $r=0$ are the output of the Uniform
process, the results for $r=1$ are the output of the DHP [30], and the other
values of $r$ correspond to other particular cases of PDHP. The darker the
better. Overall, PDHP yields good NMI values (the maximum being 1).

We are interested in varying both vocabulary and intensities overlap to
exhibit the limits of DHP and how PDHP overcomes them. Note that in the
synthetic data experiments in [30] (Figs.3a and 3b), the intensities overlap
is almost null, which makes the task easier for the Hawkes part of the
algorithm. The primary metric we use throughout the experimental section is
the normalized mutual information (NMI). During the experiments, we also
considered the Adjusted normalized rand index and the V-measure, which are
adapted to evaluate clustering results when the number of inferred clusters is
different from the true number of clusters. The observed trends in results
from these other metrics are identical to the ones observed for NMI.
Therefore, we choose to report only the results of the latter for clarity.

The NMI metric is standard when evaluating non-parametric clustering models.
It compares two cluster partitions (i.e., the inferred and the ground truth
ones); it is bounded between 0 (each true cluster is represented to the same
extent in each of the inferred ones) and 1 (each inferred partition comprises
100% of one true cluster).

######  PDHP yields better results as vocabulary overlap increases

We report our results when the intensities overlap is null, with varying $r$
and the vocabulary overlap in Fig. 34a. Because we consider ten different
datasets for each set of overlap parameters, it makes no sense to report the
absolute average NMI since it can vary greatly from one dataset to the other.
Instead, we plot the relative NMI difference between PDHP and DHP ($r=1$),
which we expect to be less dependent on the datasets we consider. However, to
give an idea of the typical performance for some parameters, we also provide
raw results for one run in Fig. 33.

Figure 34: PDHP performs better than DHP — Difference between the normalized
mutual information (NMI) of PDHP and DHP model [30] for various values of $r$,
intensities overlap and vocabulary overlap, averaged over all the datasets.
Red means PDHP performed better, blue means PDHP performed less well. Because
PDHP($r=1$)=DHP, the column $r=1$ show no difference. PDHP allows to increase
results on NMI by as much as 0.3 over DHP.

There is a clear correlation between efficiency, vocabulary overlap and $r$,
with a gain on NMI up to $+30\%$ of its maximal value over DHP. As stated at
the end of the "Model" section, this result was expected: the more vocabulary
overlap grows, the less textual content carries valuable information to
cluster the documents. This observation supports the concerns raised in [141]
about the efficiency of DHP for clustering short text documents. However,
Hawkes intensities overlap being null, the arrival time of events carries
highly valuable information when textual content does not allow to distinguish
clusters well. Therefore, PDHP provides a way to tackle the problem raised in
[141] without the need to sample observations.

Conversely, when vocabulary overlap is null, the textual content provides
enough information to distinguish clusters correctly. The temporal dimension
only allows refining the results with no significant improvement for all
values of $r$.

Finally, we can see how the Dirichlet-Uniform process (DUP, $r=0$)
consistently yields worse performances under these settings. Once again this
is expected, since in this synthetic experiment intensities overlap carries
valuable information about events clustering; DUP only considers textual
information and therefore misses valuable clues.

######  PDHP yields similar results for null vocabulary overlap

We report comparable results in Fig. 34b. Here, we consider a null vocabulary
overlap for various values of $r$ and of Hawkes intensities overlap. The
situation is now the opposite: the textual content always carries valuable
information about clusters, whereas the temporal aspect does not. We observe
the same trend as in Fig. 34a –note that the colour scale is the same. Varying
the value of $r$ does not significantly change the performance of clustering,
meaning the textual content always carries enough information. This plot shows
that PDHP can handle greater intensities overlap without collapsing into
unrealistic clustering. Since in most real-case applications, many clusters
with various dynamics may coexist simultaneously, it is comforting that the
PDHP can also handle this case.

######  PDHP yields better results in more realistic situations

We finally report the results for intermediate values of intensities and
vocabulary overlaps in Fig. 34c,d. In real-world applications, it seldom
happens that topics’ vocabularies do not overlap at all. For instance, a quick
analysis of The Gutenberg Webster’s Unabridged Dictionary by Project Gutenberg
shows that 22% of English words are associated with more than one definition.
A more detailed analysis would need to consider the usage frequency of words
to get correct statistics. Still, this number provides an estimate of the
effective vocabulary overlap in real-world situations.

In Fig. 34c, we present the results for a fixed intensities overlap of 0.5
versus various values of $r$ and vocabulary overlaps, and in Fig. 34d for a
fixed vocabulary overlap of 0.5 versus various values of $r$ and intensities
overlaps. Once again, we see that, on average, using PDHP can increase the NMI
over DHP up to +20% of the maximum possible value.

######  PDHP finds textual or temporal clusters depending on r

Figure 35: Textual (orange) and temporal (blue) NMI vs r when textual and
temporal clusters are decorrelated — From top-left to bottom-right, there are
10%, 20%, 30%, 40%, 50%, 60%, 70%, 80% and 90% of generated events that have
been randomly re-assigned a textual cluster. The orange curves are the textual
NMI vs $r$, which evaluate how well events whose vocabulary has been sampled
from the same distribution are clustered together; the blue curves are the
temporal NMI vs $r$, which evaluate how well events following the same
temporal dynamic are correctly clustered together. The values presented are
for one dataset. We see that varying $r$ allows retrieving the right temporal
($r$ large) or textual clusters ($r$ small).

We now slightly modify our experimental setup. Instead of considering that
textual clusters and Hawkes intensities are perfectly correlated, we consider
a decorrelated case. A document whose vocabulary is drawn from cluster $C_{1}$
can now follow the same temporal dynamics as cluster $C_{2}$. If we imagine a
dataset of news articles published online, it is clear why this might happen
frequently. If popular newspapers such as New York Times or Reuters publish an
article on a topic $A$ at a time t, it is likely to trigger snowball
publications of related articles from less popular journals. “Popularity” is
chosen as an indicator in this example, but it may be any other external
parameter (centrality in news networks, support of publications, etc.). In
this case, the article’s textual content allows to uncover a “story of
publication”, that is, how the article has been spread, when publication
spikes are, etc. However, the temporal information would help understand the
dynamics of publications interaction: which reduced set of articles triggered
the publication of subsequent ones.

In [30], it is assumed that every document within clusters follows a unique
dynamic. We relax this hypothesis in our datasets as follows. For null textual
and temporal overlaps, after a dataset has been generated, we resample the
textual clusters of a fraction of randomly selected events, as well as the
words associated with the event. In doing so, we decorrelate temporal and
textual clusters. Therefore, an event is now described by two cluster
indicators: its temporal cluster (which Hawkes intensity made the event appear
where it is) and a textual cluster (which vocabulary has been used to sample
the words the event contains).

For completeness, we also show the results for various decorrelations for one
run in Fig. 35. To better understand the tendency of NMIs with respect to $r$,
we plot the average difference between the NMI of textual clustering and the
NMI temporal clustering over all the datasets. Explicitly: $\Delta
NMI=NMI_{text}-NMI_{temp}$. The results are reported in Fig. 36.

Figure 36: Varying $r$ allows to choose between textual or temporal clustering
— The black line plots the difference between the NMI of textual and temporal
clustering. For small $r$, textual clustering is far better than temporal
clustering, and for large $r$, the situation is reversed. This is because $r$
determines the importance given to the temporal dimension and therefore allows
choosing between retrieving temporal or textual clusters.

As supposed at the end of the “Model” section, varying $r$ allows retrieving
one clustering or the other. Note that the value $r$ of transition from text
to time clustering depends directly on the dataset considered: number of words
sampled, vocabulary size, overlaps, etc.

######  PDHP efficiently infers the temporal dynamics of each cluster

Finally, we show that PDHP correctly infers kernels’ parameters in every
situation where events are correctly assigned to their temporal cluster. The
results are reported in Fig. 37. We looked at the mean absolute error (MAE)
and the mean Jensen-Shannon divergence (MJS) between the vector $\vec{\alpha}$
used to generate the dataset and the inferred one. We note in Fig. 37 that
when the textual overlap is small, the inferred kernel is close to the real
one and $r$ has a minor impact on the result. This is because the inferred
kernels mostly depend on the correctness of inferred clusters: when
observations are allocated to the right clusters, the Hawkes process inference
considers relevant information when inferring these clusters’ dynamics.
However, when observations are misallocated, the Hawkes process infers
dynamics also based on times that are not supposed to contribute to this
cluster’s dynamic. When the clustering task is simple and yields good results
(that is, when the textual overlap is small, see Fig. 33), the PDHP infers
correct temporal dynamics ($\sim 5\%$ MAE); this shows our method correctly
accounts for clusters dynamics given the available information.

When vocabulary overlap is large, the value of $r$ significantly influences
the kernel inference performances. However, when $r$ is chosen so that
clusters are correctly inferred, the kernel inference retrieves well the
expected kernels ($\sim 5\%$ MAE). Finally, the temporal kernel inference is
expected to be less precise when temporal overlap increases, which is what
happens in Fig. 37-bottom-right. In this case, the model does not retrieve
well the synthetic kernels even for the optimal $r$. Besides, the error bars
get wider as a consequence of the cluster allocation being more challenging.
Overall, provided the right clusters, we conclude that our method correctly
retrieves the inferred temporal kernels.

Figure 37: Varying $r$ allows to better capture the dynamics at stake — We
plot the mean average error and the mean Jensen-Shannon divergence of the
inferred kernel function $\alpha$ with respect to the kernel used to generate
the data, for various values of temporal and textual overlaps. The standard
error bars are computed over 10 independent runs. The higher the temporal
overlap, the larger the error bars; the larger the textual overlap the more
influence has $r$.

#####  .18.3.b Real-world application on Reddit

We use the PDHP prior to model real streams of textual documents. We consider
three Reddit datasets [12] about different topics. The News dataset is made of
73.000 titles extracted from the subreddits inthenews, neutralnews, news,
nottheonion, offbeat, open_news, qualitynews, truenews and worldnews, from
April 2019. We chose this month because of the wide variety of events that
happened then (for instance, Sri Lanka Easter bombings, Julian Assange arrest,
first direct picture of a black hole, Notre-Dame cathedral fire). We also
consider 15.000 post titles of the subreddit TodayILearned (TIL dataset) and
13.000 post titles of the subreddit AskScience (AskScience dataset) on January
2019. We extracted the nouns, verbs, adjectives, and symbols from the textual
data. We run the experiments using the following parameters: $\alpha_{0}=0.5$,
$\theta_{0}=0.01$, $N_{samples}=2.000$, $N_{part}=8$ and
$\omega_{thres}=\frac{1}{2N_{part}}$. The kernel vector $\vec{\kappa}$ is
chosen as in [30]. It is made of Gaussian functions, with means located at
0.5, 1, 4, 8, 12, 24, 48, 72, 96, 120, 144 and 168 hours. The variance of each
is set to 1, 1, 3, 8, 12, 12, 24, 24, 24, 24, 24 and 24 hours. The shape of
the kernel is chosen so that it can account for a dynamic that can occur at
different timescales. The algorithm will then infer the weights $\vec{\alpha}$
associated with each entry of the kernel vector $\vec{\kappa}$ for each
cluster.

Figure 38: Wordclouds, triggering kernels and intensities for clusters the
most closely related to Sri Lanka 2019 bombings for various values of $r$. The
points at the bottom of the intensity plots are individual publication events.
Note that triggering kernels are plotted on a log-log scale for visualisation
purposes because most of the intensity is focused on small times: dynamics are
bursty.

######  PDHP recovers meaningful stories

As an illustrative example, we consider the inferred clusters the most closely
related to the Sri Lanka Easter bombings of April 2019 in Fig. 38. The main
bursts in the news related to this event happened on the 21st, 22nd and 23rd
of January, and respectively correspond to the bombings themselves, the
declaration of the state of emergency, and finally their application on the
23rd. We plot the temporal kernels associated with this event on a log-log
scale because most of the intensity is focused on small times: dynamics of
information spread are bursty [68]. We see that inferred dynamics change with
$r$ as well as the cluster’s vocabulary, which is expected since clusters do
not contain the same documents. For $r=0$, the Uniform process infers clusters
based on textual information only; the triggering kernel is inferred
afterwards. For $r=2.5$ on the contrary, clusters are formed based on the
triggering kernel, and textual information follows; we see from the right-plot
that this cluster captures publications exhibiting a daily intensity cycle;
this is visible both in the intensity plot (the bump around $2.10^{2}$h which
is not present on other temporal kernels) and in the real-time axis where
publications seem to be packed around specific points in time roughly
corresponding to a daily cycle. Given the intensity spikes on 21st, 22nd, and
23rd, it is not surprising that articles about the Sri Lanka bombings are also
part of this cluster. Note that the more $r$ increases, the more intense the
triggering kernel is around 24h. We see from Fig. 38 that DHP is a specific
case of our modelling, and that tweaking the $r$ parameter allows us to
retrieve completely different results.

######  PDHP favours temporal or textual clustering depending on r

We report the values of log-likelihoods for every dataset and various values
of $r$ in Fig. 39. The textual likelihood is defined in Eq.44, and the
likelihood of a Hawkes process is defined in Eq.42. Note that $r$ does not
appear in either Eq.44 or Eq.42; the plot in Fig. 39 thus only reflects the
relevancy of the proposed textual modelling or temporal modelling
independently from PDHP. Those likelihoods evaluate how well the textual or
temporal aspect of the dataset is modelled with no consideration of the model
being used. As expected from the synthetic experiments, varying $r$ makes the
model more sensitive to either textual or temporal data –note the similarity
to Fig. 35. A low $r$ favours the textual information clustering and is thus
better at modelling documents’ textual content, whereas a high $r$ favours
temporal information which makes PDHP better at capturing the publication
dynamics.

Figure 39: $r$ allows to favour text-based of time-based clustering on real
world datasets — Textual likelihood and Hawkes process likelihood for various
values of $r$. The lower $r$ the higher textual likelihood is, and the higher
$r$ the higher Hawkes process likelihood is.

######  PDHP infers sharper textual clusters for low r

We evaluate how meaningful textual clusters are using an entropy measure. We
assume that a cluster is meaningful when it contains a reduced set of words; a
cluster talking about one topic only is more likely to have a smaller
vocabulary than a cluster about two or more topics. A way to measure this is
to see how spread the vocabulary of a cluster is using Shannon entropy. Let
$N_{c,v}$ be the count of word $v$ in cluster $c$. The normalized Shannon
entropy of a cluster $c$ is defined as:

$S(\vec{N_{c}})=\frac{1}{-\log(V)}\sum_{v}^{V}\log(\frac{N_{c,v}}{\sum_{v}^{\prime}N_{c,v^{\prime}}})\frac{N_{c,v}}{\sum_{v}^{\prime}N_{c,v^{\prime}}}$
(70)

An entropy of 0 means the vocabulary of the cluster is concentrated on a
single word among the V possible words in the vocabulary; an entropy of 1
means that every of the $V$ words is present to the same extent with
probability $\frac{1}{V}$. In Fig. 40, we plot the mean entropy for various
values of $r$ for all the datasets, along with the standard error over the
clusters. The results show that on average vocabulary is more concentrated
within clusters for low values of $r$. The inflection point of the curves
corresponds to what has been previously observed with likelihoods in Fig. 39
and Fig. 36. On the contrary, higher values of $r$ lead to clusters that
comprise a less concentrated vocabulary. This is expected because as $r$
increases, the textual information is no longer the most relevant data for
cluster formation.

Figure 40: Textual clusters are more informative for low values of $r$ —
Weighted average entropy of words distribution for every dataset. Weights
corresponds to the number of words within clusters. The error bar represents
the standard error over all the clusters.

Figure 41: PDHP allows for modelling bursty events — PDHP’s intensity for the
News dataset for two values of $r$. A lower $r$ finds globally relevant
clusters, whereas a higher $r$ allows to recover shorter bursty events.

######  PDHP controls the burstiness

In Fig. 41, we plot the intensity function associated with the News dataset on
the real-time axis for several clusters for $r=0.5$ and $r=1.5$. Note that not
each of the $\sim$300 inferred clusters are represented, but instead we
consider only the ones whose intensity went above 10 at least once, the rest
being considered as noise. First of all, both values of $r$ allow to recover
the major events of April 2019 (in order of appearance): the first direct
picture of a black hole (10/04), the arrest of Julian Assange (11/04), the
fire of Notre-Dame de Paris Cathedral (15/04), the release of the Mueller
report on Donald Trump (18/04) and the Sri Lanka Easter bombings (21/04). The
top 5 words of every cluster are reported in the legend.

When $r$ increases, PDHP retrieves new clusters associated with shorter bursty
events. For instance, the cluster associated with the release of a new episode
of Equestria Girls that went unnoticed for $r=0.5$ has been detected with
$r=1.5$. This happens because the episode has not been discussed over a long
period, and associated articles have a vocabulary significantly overlapping
with other clusters’ one. A model relying mostly on textual information might
not detect specific words (twilight, Equestria, sparkle, etc.). If detected
and a new cluster is created, it might then fail to associate subsequent
events to this new cluster if temporal information plays a lesser role. On the
other side, a model favouring temporal information is much more likely to
associate subsequent events to a new cluster despite textual information
fitting well older and more populated clusters.

This can be seen in Fig. 38, where the intensity of a kernel peaks at short
times. This results in encouraging the burstiness. When $r$ is large, a given
event is likely to be associated with subsequent ones even when the associated
vocabularies are only vaguely similar. On the other side, when $r$ is small,
older events with closer vocabularies have more chances of getting associated
with it despite their intensity not peaking at the new event time.

######  Recovering publication cycles

The limit case of encouraging events burstiness is the deterministic
allocation of documents to a cluster based only on their relative positions on
the time axis. This is achieved when $r$ is large. In this case, textual
information does not matter and only regularities in the time distributions
are detected. We illustrate such a case on the News dataset in Fig. 42.

In Fig. 42, we plot on the left the intensity associated with the events for
each cluster on the real-time axis for $r=2$. We see that the two most
populated clusters follow precise dynamics. We added on the right side of the
plot the temporal kernel corresponding to each of these clusters. On the right
plot, we retrieve the cause of the daily and weekly cycles observed for the
largest cluster on the left plot. The second most populated cluster follows
similar dynamics, except that it seems to be shifted by half a day on the
real-time axis; the peaks are in phase opposition with the largest cluster. It
is worth noting that the Notre-Dame fire cluster is still detected; this is
due to its vocabulary being different enough from the existing cluster’s ones
to trigger its own cluster, and the associated number of documents being
consequent in a short time window. Interestingly, immediately after this
cluster emerged, the dynamics on the real-time axis also follow a decaying
circadian circle over three days.

Figure 42: Limit case of encouraging bursty events clustering — We plot PDHP’s
intensity for the News dataset over the observation period for a large value
of $r$. In this case, textual information plays a marginal role, and clusters
are inferred based on the events publication dates only.

######  Heuristics

Choosing $r$ — We saw that in all the previous experiments, the optimal $r$
got determined from a grid-search-like strategy. We did not come up with a way
to automatically infer the optimal $r$ without trying several values.

However, we provide some heuristics regarding the tuning of $r$.

  * •

As $r$ increases, we usually get a smaller number of inferred clusters. This
is because considering the temporal dimension adds consistency to the language
model; the temporal intensity prior for a new observation is likely to be non-
null, which increases the probability of not opening a new cluster with
respect to a model that does not consider time.

  * •

As $r$ goes to infinity, we only infer one large cluster that comprises all
the observations. This is because even the slightest difference in the prior
intensities leads to a deterministic cluster choice.

Our leads to automatically determine $r$ involve computing an ad-hoc objective
metric to optimize jointly with the likelihood. Given there is no gold
standard for clustering in real-world processes, the choice of $r$, and
therefore the choice of such metric, is left to the user. As we showed in Fig.
36 and later in Fig. 39, the choice of $r$ simply tunes the information on
which clusters are based. The clustering objective is to be defined for each
situation, which is made possible by manual tuning of $r$. Such an objective
could consist in minimizing the clusters’ word distribution entropy, the
standard deviation of the effective triggering kernel, or the average distance
between events within a cluster. A possible procedure for such optimization
would involve a multi-arms bandit to deal with this trade-off.

Number of clusters — In previous experiments, we compared our clustering
results to the ground truth using the NMI score. We chose this metric because
it allows us to compare a different number of clusters together. Indeed, it is
seldom the case that PDHP infers exactly the right number of clusters.

Typically, in our synthetic experiments with 2 ground-truth clusters, the
number of clusters could differ significantly at the beginning of the
algorithm (up to 10 clusters at once for small values of $r$). However, as the
number of observations increases, smaller clusters are discarded as the
algorithm converges toward the 2 correct clusters.

In real-world data, the number of clusters can grow very high –even more for
small values of $r$. However, the number of observations each of these
clusters comprise seems to follow a power-law distribution. Many of the
clusters contain very few observations (5 documents or less); they are
leftovers from the process as it converges towards more robust statistics.
This is why in Fig. 38 and Fig. 41, we restrict ourselves to the study of the
largest clusters only.

####  .18.4 Conclusion

In this section, we built the Powered Dirichlet-Hawkes process as a
generalization of the Dirichlet-Hawkes process and Uniform process. We
demonstrated how it improves performance on various datasets. When textual
information conveys little information, or when temporal information conveys
little information, and when both do, our model can correctly retrieve the
original clusters used in the generation process with high accuracy.

A central consideration in document clustering is that there are no “right”
clusters. For instance, we illustrate how textual content and temporal
dynamics can be decorrelated in real-life applications. PDHP is flexible
enough to allow to choose the weight they wish to give to temporal or textual
information depending on the situation; when textual and temporal clusters are
decorrelated, the model allows to choose which of those to infer.

Many future extensions are possible for PDHP. For instance, it would be
interesting to develop its hierarchical version (PHDHP) as it has already been
done with HDHP for DHP [83]. Given several recent works have been based on the
regular Dirichlet-Hawkes process, it would be insightful to study how their
results vary when using the Powered Hawkes-Dirichlet process instead. A study
of the influence of the language model used along with PDHP would also be
interesting since the text model we used here was simple on purpose (our focus
being on the PDHP prior and not on the model it gets associated with).
Typically, we would expect online Bayesian models that account for mutations
of textual clusters over time (e.g., shifts in vocabulary, mutating words,
etc.) to bring a significant improvement in modelling real-world systems [18,
5, 11, 140]

Regarding interaction modelling, we are now close to our objective. With PDHP,
we can:

  * •

Consider entities’ content. An entity is no more described as a unique
identifier, but instead by its semantic content. Two entities that convey the
same information are now considered as such and clustered together as a more
global entity (i.e., a topic here).

  * •

Model sparse interactions. Entities are now clustered together into temporal
clusters. It makes it feasible to spot interaction terms between sets of
entities. The lifespan of entities is no more a problem since clusters can
comprise entities spanning over extended periods, which also increases the
data available for each cluster.

  * •

Model dynamic interactions. Each cluster is associated with its own intensity
function, which determines its effect on ulterior observations. Eventually,
entities ’ influence fades away as time goes by.

In addition, we can tune the importance given to the temporal and textual
dimensions when modelling interactions. By varying $r$, a cluster can self-
stimulate essentially according to its entities publication times, or
essentially according to their content.

A major shortcoming of the proposed method is that interactions can only take
place within a cluster; they are self-interactions. As stated in the
introduction, this can be interpreted as the diagonal of the interaction
matrix in Fig. 12 in its temporal version. Ultimately, we are also interested
in studying how different clusters of entities influence each other.
Therefore, we propose to extend both [30] and the PDHP (Section .18) to the
multivariate case in the next section.

###  .19 Multivariate Powered Dirichlet-Hawkes Process – Final model

####  .19.1 Introduction

#####  .19.1.a Multivariate extension of PDHP

In this section, we extend the (univariate) Powered Dirichlet-Hawkes process
introduced in Section .18 to the multivariate case as the Multivariate Powered
Dirichlet-Hawkes Process (MPDHP). The various publications in a document’s
stream will now be able to influence each other. We detail and overcome
several technical challenges that arise from considering interacting topics,
and we conduct a systematic evaluation of MPDHP on a range of synthetic
datasets. As a first step, evaluation is conducted on synthetic datasets only,
so that we can discuss the performances and limitations of MPDHP in a
completely controlled environment. By the end of this section, we want to
determine whether it is possible to use MPDHP in a real-world setting

In previous works [18, 47, 30, 83], the understanding of large flows of data
boils down to summarizing them into a composition of independent groups
–clusters. However, as discussed throughout this manuscript, the processes at
stake are more complex than that, given these clusters are not independent of
each other; they interact.

Such interaction is illustrated in Fig. 43. This figure has to be considered
in regard to Fig. 29, where each group of documents was restricted to self-
interactions only; a given topic is assumed to only trigger observations from
the same topic. In Fig. 43, it would mean that the red cluster can only
trigger observations from the red cluster, as in Fig. 29. Here instead, we
consider that the red cluster also influences the probability of triggering an
observation from the blue cluster, and conversely.

#####  .19.1.b Workflow

In this section, we extend the models discussed in Section .16 and Section .18
to account for cluster interaction mechanisms. Firstly, we detail the
challenges that arise when developing the Multivariate Powered Dirichlet-
Hawkes Process (MPDHP). We show that alleviating them makes it possible to
achieve a linear time complexity $\mathcal{O}(N)$ (as in the original [30] and
Section .18) along with getting good clustering results. In doing so, we also
relax the near-critical Hawkes process hypothesis made in [30, 83], and
correct a major flow regarding the kernel choice in previous works. Secondly,
we conduct a systematic evaluation of the MPDHP in a variety of synthetic
situations. Our goal is to clearly identify the limits of MPDHP regarding
textual and temporal overlaps, computation time, the amount of available data,
the number of co-existing clusters, etc. This section is intended as a
technical report. By the end of it, a potential user should know in which
situations MPDHP can be useful, and in which ones alternative modelling
choices are required. .

Figure 43: Illustration of the Multivariate Powered Dirichlet-Hawkes Process
prior — A new event appears at time $t=3$ from a cluster which is yet to be
determined. The a priori probability that this event belongs to a given
cluster $c_{red}$ depends on the sum of the red dotted intensity functions at
time $t=3$. Similarly, the a priori probability that this event belongs to a
cluster $c_{blue}$ depends on the sum of the blue dotted lines at time $t=3$.
In previous existing models, this prior probability depends on each cluster
self-stimulation only.

####  .19.2 The Multivariate Powered Dirichlet-Hawkes process

#####  .19.2.a Multivariate Hawkes process

As discussed earlier, a Hawkes process is a temporal point process where the
appearance of new events is conditional on the realization of previous events.
It is fully characterized by an intensity function, noted $\lambda(t)$ that
depends on the history events generated by this process up to $t$,
$\mathcal{H}_{<t}$. We recall that the term $\mathcal{H}_{<t}$ is implicit
anytime the intensity function $\lambda$ is mentioned. A multivariate Hawkes
process is an extension of the Hawkes process, where the intensity function
$\lambda(t)$ depends on the history events generated by other Hawkes
processes. It means that in the definition of $\lambda(t)$, we cannot only
consider the events that happened in the same cluster, as in Eq. 39 and Eq.
66.

Similarly to Section .18, we associate one single Hawkes process to each
cluster. However, in (P)DHP, each of them is associated with a univariate
Hawkes process, which depends only on the history of events comprised in this
cluster. In our case, instead, we associate each cluster to a multivariate
Hawkes process that depends on all the observations previous to the time
being. Let $t_{i}^{c}$ be the time of realization of the $i^{th}$ event
belonging to cluster $c$. We write the intensity function for cluster $c$ at
all times as:

$\lambda_{c}(t)=\sum_{t_{i}^{c^{\prime}}<t}\vec{\alpha}^{T}_{c,c^{\prime}}\cdot\vec{\kappa}(t-t_{i}^{c^{\prime}})$
(71)

In Eq. 71, $\vec{\alpha}_{c,c^{\prime}}$ is a vector of $L$ parameters to
infer, and $\vec{\kappa}(t-t_{i}^{c^{\prime}})$ is a vector of $L$ temporal
kernel functions depending only on the time difference between two events. As
we will see later, $\alpha_{c,c^{\prime},l}$ is readily interpreted as the
influence of $c^{\prime}$ on $c$ according to the $l^{th}$ entry of the
temporal kernel. Once again, we consider a Gaussian RBF kernel, which allows
us to model a range of different intensity functions:

$\kappa_{l}(\Delta t)=\frac{1}{\sqrt{2\pi\sigma_{l}^{2}}}e^{-\frac{(\Delta
t-\mu_{l})^{2}}{2\sigma_{l}^{2}}}\ \ \ \forall l\in L$ (72)

The log-likelihood of a multivariate Hawkes process for all observations up to
a time $T$ is identical to the univariate case:

$\log\mathcal{L}(\alpha|\mathcal{H})=\sum_{c}\int_{0}^{T}\lambda_{c}(t)dt+\sum_{t_{i}^{c}}\lambda_{c}(t_{i}^{c})$
(73)

#####  .19.2.b Multivariate Powered Dirichlet-Hawkes Process

The Multivariate Powered Dirichlet-Hawkes Process (MPDHP) arises from the
merging of the Powered Dirichlet Process (Section .17) and the Multivariate
Hawkes Process (MHP), described in the previous paragraph. As in [30, 83] and
Section .18, the counts in PDP are substituted with the intensity functions of
a temporal point-process, here MHP. The a priori probability that a new event
is associated with a given cluster no longer depends on the population of this
cluster, but on its temporal intensity at the time the new observation
appears. This is illustrated in Fig. 43, where two events from two different
clusters $c_{red}$ and $c_{blue}$ have already happened at times $t_{0}=0$ and
$t_{1}=1$. A new event appears at time $t=3$. The a priori probability that
this event belongs to the cluster $c_{red}$ depends on the sum of the
intensity functions of observations at $t_{0}$ and $t_{1}$ on cluster
$c_{red}$ at time $t=3$ –sum of the red dotted lines. Similarly, the a priori
probability that this event belongs to the cluster $c_{blue}$ depends on the
sum of the blue dotted lines at time $t=3$.

Let $t_{i}$ be the time at which the $i^{th}$ event appears. The resulting
expression reads:

$P(C_{i}=c|t_{i},r,\lambda_{0},\mathcal{H})=\begin{cases}\frac{\lambda_{c}^{r}(t_{i})}{\lambda_{0}+\sum_{c^{\prime}}\lambda_{c^{\prime}}^{r}(t_{i})}\text{
if c$\leq$K}\\\
\frac{\lambda_{0}}{\lambda_{0}+\sum_{c^{\prime}}\lambda_{c^{\prime}}^{r}(t_{i})}\text{
if c=K+1}\end{cases}$ (74)

In Eq. 74, $\lambda_{c}$ in defined as in Eq. 71, and the parameter
$\lambda_{0}$ is the equivalent of the concentration parameter described in
Eq. 54. Taking back the illustration in Fig. 43, this parameter corresponds to
a time-independent intensity function. It has a chance to get chosen typically
when the other intensity functions are below it (meaning they do not manage to
explain the dynamic aspect of a new event). In this case, a new topic is
opened, and gets associated with its own intensity function.

#####  .19.2.c Language model

Similarly to what has been done in the previous sections and in [30], the
MPDHP must be associated with a Bayesian model given it is a prior on
sequential data. Since we study applications on textual data, we choose to
side the MPDHP prior with the same Dirichlet-Multinomial language model as in
previous sections. We recall the likelihood of the $i^{th}$ document belonging
to cluster $c$ reads (see Section. .16.2.d):

$\begin{split}\mathcal{L}(C_{i}=c|N_{<i,c},n_{i},\theta_{0})&=P(n_{i}|C_{i}=c,N_{<i,c},\theta_{0})\\\
&=\frac{\Gamma(N_{c}+\theta_{0})}{\Gamma(N_{c}+n_{i}+\theta_{0})}\prod_{v}\frac{\Gamma(N_{c,v}+n_{i,v}+\theta_{0,v})}{\Gamma(N_{c,v}+\theta_{0})}\end{split}$
(75)

where $N_{c}$ is the total number of words in cluster $c$ from observations
previous to $i$, $n_{i}$ is the total number of words in document $i$,
$N_{c,v}$ the count of word $v$ in cluster $c$, $n_{i,v}$ the count of word
$v$ in document $i$ and $\theta_{0}=\sum_{v}\theta_{0,v}$.

####  .19.3 Implementation

#####  .19.3.a Base algorithm

The Sequential Monte Carlo (SMC) algorithm used for the optimization has
already been described in Section .18.2.f and in Fig. 30. We briefly review it
in this section, before discussing the challenges that arise when using it in
the multivariate case.

######  SMC algorithm

The goal of the SMC algorithm is to jointly infer textual documents’ clusters
and the dynamics associated with them. It runs as follows. First, the
algorithm computes each cluster’s posterior probability for a new observation
by multiplying the temporal prior on cluster allocation (see Eq. 74,
illustrated Fig. 43) with the textual likelihood (see Eq. 75). It results in
an array of $K+1$ probabilities, where $K$ is the number of non-empty
clusters. A cluster label is then sampled from this probabilities vector. If
the empty $(K+1)^{th}$ cluster is chosen, the new observation is added to this
cluster, and its dynamics are randomly initialized (i.e., a $(K+1)^{th}$ row
and a $(K+1)^{th}$ column are added to the parameters matrix $\alpha$). If a
non-empty cluster is chosen, its dynamics are updated by maximizing the new
likelihood Eq. 73. The process then goes on to the next observation.

This routine is repeated $N_{part}$ times in parallel. Each parallel run is
referred to as a particle. Each particle keeps track of a series of cluster
allocation hypotheses. After an observation has been processed, we compute the
particles likelihood given their respective cluster allocations hypotheses.
Particles that have a likelihood relative to the other particles’ one below a
given threshold $\omega_{thres}$ are discarded and replaced by a more
plausible existing particle.

######  Sampling the temporal kernel

The parameters $\alpha$ are inferred using a sampling procedure. A number
$N_{sample}$ of pre-computed vectors are drawn from a Dirichlet distribution
with probability $P(\alpha|\alpha_{0})$, with $\alpha_{0}$ a concentration
parameter. As the SMC algorithm runs, within each existing cluster, each of
these candidate vectors is associated with a likelihood computed from Eq. 73,
noted $P(\mathcal{H}|\alpha)$, where $\mathcal{H}$ represents the data. The
sampling procedure returns the average of each of the $N_{sample}$ pre-
computed $\alpha$, weighted by the posterior distribution associated with them
${P(\alpha|\mathcal{H})\propto P(\mathcal{H}|\alpha)P(\alpha|\alpha_{0})}$.
The so-returned matrix is guaranteed to be a good statistical approximation of
the optimal matrix, provided the number of sample matrices $N_{sample}$ is
large enough.

######  Limits

This algorithm described here works well for the univariate case but fails for
the multivariate case. In particular, updating the multivariate intensity
function of each cluster requires knowing the number of already existing
clusters, which vary over time. Therefore, we cannot pre-compute the sample
matrices in advance –they must be updated as the algorithm runs to account for
the right number of non-empty clusters. Moreover, the number of parameters to
estimate evolves linearly with the number of active clusters $K$, instead of
remaining constant as in DHP and variants [30, 83]. Because the number of
parameters is not constant anymore, their candidate values cannot be sampled
from a Dirichlet distribution anymore. In the following, we review these
challenges and propose our solutions to overcome them. Eventually, we manage
to preserve a constant time complexity for each observation.

#####  .19.3.b Optimization challenges

######  Updating the triggering kernels

In the univariate case [30, 83] and Section .18, the coefficients
$\alpha_{c}\in\mathbb{R}^{L}$ are sampled from a collection of existing sample
vectors computed at the beginning of the algorithm (where $L$ is the size of
the kernel vector). However, we must now infer a matrix instead. We recall
that matrix $\alpha_{c}$ represents the weights given to the temporal kernel
vector of every cluster influence on $c$ –see Eq. 71. The likelihood Eq. 73
can be updated incrementally for each sample matrix. A given cluster $c$ has a
likelihood value associated with each of those $N_{sample}$ sample matrices,
which represents how fit one sample matrix is to explain one cluster’s
dynamics. The final value of the parameters matrix is then computed by
sampling, simply averaging the sample matrices weighted by their likelihood
for a given cluster times the prior probability of these vectors being drawn
in the first place.

However, such sampling was possible in the univariate case, where each sample
matrix was in fact a vector of fixed length. In our case, because Hawkes
processes are multivariate, each entry $\alpha_{c}\in\mathbb{R}^{K\times L}$
is now a matrix. Moreover, the number of existing clusters $K$ increases over
time and can grow large. Each time a new cluster is added to the computation,
a row is appended to the $\alpha_{c}$ matrix –it accounts for the influence of
this new cluster regarding $c$.

Our solution is that some events can be discarded from the computation so that
some old clusters can also be discarded. Clusters whose last observation has
exceeded a certain age has a near-zero chance to get sampled once again. It
means these clusters’ contribution to the likelihood Eq. 73 is fixed.
Therefore, they do not have a role in the computation of the parameters matrix
$\alpha_{c}$. The row corresponding to each of these clusters in the
parameters matrix can then be discarded in every sample matrix. Put
differently, the last sampled value for their influence on $c$ will remain
unchanged for the remaining of the algorithm. The dimension of $\alpha_{c}$
only depends on the number of active clusters, whose intensity function has
not faded to zero. For a given dataset, the number of active clusters
typically fluctuates around a constant value, making one iteration running in
constant time $\mathcal{O}(1)$.

######  A beta prior on parameters

Another problem inherent to the proposed multivariate modelling is the prior
assumption on sample vectors. In [30, 83], each sample vector is sampled from
a Dirichlet distribution. This choice is to infer Hawkes processes that are
nearly unstable: the spectral radius of their temporal kernel function
$\lambda_{c}(t)$ is close to 1. However in our case, such an assumption is not
possible because the size of each sample matrix can vary as the number of
active clusters evolve. Drawing one Dirichlet vector of size $L$ for each
entry $\alpha_{c,c^{\prime}}$ would force the spectral radius of $\alpha_{c}$
to equal $K$, which transcribes a highly unstable Hawkes process. Our solution
is to consider the parameters as completely independent of each other. Each
entry of the matrix $\alpha$ is drawn from an independent $\beta$ distribution
of parameter $\beta_{0}$. In this way, we make no assumption on the spectral
radius of the Hawkes process, and sample rows/columns corresponding to new
clusters can be generated one after the other.

######  On the temporal concentration parameter $\lambda_{0}$

While not specifically related to the implementation of the multivariate case,
we discuss in this paragraph an important consideration when designing DHP-
based models. In most recently published works on the topic [30, 83] and in
the experiments conducted in Section .18, inference on real-world processes is
done using an RBF temporal kernel. It means that time is paved with Gaussian
functions centered at various points in time; the parameter $\alpha$ in DHP-
based models accounts for the weights given to each of these Gaussian
functions.

In these works, the kernel is chosen so that it accounts for different time
scales by centering Gaussian functions on unevenly spaced points in time. The
standard deviation of each of these entries varies to account for larger time
ranges. However, all values of a Gaussian function are small when their
standard deviation is large, for normalization reasons –the maximum value of a
Gaussian function whose standard deviation is $\sigma$ is
$\frac{1}{\sqrt{2\pi\sigma^{2}}}$.

In the SMC algorithm, this RBF kernel is evaluated at a single point in time
and confronted with the temporal concentration parameter $\lambda_{0}$ (see
Eq.74) to determine whether to open a new cluster. In [30], such values are
compared to $\lambda_{0}$ constant in time. It means that, mechanically, these
methods cannot detect observations triggered by such Gaussian functions as
their value is systematically lower than $\lambda_{0}$ –typically at long time
ranges in [30, 83], which can be seen from these articles’ kernel plots that
fade as time goes. We explicitly illustrate how this leads to a bias in the
modelling in Fig. 44.

Figure 44: Choosing the right temporal concentration parameter $\lambda_{0}$ —
The choice of the temporal concentration parameter $\lambda_{0}$ can lead to
bias. (Left) The problem with its choice in [30] is that events happening at
large time ranges are likely to go undetected, as the Hawkes intensity at
these ranges cannot be larger than $\lambda_{0}$. (Middle) A first solution
consists in paving the space with evenly spaced Gaussian functions that all
share the same standard deviation. (Right) A second solution is to make
$\lambda_{0}$ a function of time so that its ratio with the temporal kernel
remains constant.

Consider for instance the RBF kernel used in [30, 83], plotted in Fig. 44,
with the Gaussian means equal to 0.5, 1, 8, 12, 24, 48, 72, 96, 120, 144 and
168 hours, and the corresponding deviations equal to 1, 1, 8, 12, 12, 24, 24,
24, 24, 24, and 24 hours. The authors used $\lambda_{0}=0.01$. For the last
entry of their RBF kernel, the maximum value of the Gaussian function
$\mathcal{G}(\mu=168;\sigma=24)$ is about $3.10^{-4}$, which is much smaller
than $\lambda_{0}=1.10^{-2}$. It means that even for a cluster whose intensity
function only acts at long-ranges, the chances of spotting events triggered by
such clusters are about 3%. This makes the models presented in [30, 83] unfit
to spot long-range interactions.

There are two ways to overcome this problem (Fig. 44-Left) so that
$\lambda_{0}$ can be consistently confronted to the clusters’ temporal kernels
(see Eq. 71):

  * •

Fig. 44-Middle – To consider an RBF kernel whose Gaussian functions all share
the same deviation, while keeping $\lambda_{0}$ constant. We choose this
solution in the follow-up experimental section.

  * •

Fig. 44-Right – To consider a $\lambda_{0}$ that can vary in time according to
the maximum value of the RBF kernel at different time points –which depends on
their standard deviation.

####  .19.4 Experiments

#####  .19.4.a Setup

We now design a series of experiments to explore the possible use domain for
the Multivariate Powered Dirichlet-Hawkes Process. We list the parameters we
consider in our experiments. When a parameter does not explicitly vary, it
takes a default value given between parentheses. These parameters are: the
textual overlap (0) and the temporal overlap (0) discussed further in the
text, the temporal concentration parameter $\lambda_{0}$ (0.01), the strength
of temporal dependence $r$ (1), the number of synthetically generated clusters
$K$ (2), the number of words associated with each document $n_{words}$ (20),
the number of particles $N_{part}$ (10) and the number of sample matrices used
for sampling $\alpha$, noted $N_{sample}$ (2,000). For the detail of these
parameters, please refer to Eq. 74.

Note that the overlap $o(f_{1},...,f_{N})$ between $N$ functions is defined as
the sum over each function $f_{i}$ of its intersecting area with the largest
of the $N-1$ other functions, divided by the sum of each function’s total
area. This value is bounded between 0 (perfectly separated functions) and 1
(identical functions). Mathematically:

$o(f_{1},...,f_{N})=\frac{1}{\sum_{i}\int_{\mathbb{R}}f_{i}(x)dx}\sum_{i}\int_{\mathbb{R}}min(f_{i}(x),max(\\{f_{j}(x)\\}_{j\neq
i}))dx$ (76)

For each combination of parameters considered, we generate 10 different
datasets. In all datasets, we consider a fixed size vocabulary $V=1000$ for
each cluster. All datasets are made of 5,000 observations. Observations for
each cluster $c$ are generated using an RBF temporal kernel $\vec{\kappa}(t)$
weighted by a parameter matrix $\alpha_{c}$. Explicitly, we set
$\vec{\kappa}(t)=\left[\mathcal{G}(3;0.5);\mathcal{G}(7;0.5);\mathcal{G}(11;0.5)\right]$
where $\mathcal{G}(\mu;\sigma)$ is a Gaussian function of mean $\mu$ and
deviation $\sigma$ –following the discussion raised in Section .19.3.b.3. We
note $L=3$ the number of entries of $\vec{\kappa}$. The inferred entries of
$\alpha$ determine the amplitude (or weight) of each Gaussian kernel function.

The generation process is as follows. First, we draw a random matrix
$\alpha\in\mathbb{R}^{K\times L}$ and normalize it so that its spectral radius
equals 1 –near unstable Hawkes process. We repeat this process until we obtain
the wanted temporal overlap.444Note that the overlap as defined here is
different from the one used in Section. .18. In the latter, we considered the
overlap of the intensity function plot on the real-time axis. Here instead we
consider the overlap between the kernel intensity functions. This is because
the temporal overlap as defined in Section .18 is always close to 1 in the
multivariate case, because different clusters’ events strongly interact with
each other. Then, we simulate the multivariate Hawkes process using the
triggering kernels $\vec{\alpha}\cdot\vec{\kappa}(t)$, where $\vec{\kappa}(t)$
is the RBF kernel as defined earlier. Given the Hawkes process is
multivariate, each event is associated with its class it has been generated
from among $K$ possible classes. For each so generated event, we draw
$n_{words}$ words from a vocabulary of size $V$. The vocabularies are drawn
from a multinomial distribution and shifted over this distribution so that
they overlap to a given extent (see Eq.76).

#####  .19.4.b Baselines

We compare our approach to 3 closely related baselines.

  * •

Dirichlet-Hawkes process (DHP) [30] – In this model, clusters can only self-
replicate. It means that the intensity function of a cluster $c$ Eq. 71 only
considers past events that happened in the same cluster $c$.

  * •

Dirichlet process (DP) – This prior is standard in clustering problems. It
corresponds to a special case of Eq. 54 where $r=1$. It assumes that the prior
probability for an observation to belong to a cluster depends linearly on the
population of this cluster.

  * •

Uniform process (UP) [124] – This prior corresponds to a special case of Eq.
54 where $r=0$. It assumes that the prior probability for an observation to
belong to a cluster does not depend on any information about this cluster
(neither population nor dynamics).

As in previous experiments, we evaluate our results in terms of normalized
mutual information (NMI) score. We recall that this metric is standard when
evaluating non-parametric clustering models. It compares two cluster
partitions (i.e., the inferred and the ground truth ones); it is bounded
between 0 (each true cluster is represented to the same extent in each of the
inferred ones) and 1 (each inferred partition comprises 100% of one true
cluster).

#####  .19.4.c Results

Figure 45: Experimental results on synthetic data — MPDHP consistently
outperforms other baselines by considering both textual information and
temporal information.

######  MPDHP outperforms state-of-the-art

In Fig. 45, we plot our results for datasets that has been generated using a
Multivariate Hawkes process (clusters have an influence on each other) and a
Univariate Hawkes process (clusters can only influence themselves). We compare
MPDHP to the proposed baselines for various values of textual overlap. We draw
the following conclusions:

  * •

MPDHP systematically outperforms the proposed baselines on multivariate data
–when clusters interact with each other. Considering that clusters interact
with each other improves our description of the datasets.

  * •

MPDHP performs at least as good as PDHP on univariate data –when clusters can
only self-stimulate. The complexity induced by MPDHP does not make it unfit
for simpler tasks.

  * •

PDP performs well for small textual overlaps, but rapidly fails when the
textual overlap increases. This is expected since only the textual information
is considered by the PDP. It highlights the importance of also considering the
temporal information.

  * •

PDHP performs better than MPDHP when the textual overlap is large (textual
overlap of 0.8). This is due to the increased complexity of MPDHP over PDHP.
In challenging situations such as this one, a simpler model takes fewer
efforts to overcome initialization mistakes, as there are fewer parameters to
put back on the right track.

  * •

All models fail when the textual overlap is complete; clusters cannot be
inferred from temporal information only.

######  Uninformative textual content and entangled dynamics

In Fig. 46, we plot the results of MPDHP for different values of textual and
temporal overlap. Textual overlap is defined as in Eq. 76. According to Eq.
71, the influence kernel of cluster $c^{\prime}$ on cluster $c$ can be written
$\vec{\alpha}_{c,c^{\prime}}\cdot\vec{\kappa}(t)$. For each cluster $c$, we
generate values of $\alpha_{c,c^{\prime}}$ so that the overlap between all the
functions in the set
$\\{\vec{\alpha}_{c,c^{\prime}}\cdot\vec{\kappa}(t)\\}_{c^{\prime}}$ equals a
given value. The idea is to evaluate whether MPDHP is robust when clusters
have similar dynamics.

Overall, we see that when the textual overlap is small, MPDHP yields good
results independently from the temporal overlap. It means that in this case,
the textual content is enough to differentiate clusters despite their dynamics
being similar. However, as textual content gets less informative (textual
overlap $\geq$ 0.6), results are better when the temporal overlap is low. In
these cases, textual information is not enough and MPDHP relies more on
temporal data. Overall, MPDHP handles challenging cases provided either
textual or temporal information is informative enough – for instance temporal
overlap of 0 and textual overlap of 0.7, or temporal overlap of 1 and textual
overlap of 0.4. It fails when both are uninformative – for instance, temporal
and textual overlaps of 1.

Figure 46: MPDHP handles scarce textual or temporal information — MPDHP
handles challenging cases provided either textual of temporal information is
informative enough (temporal overlap of 0 and textual overlap of 0.7; temporal
overlap of 1 and textual overlap of 0.4) and fails when both are uninformative
(overlaps of 1).

######  Highly interacting processes

Next, we assess whether MPDHP works when a large number of clusters coexist
simultaneously. The rate at which new clusters get opened is mainly controlled
by the $\lambda_{0}$ hyperparameter (see Eq. 74), which we vary to see whether
MPDHP is robust against it. In Fig. 47, we plot the performances of MPDHP
according to these two parameters. We can draw two conclusions:

  * •

MPDHP can handle a large number of coexisting clusters and still correctly
identify to which one each document belongs.

  * •

MPDHP is robust versus variations of $\lambda_{0}$. In this case, results are
similar for $\lambda_{0}$ varying over 5 orders of magnitude. It means MPDHP
does not have to be fine-tuned according to the number of expected clusters in
cases where this number is not known in advance.

Figure 47: MPDHP can handle a large number of coexisting clusters — MPDHP
yields good results when a large number of clusters coexist simultaneously. It
is also robust against variations of $\lambda_{0}$ over 5 orders of magnitude.

######  Handling scarce textual information

In this paragraph, we determine how much data should be provided to MPDHP to
get satisfying results. In Fig. 48, we plot the performances of MPDHP with
respect to the number of words generated by each observation and to the
clusters’ vocabulary overlap. MPDHP needs a fairly small number of words to
yield good results over 5,000 observations. For reference, the overlap between
topics can be estimated at around 0.25 ([95], in Spanish). Similarly, we can
estimate an average of $\sim$10-20 named entities per Twitter post (240
characters). These results support the application of MPDHP to model real-
world situations.

Figure 48: How much data to use MPDHP — MPDHP performances according to the
number of words associated with each event and the overlap between clusters’
vocabulary. Overall, MPDHP needs little data. This plot provides a map of how
much data must be provided to MPDHP to make it work. For reference, topic
overlaps in the real world can be estimated at around 0.25 and the number of
named entities per Twitter post around 20.

######  Computational needs

Finally, we investigate how much computational resources should be allocated
to MPDHP’s sequential Monte-Carlo (SMC) inference algorithm to obtain good
results. In Fig. 49, we plot the model’s performance against the two main
optimization parameters –the number of sample matrices and the number of
particles. We recall that the sample matrices are used to infer each value of
the kernel weights matrices for cluster $c$, noted $\alpha_{c}$; the more
sample matrices, the better the estimation. The number of particles represents
the number of different cluster allocations hypotheses explored by the SMC
algorithm at each step; the more particles, the more hypotheses are tested
simultaneously. Overall, we see that MPDHP works well with few resources. In
our experiments, results do not seem to improve significantly when using more
than 20 particles, and when using more than 1000 sample vectors.

Figure 49: How complex should the algorithm be — Performance of MPDHP using
different versions of the Sequential Monte-Carlo algorithm. Here, we plot the
model’s performance with respect to the number of sample matrices used to
estimate the kernel’s weights $\alpha_{c}$ (see Eq. 71) and the number of
particles $N_{part}$ used for the inference. Overall, MPDHP functions well
with few computational resources.

####  .19.5 Conclusion

In this section, we extended existing priors (the Dirichlet-Hawkes process and
the Powered Dirichlet-Hawkes process) so that they can consider multivariate
Hawkes processes, resulting in the resulting Multivariate Powered Dirichlet-
Hawkes process (MPDHP). This new process is used as a Bayesian prior coupled
to a textual model to infer clusters temporal interaction network from textual
data flow. Along with its derivation came several optimization challenges,
that we overcome to preserve a computational time that scales linearly with
the size of the dataset $\mathcal{O}(N)$.

Through systematic experiments, we tested our approach against state-of-the-
art models and explored its limits by varying the parameters used for
synthetic data generation. We show that MPDHP outperforms existing baselines
when clusters are allowed to interact with each other, and performs at least
as well as the PDHP baseline when clusters are only allowed to self-interact
(which PDHP is designed to model). We show that MPDHP can handle cases where
textual content is uninformative better than other baselines. Besides, it
handles cases where temporal dynamics are similar across clusters. We also
show that MPDHP is robust against tuning of the temporal concentration
parameter $\lambda_{0}$, which allows it to handle highly intricate processes
where a lot of clusters coexist simultaneously. We evaluated the robustness of
MPDHP against the number of words observed for each event and the overlap
between various clusters and showed that it performs well with few textual
data when vocabularies overlap is not total. Finally, we discussed the
computational needs of PDHP, and show that it works correctly with minimal
computational resources.

The present section was intended as a report on what can and what cannot be
achieved using MPDHP. Our various results suggest that this prior can be
applied in a robust way to a broad range of problems for a minimal
computational cost. In particular, the results from these extensive
experiments support the possibility of applying MPDHP to real-world
situations.

Regarding the task at hand, it appears that MPDHP provides a robust way to
model interactions in information spread. The PDHP allowed to model sparse and
dynamic self-interactions between semantic entities; these interactions can
now take place between different entity clusters.

###  .20 Case study on a real-world dataset – Reddit news

####  .20.1 Introduction

Throughout the previous sections, we developed a plethora of models to
investigate interactions in information spread. A first approach underlined
the necessity of considering clustering and a second one the necessity of
considering the temporal dimension. To answer these conclusions, we extended a
promising class of models, the Dirichlet-Hawkes process, so that it becomes
possible to spot temporal interactions between clusters of textual documents
in large-scale settings. The resulting Multivariate Powered Dirichlet-Hawkes
Process (MPDHP) answers all these constraints and yields interpretable results
on interacting processes.

As a closure of our work on interactions modelling in information spread, we
propose to conduct a large-scale experiment on a real-world dataset using
MPDHP. In this section, we first describe and justify the choice of a large-
scale real-world dataset from Reddit. In a second time, we conduct an in-depth
analysis of the results yielded by MPDHP. We finally conclude on the role of
interactions in the spread of news headlines on Reddit.

####  .20.2 Dataset

#####  .20.2.a Origin and raw data

######  Why Reddit?

Reddit is a social network platform that gathers 48M of monthly users, as of
today. Reddit is composed of a galaxy of forums (or subreddits) that are often
specific to a topic or an organisation. Within each forum, users can open a
new post made of a title and a content (typically textual or visual). Users
can then comment on the post and get answered to. The contents published on
the platform are therefore user-generated and highly dynamic, with an
estimated 1M of new publications per day. Each post can be upvoted (score +1)
or downvoted (score -1); the popularity of the post is defined as the
difference between upvotes and downvotes.

For these reasons, this corpus fits our demonstration: information emerges
from a large user-generated data flow, its contents are formatted, textual,
timestamped and user-generated, and the number of daily publications makes it
likely that some of them interact with each other.

######  Data

We collected our dataset from the Pushshift Reddit repository [12]. As its
authors describe it, “Pushshift is a social media data collection, analysis,
and archiving platform that since 2015 has collected Reddit data and made it
available to researchers. Pushshift’s Reddit dataset is updated in real-time
and includes historical data back to Reddit’s inception. […] The Pushshift
Reddit dataset makes it possible for social media researchers to reduce time
spent in the data collection, cleaning, and storage phases of their projects.”

In practice, the dataset comprises the entirety of the content posted on
Reddit up to June 2021. In particular, for each Reddit post, we can retrieve
the subreddit it came from, the title of the publication, its publication date
and its score.

#####  .20.2.b Preprocessing

######  Selecting the news subreddits

For the need of our study, we restrict ourselves to consider only popular
English news subreddits. Namely, we select only posts from the following
subreddits: inthenews, neutralnews, news, nottheonion, offbeat, open news,
qualitynews, truenews, worldnews. This first routine leaves us with 867,328
headlines, which makes a total of 1,111,955 words drawn from a vocabulary of
size 36,284.

######  Cleaning the textual data

As it is common in natural language processing, we must clean the raw text
extracted from the Reddit posts so it becomes usable. To do so, we conduct the
following routine:

  1. 1.

Remove the web addresses

  2. 2.

Put the text in lowercase

  3. 3.

Remove punctuation signs

  4. 4.

Remove extra white spaces

  5. 5.

Remove all English stopwords (imported from nltk)

  6. 6.

Remove all words whose length is lesser than 4

  7. 7.

Remove all words that appear less than 3 times in the original dataset

######  Removing uninformative documents

Next, we remove publications that carry lesser textual or temporal
information.

Firstly, we choose not to consider the publications that have a popularity
lesser than 20 – meaning that they received less than 20 positive votes more
than negative votes. We make this choice so that we consider publications that
are visible enough to have any influence on the data generation process.

Secondly, we remove the publications that comprise less than 3 words. The
semantic information carried is expected to be poor and is not considered in
our analysis.

######  Final dataset

After curating the dataset in the way described above, we are left with
102,045 news headlines (one-eighth of the original data), which makes a total
of 875,334 tokens (named entities, verbs, numbers, etc.) drawn from a
vocabulary of size 13,241 (one-third of the original vocabulary). The
characteristics of this dataset are shown in Fig. 50.

Figure 50: Characteristics of the News dataset — For $\sim$100,000 headlines
and $\sim$13,000 different words (Top-Left) Distribution of the words count
(Top-Right) Distribution of headlines popularity (Bottom-Left) How many
headlines per subreddit (Bottom-Right) How publications spread over time

####  .20.3 Experimental setup

As announced earlier, we will apply the MPDHP of Section .19 to the News
dataset detailed above. First, we must determine which hyper-parameters to
use.

######  Temporal kernel

We run our experiments using three different RBF kernels, which should account
for publication dynamics at three different timescales: minute, hour, and day.
The “Minute” RBF kernel is made of Gaussian functions centered at the
following times: $\left[0,10,20,30,40,05,60,70,80\right]$ minutes; each entry
shares a same standard deviation $\sigma$ of 5 minutes, and
$\lambda_{0}=0.01$. The “Hour” RBF kernel has Gaussians centered around
$\left[0,2,4,6,8\right]$ hours, with a standard deviation $\sigma$ of 1 hour,
and $\lambda_{0}=0.001$. The “Day” RBF kernel is centered around
$\left[0,1,2,3,4,5,6\right]$ days, with a standard deviation $\sigma$ of 0.5
days, and $\lambda_{0}=0.0001$.

For each of these kernels, we set the concentration parameter $\lambda_{0}$ so
that it reaches roughly the value of one Gaussian function evaluated at
$2\sigma$. It means that an event which is $2\sigma$ away from the center of
the Gaussian kernel of a single observation has 50% chances of getting
associated with this Gaussian kernel entry, and 50% chances of opening a new
cluster. It also translates that the chances to open a new cluster given an
observation right in the center of a single RBF kernel entry are roughly
$\frac{1}{8}$.

This allows us to spot interactions that might occur at different timescales.
However, one of these timescales is likely to explain the data much better
than the others. If that is the case, it would mean that the dynamics at stake
are not scale-free –that interactions do not have noteworthy influence at
every temporal scale.

######  Language model

We use the same Dirichlet-Multinomial language model as in [30] and in Section
.18 and Section .19. We expect the vocabulary to show enough variations across
the range of topics discussed in the News dataset. We vary the concentration
parameter of this model between $\theta_{0}=0.001$ and $\theta_{0}=0.01$. The
choice of these values is supported by other works using a similar model [30]
and by our own observations in Section .18. A larger value of $\theta_{0}$
makes the inferred clusters cover a broader range of document types, whereas a
small value makes the inferred clusters more specific to a topic.

######  SMC algorithm

Finally, we detail and justify the parameters used for running the SMC
algorithm. The procedure to update the optimal parameters matrix $\alpha$
relies on the average of a set of sample vectors weighted by these vectors’
likelihood. The number of active clusters can grow large, and the number of
parameters scales as the square of this number. We therefore need to have
enough sample matrices to guarantee a good approximation of the optimal
$\alpha$ –the value recommended for 2 clusters in Fig. 49 is likely to be too
small for the task.

We set $N_{samples}=100000$. From our observations, the number of coexisting
clusters can go as large as 80 coexisting ones (roughly 40,000 parameters to
estimate), that is 2.5 sample values per parameter. In practice however, the
number of coexisting clusters remains fairly low, around 10 coexisting
clusters (roughly 1,000 parameters), allowing sampling each parameter from 100
candidate values.

Each sample parameter is drawn from an identical Beta distribution of
concentration parameter $\alpha_{0}=2$. We set this value so that extreme
values (0 and 1) for these parameters are rare, and so that the sampled
parameters matrices can show great variations from one another.

Finally, we consider 8 particles for the SMC algorithm, similarly to what is
done in [30] and in Section .18, and recommended in Fig. 49.

####  .20.4 Results

#####  .20.4.a Overview of the experiments

In our experiments, we used three temporal kernels that account for dynamics
at different time scales, and tested diverse values for the language model
concentration parameter $\theta_{0}$ and for the exponent $r$. In Table 11, we
represent the main characteristics for each of these runs in terms of number
of inferred clusters $K$, the average cluster population $<N>$ (where
$<\cdot>$ denotes the average), the average normalized entropy of the
vocabulary of the top 20 clusters $S_{text}^{(20)}$, the average normalized
entropy of the subreddits partition of the top 20 clusters $S_{sub}^{(20)}$.

We recall that the normalized entropy is defined as:

$S(\vec{x})=-\frac{1}{\ln|\vec{x}|}\sum_{i}^{|\vec{x}|}x_{i}\ln x_{i}$ (77)

where $\vec{x}$ is a vector that sums to $1$ and $|\vec{x}|$ its cardinal
(length). Each entry $x_{i}$ represents the probability of $i$. When
considering counts, $\vec{x}$ can be set equal to the frequency of each
observation. The entropy is normalized between 0 (minimal spread,
$\vec{x}_{i}=\delta_{ij}\,\forall i$) and 1 (maximal spread,
$\vec{x}_{i}=\frac{1}{|\vec{x}|}\,\forall i$). In our case, a low entropy
$S_{text}^{(20)}$ (resp. $S_{sub}^{(20)}$) means that clusters contain
documents that are concentrated around a reduced set of words (resp. of
subreddits); conversely, a large entropy means that clusters do not account
for documents concentrated around a specific vocabulary (resp. set of
subreddits).

$\vec{\kappa}(t)$ | $\theta_{0}$ | $r$ | $K$ | $<N>$ | $S_{text}^{(20)}$ | $S_{sub}^{(20)}$
---|---|---|---|---|---|---
Minute |  0.01 | 0.0 | 16150 | 6 | 0.744(68) | 0.400(36)
0.5 | 8498 | 12 | 0.796(55) | 0.441(24)
1.0 | 5069 | 20 | 0.790(49) | 0.476(53)
1.5 | 2730 | 37 | 0.808(43) | 0.490(52)
2-7 |  0.001 | 0.0 | 46304 | 2 | 0.475(48) | 0.239(56)
| 0.5 | 37277 | 3 | 0.485(40) | 0.256(60)
| 1.0 | 26858 | 4 | 0.501(37) | 0.266(72)
| 1.5 | 19275 | 5 | 0.506(39) | 0.280(58)
Hour |  0.01 | 0.0 | 3792 | 27 | 0.798(52) | 0.474(106)
0.5 | 1735 | 59 | 0.791(45) | 0.469(77)
1.0 | 825 | 124 | 0.803(47) | 0.484(75)
1.5 | 426 | 240 | 0.795(37) | 0.489(69)
2-7 |  0.001 | 0.0 | 18012 | 6 | 0.760(86) | 0.397(63)
| 0.5 | 11923 | 9 | 0.784(72) | 0.425(38)
| 1.0 | 4837 | 21 | 0.821(50) | 0.497(30)
| 1.5 | 2368 | 43 | 0.814(42) | 0.481(91)
Day |  0.01 | 0.0 | 609 | 168 | 0.713(34) | 0.413(103)
0.5 | 326 | 313 | 0.728(33) | 0.429(98)
1.0 | 172 | 593 | 0.743(31) | 0.461(94)
1.5 | 96 | 1063 | 0.755(36) | 0.464(86)
2-7 |  0.001 | 0.0 | 4349 | 23 | 0.705(49) | 0.396(103)
| 0.5 | 2654 | 38 | 0.721(59) | 0.404(104)
| 1.0 | 1399 | 73 | 0.734(57) | 0.431(106)
| 1.5 | 764 | 134 | 0.741(54) | 0.442(98)
Table 11: Results for each experiment — We tried various combinations of
parameters $\vec{\kappa}(t)$, $\theta_{0}$ and $r$ and observe how they result
in a variety of outputs. We characterize these outputs in terms of clusters
(number, size, textual entropy, subreddits entropy). The standard deviation on
the last digits is given in standard notation – $0.123(12)\Leftrightarrow
0.123\pm 0.012$.

We can make several observations from Table 11:

  * •

We recover the fact that textual clusters have a lower entropy for small
values of $r$; this is because their creation is based more on textual
coherence than on temporal coherence.

  * •

The subreddit entropy seems to grow with $r$, but no clear trend is visible
due to large error bars. A possible interpretation is that favouring the
temporal information for cluster creation results in larger clusters. They
would be too large to account for subreddit-specific dynamics. However, the
entropy remains fairly low. The entropy of the distribution Fig. 50-top-left
is equal to 0.51.

  * •

The number of inferred clusters decreases with $r$, and their average
population increases.

  * •

The number of clusters grows large for the “Minute” kernel. This is because
the short time range considered does not allow for clusters to last in time. A
cluster that does not replicate within 1h30 is forgotten.

######  Choosing a timescale

From Table 11, we see there can be a large variety of outputs to analyse,
depending on the modelling choices. If we are interested in the micro-
interactions that happen at tiny time scales, the “Minute” RBF kernel should
be considered. However, the short time range it allows for interactions makes
the number of clusters large, and the average cluster population small. On the
other hand, the “Day” RBF kernel spans over larger periods, which prevents
discarding clusters too soon. For instance, it will not discard clusters that
follow a circadian publication dynamic, unlike the “Minute” kernel that cannot
account for such time ranges.

######  Choosing $\theta_{0}$

In a Dirichlet-Multinomial textual model, such as written in Eq. 44, the
hyperparameter $\theta_{0}$ controls the concentration of topics’ vocabulary.
A small value of $\theta_{0}$ makes it so that a new document should have an
almost identical word distribution as a given cluster to enter it; there are
more chances that small discrepancies lead to opening a new cluster.
Conversely, larger values of $\theta_{0}$ allow documents to belong to
clusters even if their word distribution does not fit exactly the cluster’s
content. We tested a small value $\theta_{0}=0.001$ and a large value
$\theta_{0}=0.01$, which are standard in usual text modelling [19, 18, 30].