Nonnegative Bayesian nonparametric factor models with completely random measures for community detection
Fadhel Ayed François Caron






Bayesian Nonparametric Poisson factorization


  
  

  Department of Statistics, University of Oxford

           e1,e2



  University of Oxford


  F. Ayed and F. Caron




We present a Bayesian nonparametric Poisson factorization model for modeling network data with an unknown and potentially growing number of overlapping communities. The construction is based on completely random measures and allows the number of communities to either increase with the number of nodes at a specified logarithmic or polynomial rate, or be bounded. We develop asymptotics for the number of nodes and the degree distribution of the network and derive a Markov chain Monte Carlo algorithm for targeting the exact posterior distribution for this model. The usefulness of the approach is illustrated on various real networks.





Poisson factorization
Community detection
Generalized gamma process












§ INTRODUCTION



Non-negative matrix factorization (NMF) methods <cit.> aim to find a latent representation of a positive n× m matrix A as a sum of K non-negative factors. For integer-valued data, Poisson factorization models <cit.> offer a flexible probabilistic framework for non-negative matrix factorization, and have found wide applicability in signal processing <cit.> or recommender systems <cit.>. In this paper, we focus on the application to network analysis, where m=n and the n× n count matrix A represents the number of directed or undirected interactions between n individuals; the latent factors may be interpreted as latent and potentially overlapping communities <cit.>, such as sport team members or other social activities circles. We also consider binary data where the matrix represents the existence or absence of a directed or undirected link between individuals. The estimated latent factors can be used for the prediction of missing links/interactions, or for interpretation of the uncovered latent community structure.


Poisson factorization approaches require the user to set the number K of latent factors, which is typically assumed to be independent of the sample size n. To address this problem, <cit.>, <cit.> and <cit.> proposed Bayesian nonparametric approaches that allow the number of latent factors to be estimated from the data, and to grow unboundedly with the size n of the matrix. In particular, <cit.> and <cit.>, considered a Poisson factorization model

    A_ij∼  (∑_k=1^∞ r_k v_ikv_jk ), 1≤ i,j≤ n

where the positive weights (r_k)_k≥ 1 represent the importance of community k, and v_ik>0 represents the level of affiliation of individual i to community k. <cit.> and <cit.>, extending work from <cit.>, assume that the weights (r_k) are the jumps of a gamma process, ensuring the sum in equation (<ref>) is almost surely finite. Using properties of Poisson random variables, the model (<ref>) can be equivalently represented as

    A_ij   =∑_k=1^∞ Z_ijk
    
    Z_ijk   ∼(r_k v_ikv_jk), k=1,2,…

for 1≤ i,j≤ n. The latent count variables Z_ijk may be interpreted as the number of latent interactions between two individuals i and j via community k, the overall number A_ij of interactions being the sum of those community interactions. For example, two members of the same company who also play sport together may meet five times at the company, and twice at the sport center, resulting in seven interactions overall. 
The overall number

    K_n=∑_k=1^∞∑_1≤ i, j≤ n Z_ijk>0

of communities k that generated at least one interaction between the n individuals is termed the number of active communities. For the gamma process Poisson factor model <cit.>, the number of active communities K_n grows logarithmically with the number n of individuals. The logarithmic growth assumption may be too restrictive. For example, the number of active communities may actually be unknown but bounded above; alternatively, it may increase at a rate faster or slower than logarithmic. 

In this paper, we consider generalizations of the gamma process Poisson factorization model, using completely random measures (CRM)  <cit.>. CRMs offer a flexible and tractable modeling framework <cit.>. The proposed models fit in the class of multivariate generalized Indian Buffet process priors recently developed by <cit.> and are also related to compound completely random measures <cit.>. We consider that (r_k) are the points of Poisson point process with mean measure ρ. Depending on the properties of this measure, the number of active communities K_n is either (i) bounded, with a random upper bound, (ii) unbounded and grows sub-polynomially (e.g. log n or loglog n) or (iii) unbounded and grows as n^2 σ, for some σ∈(0,1). For the implementation, we focus in particular on the generalized gamma process <cit.> where a single parameter flexibly controls all three behaviors.

The article is organized as follows. In Section <ref>, we describe the statistical model for count and binary matrices. The asymptotic properties of the model are derived in Section <ref>. In particular, we relate the asymptotic growth of the number of active features to the regular variation properties of the measure ρ. In Section <ref> we derive a Markov chain Monte Carlo algorithm for posterior inference that does not require any approximation to the original model. In Section <ref> we consider applications of our approach to overlapping community detection and link detection in networks, considering real network data with up to tens of thousands of nodes.



§ STATISTICAL MODEL FOR COUNT AND BINARY DATA





 §.§ General construction


We present here the model for directed count or binary observations, but the model can be straightforwardly adapted to undirected interactions. Let (r_k)_k=1,2…, be the points of a Poisson point process with σ-finite mean measure ρ on (0,∞), and assume that v_ik, i=1,…,n, k≥ 1, are independent and identically distributed from some probability distribution F on ℝ_+=[0,∞). The variable v_ik can be interpreted as the level of affiliation of an individual i to community k, and r_k to the importance of that community.

 For count data (A_ij), where A_ij denotes the number of directed interactions from node i to node j,  we consider the Poisson factorization model

    A_ij| (r_k,v_ik)∼  (∑_k=1^∞ r_k v_ikv_jk ), 1≤ i,j≤ n.

Denoting Λ_ij=∑_k=1^∞ r_k v_ikv_jk the Poisson rate for A_ij, the n× n rate matrix Λ^(n)=(Λ_ij)_1≤ i,j≤ n admits the following factorization as an infinite sum of rank-1 matrices

    Λ^(n)=∑_k=1^∞ r_k v_1:n,kv^⊺_1:n,k

where v_1:n,k=(v_1k,…,v_nk)^⊺. For the model to be well specified, the sum in the right-handside of Equation (<ref>) needs to be almost surely finite. A necessary and sufficient condition is

    ∬(1-e^-rv^2)ρ(dr)F(dv)<∞   and   ∬(1-e^-rv_1 v_2)ρ(dr)F(dv_1)F(dv_2)<∞.

A sufficient set of conditions[The sufficientness follows from the bound (<ref>) given in Appendix.], which we will assume to hold in the rest of this article, is that ρ is a Lévy measure and F has finite second moment, that is

    ∫_0^∞(1-e^-r)ρ(dr)   <∞   and A1
    ∫_0^∞ v^2F(dv)   <∞.A2

In this case, the community affiliations and weights for n nodes can be conveniently represented by a completely random measure

    G=∑_k≥ 1 r_kδ_v_1:n,k

on ℝ_+^n with mean measure ρ(dr)F^⊗^n(dv_1,…,dv_n) where F^⊗^n denotes the nth product measure of F; see <cit.> and <cit.> for background on CRMs and their applications. If the Lévy measure is finite, that is, if

    ∫_0^∞ρ(dr)<∞

then the number of points (r_k), and therefore the number of communities, is almost surely finite. Otherwise, when ∫ρ(dr)=∞, the number of communities is infinite. 

When we have binary observation (Y_ij), we treat the count matrix (A_ij) as a latent variable, and consider that Y_ij=1_A_ij>0 as in <cit.>. Integrating out (A_ij), this leads to the following model for binary observations

    Y_ij| (r_k,v_ik)∼ (1-exp  [∑_k=1^∞ r_k v_ikv_jk ] ), 1≤ i,j≤ n.




 §.§ Specific model


In the inference and experimental part, we use the following choice for the ρ and F. The Lévy measure ρ is taken to be that of a generalized gamma process (GGP, see <cit.>, <cit.>)

    ρ(dr) = κ/Γ (1-σ_0) r^-1-σ e^-τ r dr

where σ_0∈ (-∞, 1), κ>0 and τ>0. When σ_0=0, we obtain a gamma process, and the model corresponds to that of <cit.>. When σ_0 < 0, the Lévy measure is finite, while when σ_0 ≥ 0, the Lévy measure is infinite.

Concerning the affiliations, we will assume that F is a gamma distribution with parameters α>0 and β>0. That is, the probability density function (pdf) f is given by


    f(v) = β^α/Γ(α) v^α-1 e^-β v

where Γ denotes the usual gamma function. The hyperparameters (κ,σ_0,τ,α,β) and (κ^'=κ/ β^2σ_0, σ_0, τ^'=τβ^2,α,1) induce the same distribution for the latent factors (Λ_ij). In order to guarantee the identifiability of the hyperparameters, we therefore set β=1.






 §.§ Related work


The model introduced in this section can be seen from different perspectives that nicely connect it to the existing literature. First, the model can be seen as obtained from a functional of a CRM. Recall the definition of the CRM G in Eq. (<ref>). Define the n× n matrix Λ^(n) as the following functional of G

    Λ^(n)=∫_(0,∞)^n h(u)G(du)=∑_k≥ 1 r_k v_1:n,k v^⊺_1:n,k

where h(u)=u u^⊺.
Alternatively, this can be interpreted in the framework of compound completely random measures <cit.>. For each 1≤ i,j≤ n, denote G_ij=∑_k≥ 1 r_k v_ikv_jkδ_ζ_k where ζ_k are some community locations in some domain Θ, iid from some distribution H, irrelevant here. Then (G_ij)_1≤ i,j≤ n are compound CRMs on Θ and Λ_ij=G_ij(Θ). In the same vein, the model can also be interpreted as an instance of the class introduced by <cit.>. Denote Z^(n)_k  the n× n matrix with entries Z_ijk. Then the matrix-valued process ∑_k≥ 1 Z_k^(n)δ_ζ_k is a draw from a multivariate Indian buffet process.

Finally, as mentioned in the introduction, the model admits as a special case the Poisson factorization based on the gamma process of <cit.>.






§ ASYMPTOTIC PROPERTIES



In this section we study the asymptotic properties of the proposed class of models, and in particular the growth rate of the number of active communities as the sample size n grows, and the asymptotic proportion of communities of a given size. For a given sequence (r_k)_k≥ 1 and (v_ik)_i≥ 1,k≥ 1, denote A^(n)_ij and Z^(n)_ijk where n≥ 1, 1≤ i,j≤ n, k≥ 1 respectively the number of directed interactions and the number of community directed interactions distributed from Equations (<ref>) and (<ref>). We consider two different asymptotic settings


  * Constrained setting. For any 1≤ m ≤ n, and 1≤ i,j ≤ m, A_ij^(n)=A_ij^(m). In this setting, we suppose that the connections between the already observed nodes remain unchanged. It is equivalent to assuming that there is an infinite but  fixed graph and A^(n) represents the connections between the n first nodes of that graph.

  * Unconstrained setting. This setting is more general, and we only assume that A^(n)_ij and Z^(n)_ijk are marginally sampled from Equations (<ref>) and (<ref>).


All the results of this section, otherwise stated, hold for the unconstrained setting. We indicate when a stronger result holds in the constrained setting. All proofs are given in Appendix <ref>.



 §.§ General model


Let d^(n)_k be the degree of the community/feature k, corresponding to the number of interactions amongst n individuals due to community k, and defined as

    d^(n)_k = ∑_1≤ i,j≤ n Z^(n)_ijk.


A community is active if d^(n)_k≥ 1. The number of active communities is therefore defined as

    K_n = ∑_k=1^∞d^(n)_k ≥ 1

Denote K_n,j the number of communities with degree j≥ 1

    K_n,j = ∑_k=1^∞d^(n)_k = j


Note that under the constrained setting, d^(n)_k, K_n and ∑_ℓ≥ jK_n,ℓ are all almost surely increasing with the sample size n, whereas this is not necessarily the case for the unconstrained setting.


 
Under Assumptions (<ref>) and (<ref>), the number of active communities K_n is a Poisson random variable with mean

    Ψ(n) = ∫ (1-e^-r (∑_i=1^n v_i)^2  )  [∏_i=1^n F(dv_i)  ] ρ(dr) < ∞.

  The number K_n,j of communities with degree j is also Poisson distributed, with mean
  
    Ψ_j(n)=1/j!∫ r^j  ( )^2j e^-r  ( )^2  [∏_i=1^n F(dv_i)  ] ρ(dr).


  Finally, for j≥ 1, ∑_k ≥ j K_n,k, the number of communities with degree at least ℓ, is also Poisson distributed with mean ∑_k ≥ jΨ_k(n).







In the rest of the section we relate the asymptotic behavior of quantities of interest to the properties of the mean measure ρ. Let consider the tail Lévy intensity defined as

    ∀ x > 0, ρ(x) = ∫_x^∞ρ(dr).

We assume that ρ is a regularly varying function at 0, that is

    ρ(x) ∼ x^-σℓ(1/x) as x→ 0 A4

where σ∈[0,1) and ℓ is a slowly varying function verifying lim_t → +∞ℓ(at)/ℓ(t)  = 1 for all a>0. Examples of slowly varying functions include functions converging to a constant, log^a t for any t, loglog t, etc. Note that the CRM is finite activity if and only if σ=0 and ℓ(t)→ C<∞.


Now, let us consider the asymptotic behavior of the number of active communities K_n.


Let K_n be the number of active communities.
Then for 0≤σ < 1,

    𝔼[K_n] ∼Γ(1-σ) m_f^2σ n^2σℓ(n^2)

as n tends to infinity, where m_f = ∫ v F(dv). Additionally, for 0 < σ < 1,

    K_n ∼𝔼[K_n]   a.s.

If we further assume that the sequence (K_n)_n≥ 1 is almost surely non-decreasing (as in the constrained setting), then (<ref>) holds for σ = 0 and ℓ(t)→∞ as well. In the finite activity case, that is σ=0 and ℓ(t)→ρ(0)=∫_0^∞ρ(dr)<∞, we have

    K_n→ K_∞

as n tends to infinity, where K_∞ is a Poisson random variable with mean ρ(0). The above convergence holds in distribution for the unconstrained setting and almost surely for the constrained setting.





Let K_n,j be the number of communities of degree j. Then for 0 < σ < 1 and any j≥ 1,


    K_n,j∼σΓ(j-σ) /j! m_f^2σ n^2 σℓ(n^2)  a.s.

as n tends to infinity. Therefore,

    K_n,j/K_n→σΓ(j-σ)/Γ(1-σ) j!  a.s.

as n tends to infinity. This corresponds to a power-law behavior as

    σΓ(j-σ)/Γ(1-σ) j!∼σ/j^σ+1

for large j. If we further assume that for all k ≥ 1, (∑_j ≥ k K_n,j )_n≥ 1 is non-decreasing (constrained setting), then (<ref>) holds also for σ = 0 and ℓ(t)→∞.


Finally, let c^(n)(k,k') denote the cosine between the corresponding affiliation vectors

    c^(n)(k,k') = ∑_i=1^n v_ikv_ik'/√(∑_i v_ik^2)√(∑_i v_ik'^2) .

This coefficient gives a measure of the overlap between two communities k and k'. By the law of large numbers, for any k≠ k',

    c^(n)(k,k')∼(∫ v F(dv))^2/∫ v^2 F(dv)   a.s.  as n→∞.





 §.§ Specific case of the GGP


In the case of the GGP, we have

    ρ(x)=κτ^σ_0Γ(-σ_0,τ x)/Γ(1-σ_0)∼{[          -κτ^σ_0/σ_0             if σ_0<0;            κlog(1/x)             if σ_0=0; κ x^-σ_0/σ_0Γ(1-σ_0)             if σ_0>0 ] .

as x tends to 0, where Γ(a,x) is the incomplete gamma function. Note that ρ(x) is of the form x^-σℓ(1/x) where σ=max(0,σ_0) and

    ℓ(t)={[   -κτ^σ_0/σ_0      if σ_0<0;       κlog(t)      if σ_0=0; κ/σ_0Γ(1-σ_0)      if σ_0>0 ] .

is a slowly varying function at infinity. The results of the previous subsection therefore apply. For simplicity, we state the results for the constrained setting. We have, almost surely as n→∞

    K_n∼{[                K_∞           if σ_0<0;           2κlog(n)           if σ_0=0; κα^2σ_0n^2σ_0 /σ_0           if σ_0>0 ] .

where K_∞∼(-κτ^σ_0/σ_0). Additionally, for σ≥ 0,

    K_n,j/K_n→σ_0 Γ(j-σ_0)/Γ(1-σ_0) j!

almost surely as n→∞. Finally,

    c^(n)(k,k')→α/α +1.


Therefore, σ_0 governs the asymptotic behavior of the number of active communities. K_n is bounded with a random upper bound (σ_0<0), increases logarithmically (σ_0=0) or polynomially (σ_0>0). In the polynomial case, σ_0 also controls the power-law exponent of the proportion of communities of a given size. The parameter κ is an overall linear scaling parameter. Finally, the parameter α governs the amount of  overlapping between two communities.



§ SIMULATION, POSTERIOR CHARACTERIZATION AND INFERENCE




In this section we describe the marginal distribution and conditional characterization of the model. Building on these, we derive an exact sampler for simulating from the model, and a Markov chain Monte Carlo algorithm to approximate the posterior distribution. Importantly, the sampler targets the distribution of interest and does not require any truncation or approximation. For simplicity of exposition, we assume that ρ and F are absolutely continuous with respect to the Lebesgue measure, with ρ(dr)=ρ(r)dr and F(dx)=f(x)dx.


 §.§ Marginal distribution and simulation


For a fixed n, recall that K_n denotes the number of active communities. Let ((r_1,v_1:n,1),…, (r_K_n,v_1:n,K_n)) be the subsequence of (r_k,v_1:n,k) such that community k is active, meaning that ∑_1≤ i,j≤ n Z_ijk≥ 1, arranged in random order. Let Z_ijk be the number of community interactions corresponding to the active community (r_k, v_1:n,k). Note that

    A_ij=∑_k=1^K_nZ_ijk.

Let Z_k=(Z_ijk)_1≤ i,j≤ n. Using Proposition 5.2 of <cit.>, we obtain the following lemma.



The joint distribution of (K_n, (r_1:K_n, v_1:n,1:K_n), (Z_k)_k=1,…,K_n) is given by

    K_n    ∼(Ψ(n))

where Ψ(n) is defined in Eq.(<ref>), and

    p((r_1:K_n, v_1:n,1:K_n)|K_n)   =∏_k=1^K_n p(r_k,v_1:n,k|K_n)

where

    p(r_k,v_1:n,k|K_n)∝ (1-e^-r_k (∑_i=1^n v_ik)^2)ρ(r_k)∏_i=1^n f(v_ik).

Finally, for each k=1,…,K_n,

    Z_k|(r_1:K_n, v_1:n,1:K_n)   ∼(r_k v_1:n,kv^⊺_1:n,k )

where (Λ) denotes the distribution of a integer-valued matrix with Poisson entries with mean values Λ_ij, conditionally on the sum of the entries being strictly positive. This has probability mass function

    p(A)=
    {[ (1-e^-∑_ijΛ_ij)^-1∏_1≤ i,j≤ nΛ_ij^A_ije^-Λ_ij/A_ij!                                      if ∑_ij A_ij>0;                                                   0                                           otherwise;                                                     ].



The model has an infinite number of parameters, but Lemma <ref> allows us to derive an algorithm to exactly sample from it, by successively simulating K_n, (r_1:K_n, v_1:n,1:K_n), (Z_k)_k=1,…,K_n and A using Equations (<ref>), (<ref>), (<ref>) and (<ref>).


Sampling from the conditional distribution (<ref>) can be done efficiently by first sampling the number of multiedges ∑_i,jZ_i,j,k from a truncated Poisson with mean r_k (∑_i v_i,k)^2, then sampling iid the end nodes of the edges proportionally to the affiliation vector. Simulating from the conditional distribution  (<ref>) can be more challenging since it requires sampling a n+1 dimensional vector. However, if we suppose that the affiliations are Gamma distributed,   the problem reduces to sampling (r_k,∑_i v_i,k), which is a two dimensional vector, and independently sample the normalized affiliations from a Dirichlet distribution. Indeed, if the affiliations are Gamma distributed, we consider the following change of variable.

    ς_k   =∑_i=1^n v_ik,   k=1,…,K_n
    φ_ik   =v_ik/ς_k,    k=1,…,K_n; i=1,…,n

This gives the following algorithm for exact simulation from the model.


  * Sample K_n from Eq. (<ref>)

  * For k=1,…,K_n


  * Sample (φ_1k,…,φ_nk)∼(α,…,α)

  * Sample ς_k from

    p(ς)∝ψ(ς^2)(ς;nα,β)


  * Sample r_k|ς_k from

    p(r|ς)∝ (1-e^-rς^2) ρ(r)


  * Sample Z_k^(n) from Eq. (<ref>)


  * For 1≤ i,j≤ n, set A_ij=∑_k=1^K_nZ_ijk

where ψ(t)=∫_0^∞ (1-e^-wt)ρ(dw) is the Laplace exponent, (α,…,α) denotes the standard Dirichlet distribution and (x;a,b) denotes the probability density function of a Gamma random variable with parameters a and b, evaluated at x. In the case of the GGP, the Laplace exponent is

    ψ(t)=κ/σ((t+τ)^σ - τ^σ).

One can sample from Eq. (<ref>) and (<ref>) using rejection.



 §.§ Posterior characterization


Using Proposition 5.1 in <cit.>, one can characterize the conditional distribution of the CRM G given the latent community counts Z_ijk.


Conditionally on (Z_k^(n))_k=1,…,K_n, the CRM G has the same distribution as

    G^' + ∑_k = 1^K_nr̃_k δ_ṽ_1:n,k

where G^' is an inhomogeneous CRM on ℝ_+^n with mean intensity 
    e^-r (∑_i=1^n v_i)^2ρ(r)∏_i=1^n f(v_i)
 and (r̃_k,ṽ_1:n,k)_k=1,…,K_n are independent of G^' and iid with density

    p(r̃_k,ṽ_1:n,k|Z_k^(n))=e^-r_k (∑_i v_ik)^2 r_k^d_kρ(r_k) ∏_i=1^n v_ik^m_ikf(v_ik)

where m_ik=∑_j Z_ijk+Z_jik and d_k=∑_i,jZ_ijk.




In the case where f is a gamma pdf, we can use the same reparameterization as in the previous subsection with (ς_k, φ_1:n,k) in place of v_1:n,k. This leads to the following conditional distributions.

    ϕ_1:n,k|Z_k^(n)   ∼(α+m_1k,…,α+m_nk)
    
    p(ς_k|Z_k^(n))   ∝ϰ(d_k,ς^2) (ς;nα+2d_k,β) 
    
    p(r_k|ς_k,Z_k^(n))   ∝ e^-r_kς_k^2r_k^d_kρ(r_k)

where ϰ(m,t)=∫_0^∞ r^m e^-rtρ(r)dr. In the GGP case, we have

    ϰ(m,t)=κΓ(m-σ)/Γ(1-σ)(t+τ)^σ-m

and

    r_k|ς_k,Z_k^(n)∼(d_k-σ, ς_k^2+τ).




 §.§ Slice sampler for posterior inference


We recall that θ denote the set of hyperparameters of the mean measure ρ and pdf f. To simplify the presentation, here we suppose that we observe the complete adjacency matrix A, which means that we observe a directed and weighted graph with no missing (hidden) edge. The objective is obtain samples distributed from the conditional distribution

    p(K_n,(r_k,v_1:n,k)_k=1,…,K_n,θ| A).

In the Appendix, we show how to do inference when we only observe a partial graph (with missing edges to predict) that can be directed or undirected, weighted or binary. In order to leverage the Poisson factorization construction, we augment the model with the latent community counts Z_k. Additionally, to deal with the unknown number of active communities K_n, we use auxiliary slice variables, similarly to other Gibbs sampler for Bayesian nonparametric models <cit.>. For each directed pair (i,j) such that A_ij≥ 1 consider the scalar latent variable

    s_ij|(r_k,Z_ijk)_k=1,…,K_n∼ (0, min_{k|Z_ijk≥ 1}r_k )

and denote s=min_ij s_ij. Note that by definition, r_k≥ s for all k=1,…,K_n. Let

    G   =∑_k r_kδ_v_1:n,kr_k≥ s:=∑_k=1^K_nr_kδ_v_1:n,k

 be the CRM corresponding to the set of active or inactive communities with weight r_k≥ s, of (almost surely finite) cardinality K_n≥ K_n. Denote Z_ijk≥ 0 the associated community interactions, and Z_k=(Z_ijk). The data augmented slice sampler draws samples asymptotically distributed from

    p((Z_k)_k=1,…,K_n,G,θ, s | A).


The main steps of the algorithm are as follows.


  * For each directed pair (i,j) such that A_ij≥ 1, Update (Z_ijk)_k=1,…,K_n given the rest of the variables,

  * Update the hyperparameters θ given the rest  of the variables,

  * Update (G, s) given the rest of the variables.

The details of each step are given in Appendix <ref>. Each iteration of the Gibbs sampler has a time complexity scaling in K_n S where S is the number of nonzero entries of the matrix. Therefore, the algorithm takes advantage of the sparsity of the networks. Additionally, each entry of the sparse graph can be dealt with independently, making the algorithm straightforwardly parallelizable.





§ EXPERIMENTS



We implement the algorithm described in the previous section with the GGP-Gamma scores model. We assign Gamma priors on the hyperparameters κ, τ, α with parameters (0.1,0.1). We fix β = 1. We allow up to a linear growth of the number of communities, corresponding to σ<  0.5, for small datasets and use a Gamma prior with parameter (0.1,0.1) on 1-2σ. For larger datasets, we restrict σ < 0.25, meaning that the number of communities cannot grow at a faster rate than √(n). This is obtained by using a Gamma prior with parameter (0.1,0.1) on 1-4σ.



 §.§ Synthetic dataset

We first run the algorithm on a synthetic dataset simulated from our model, to check that the algorithm can recover the true parameters. We sample a directed and unweighted graph from the GGP-gamma model with size n=800 and σ = 0.2,κ=1,τ=0.15,α=0.05,β=0.2. The number of edges of the obtained graph is 20198, and the true number of active communities is 42.
We run three chains in parallel with 500,000 iterations, with 250,000 iterations for burn-in. We show in Figure <ref> trace plots of the number of active communities K_n and parameter σ showing the MCMC algorithm can recover these parameters.







 §.§ Political blogs

The polblogs network <cit.> is the network of the American political blogosphere in February 2005. It is a directed unweighted graph, where there is an edge (i,j) if blog i cites blog j. It is composed of 1490 nodes and 19025 edges. For each node, some ground truth information about its political affiliation (republican/democrat) is known. 

We will use this dataset in order to illustrate the role of the parameter α in the model. As indicated in Section <ref> this parameter tunes the amount of overlapping between the communities. A smaller value enforces less overlap between communities. We run three chains with 500,000 iterations. The posterior samples of σ for three different values of α are in also shown in Figure <ref>. The model allows overlapping communities but, for visualization purposes, it is useful to obtain an associated partition of the nodes. For each iteration, one can cluster the nodes by assigning each node to the community where it is most active. That is, at iteration t of the MCMC algorithm, define for i=1,…,n

    c_i^(t)=_k{√(r^(t)_k) v^(t)_ik}

the cluster membership of node i. We then compute an approximate Bayesian point estimate c=(c_1,…,c_n) of the partition of the nodes, using Binder's loss function <cit.>. Nodes are reordered according to their estimated membership c, and Figure <ref> shows the densities of connection between and within clusters for three different values of α.  Depending on the amount of overlapping, we obtain two (α=0.8), three (α=0.4) or four (α=0.2) communities. In order to interpret those communities, we calculate in Table <ref> for each community the proportion of interactions between democrat blogs, between a democrat and a republican blog, and between two republican blogs. For α=0.8, there are two estimated communities which can clearly be identified as democrat (community #1) and republican (community 2). For α=0.4, we have three communities. One is mostly associated to democrat blogs (#1) while the other two correspond to a split of the republican blogs into right (#2) and center-right (#3) groups. For α=0.2, we obtain a further split of the democrat blogs into left (#1) and center-left (#2) groups. Increasing the value of α therefore leads to a finer and finer partition of the nodes.












 §.§ Wikipedia topcast

The network is a partial web graph of Wikipedia hyperlinks collected in September 2011 <cit.>. It is a directed unweighted graph where an edge (i,j) corresponds to a citation from a page i to page j. We restrict it to the first 3000 nodes, and the associated 5687 edges.  We run three MCMC chains for  200,000 iterations. Trace plots of the number of active communities and parameter σ are given in Figure <ref>. Figure <ref> shows the adjacency matrix reordered by communities, as explained in the previous section. In order to check that the learnt communities/features are meaningful, we report in Figures the proportions of webpages associated to a given category within a given community/feature (note that a webpage can be associated to multiple categories hence the proportion do not sum to 1).




Note that, while the approach is able to estimate the latent block-structure, this dataset has the particularity of having star nodes, a feature that is not captured by our model.










 §.§ Deezer

The dataset was collected from the music streaming service Deezer in November 2017 <cit.>. It represents the friendship network of a subset of Deezer users from Romania. It is an undirected unweighted graph where nodes represent the users and edges are the mutual friendships. There are 41 773 nodes and 125 826 edges. We run three chains with 100 000 iterations each. Posterior histograms of the number of active communities and σ are given in Figure <ref>. The algorithms finds around 45 communities/features for this dataset. The reordered adjacency matrix and block densities based on the point estimate of the partition are given in Figure <ref>.



Now we can reorder the nodes using approximate MAP clustering as previously. We obtain the following adjacency matrix


For each individual in the network, a list of musical genres liked by that person are available. There are in total 84 distinct genres. We represent in Figure <ref> the proportion of individuals who liked a subset of the 84 genres for three different communities where the interpretation in terms of genres is quite clear. The overall proportion of individuals liking a given genre is shown at the bottom of Figure <ref>. If the bar is red, this indicates that the proportion is 10% higher in the community than in the population. If the bar is blue, this means the proportion is 10% lower. Community 11 can be interpreted as R&B, Community 8 as Dance, and Community 3 as Rock music. For some of the communities, not reported here, the interpretation in terms of the liked genres is less clear, and may be due to other covariates. 







§ DISCUSSION


The model presented in this paper assumed the same parameter β for each node. We can also consider a degree corrected version of the model, similarly to <cit.>, where each node is assigned a different parameter β_i > 0 and then defining Z_ijk∼(r_k v_ik v_jk/β_i β_j). It is unclear however if a MCMC sampler targeting the exact posterior distribution could be implemented, and one may need to resort to some truncation approximation as in <cit.>.

The count matrix (A_ij) is infinitely exchangeable, hence the model presented in this article lead to asymptotically dense graphs. That is, ∑_1≤ i,j≤ n A_ij≍ n^2 as n tends to infinity. In order to obtain sparse graphs, we could consider two different strategies. The first solution consists in dropping the infinite exchangeability property and take β^(n)_i →∞ with n, then the number of edges will behave as (n/β^(n))^2 (we can for instance take β^(n)_i = √(n) for any node i to obtain a linear growth of the number of edges). The  model would still be finitely exchangeable for any fixed n, but not projective anymore. The second solution would be to consider the different notion of infinite exchangeability developed in <cit.> and consider (β_i)_i as a realization of a Poisson point process.

Finally, we presented a model for count (and binary) data. The results build on the additive contributions of the communities, which is why we chose the Poisson distribution on the entries of the adjacency matrix (A_ij). We can generalize to non count data using other probability distributions which are closed under convolution. For example, one could consider A_ij∼(∑_k r_k v_ik v_jk,1) for A_ij∈ℝ_+ or A_ij∼𝒩(∑_k r_k v_ik v_jk,1) for A_ij∈ℝ.






plainnat











§ PROOFS




 §.§ Technical Lemmas



(<cit.>, Propositions 17 and 19) Let ρ be a Lévy measure, let ρ(x) = ∫_x^∞ρ(r) dr be the tail Levy intensity and ψ (t) = ∫ (1-e^- r t) ρ(dr) its Laplace exponent. Then the two following conditions are equivalent:

    ρ(x)    x → 0_+∼   ℓ(1/x) x^- σ
    ψ(t)    x → +∞∼   Γ (1-σ)  t^σℓ(t)

with ℓ a slowly varying function and 0 ≤σ < 1.


Besides, if we let ψ_d(t) = t^d/d!∫ r^d e^-rtρ(dr)


  * if σ > 0, then (<ref>) implies that ψ_d(t) t → +∞∼σΓ (d-σ)/d! t^σℓ(t)

  * if σ = 0, then (<ref>) implies that ψ_d(t) = o(ℓ(t))




<cit.>
Let X be a Poisson random variable with parameter λ . For any t>0

    ℙ(|X-λ| ≥λ t) ≤ 2 e^-λ t^2/2(1+t).











































 
Let (X_n)_n≥ 1 be a sequence of Poisson random variables with mean (μ_n)_n≥ 1.  If log n = o(μ_n) then X_n ∼μ_n almost surely as n tends to infinity.









Let 0< ϵ < 1/2. Using Lemma <ref>, we have

    ℙ (  |X_n/μ_n -1  | ≥ϵ )    ≤ 2 e^-ϵ^2 μ_n/4
       =2 n^-ϵ^2 μ_n/4 log n





Using the assumption, we have that -ϵ^2  μ_n/4 log n→ -∞. Therefore, the RHS of (<ref>) is summable. The almost sure result follows from Borel-Cantelli lemma.



For any x,y≥ 0, we have the following bound

    1-e^-xy≤max(1,y)(1-e^-x)




The bound is trivial when y≤ 1. Consider the case y>1. For all x, the function y→ e^-xy is convex hence f_x(y)=e^-xy-1/y is a monotonically non-decreasing function of y therefore

    e^-xy-1/y≥ e^-x-1

for y≥ 1.



























 §.§ Proofs of Section 3





The result for K_n is proved in <cit.> in the general context of GIBP. We provide here the details of the proof for K_n, which can be straightforwardly adapted to K_n,j.

First, let us remark that the bound (<ref>) together with assumptions (<ref>) and (<ref>) imply Ψ(n)<∞. For s < 0,

    𝔼 [e^s K_n]    =   𝔼 [ ∏_k=1^∞𝔼 [e^s 1_∑_1≤ i,j ≤ n Z_ijk≥ 1 | G]  ] 
       =   𝔼 [ ∏_k=1^∞ [ e^s + (1-e^s) e^-r_k (∑_i=1^n v_ik)^2 ] 
       =   𝔼 [ e^∑_k log [e^s + (1-e^s) e^-r_k (∑_i=1^n v_ik)^2] ]

Then, since

    log [e^s + (1-e^s) e^-r_k (∑_i=1^n v_ik)^2]    =   log [e^-r_k (∑_i=1^n v_ik)^2 + e^s (1- e^-r_k (∑_i=1^n v_ik)^2)] 
       ≤   log [1 + e^s (1- e^-r_k (∑_i=1^n v_ik)^2)] 
       ≤    1 - e^-r_k (∑_i=1^n v_ik)^2

and the last part is integrable, we can use Campbell's theorem <cit.> to get:

    𝔼 [e^s K_n]    =   exp [ ∫ (e^log [e^s+ (1-e^s)e^-r (∑_i=1^n v_i)^2]  -1 ) ∏_i=1^n F(dv_i) ρ(dr)   ] 
       =   exp [ (e^s-1) Ψ(n)  ]

We can prove similarly that K_n,d is a Poisson random variable with mean Ψ_d(n) and that ∑_d ≥ D K_n,d is Poisson distributed with mean ∑_d ≥ DΨ_d(n). The assumption Ψ_d(n) < ∞ is also sufficient in this case to apply Campbell's theorem.


























From Proposition <ref>, we get that

    𝔼(K_n) = ∫ (1-e^-r (∑_i=1^n v_i)^2 ) ∏_i=1^n F(dv_i)  ρ(dr)  = Ψ(n)

Let (V_i)_i ∈ℕ be i.i.d random variables with distribution F. By assumption, 0 < 𝔼[V_i] = μ < +∞ and 𝕍ar [V_i] = τ^2 < +∞. Let ϵ > 0. Let A(r) be defined for r > 0 by

    A(r) = 𝔼 [ 1 - e^-r (∑_i=1^n V_i)^2 ]

Since v↦ 1 - e^-r v is concave, using successively Jensen's inequality and the independence of (V_i), we obtain

    A(r)    ≤    1 - e^-r 𝔼[(∑_i=1^n V_i)^2]
       ≤    1 - e^-r (n^2 μ^2 + n τ^2)
       ≤    1 - e^- (1+ϵ) r n^2 μ^2

where the last inequality holds for any ϵ>0 when n>τ^2/ϵμ^2. Therefore, for ϵ>0 and n>τ^2/ϵμ^2

    Ψ(n) ≤∫ (1-e^-(1+ϵ) r n^2 μ^2) ρ(dr).


Besides, since v↦ 1 - e^-r v is increasing, by Markov's inequality we have for any ϵ>0

    A(r)    =   𝔼 [ 1 - e^-r n^2 (1/n∑_i=1^n V_i)^2 ] 
       ≥   ℙ (1/n∑_i=1^n V_i ≥μ/(1+ϵ) )  (1-e^-n^2 μ^2/(1+ϵ)^2 r  )


Hence

    ℙ (1/n∑_i=1^n V_i ≥μ/(1+ϵ) ) ψ (n^2 μ^2/(1+ϵ)^2 ) ≤ Ψ(n) ≤ψ  ( (1+ϵ)  n^2 μ^2 )

where ψ (t) = ∫ (1-e^- r t) ρ(dr) is the Laplace exponent. Furthermore, by the law of large numbers,

    ℙ (1/n∑_i=1^n V_i ≥μ/(1+ϵ) ) → 1 .

Therefore,  under Assumption (<ref>), Lemma <ref> implies

    𝔼(K_n) ∼Γ(1-σ)μ^2σ n^2σℓ(n^2)

as n tends to infinity.

In the finite-activity case, that is  σ=0 and ℓ(t)→ρ(0)<∞, we have 𝔼(K_n) →ρ(0) hence K_n tends in distribution to (ρ(0)).

 Now, for σ > 0, the almost sure result (<ref>) follows from Lemma <ref> and the fact that for every slowly varying function ℓ_0 and every ϵ>0

    lim_x →∞ℓ_0(x) x^-ϵ = 0.






Finally, assume that (K_n)_n≥ 1 is non-decreasing. We only need to prove the asymptotic behavior for σ = 0. In that setting, Ψ(n) ∼ℓ(n^2). Using the assumption that ∫ρ(dr) = ∞, we therefore have lim_n→∞Ψ(n) = ∞. Let n ≥ 1,

    Ψ(n+1) - Ψ(n) = ∫𝔼 [ e^-r (∑_i=1^n V_i)^2 - e^-r (∑_i=1^n+1 V_i)^2 ] ρ(r) dr

Since V_i ≥ 0 a.s, it comes that (Ψ(n))_n is non-decreasing.
Now extend the sequence (Ψ(n))_n to a non-decreasing and continuous function Ψ on ℝ_+ (by linear interpolation for instance). Let t > 1, then

    1 ≤Ψ(t+1)/Ψ(t)≤Ψ(3⌊ t ⌋)/Ψ(⌊ t ⌋)∼ℓ(9⌊ t ⌋^2)/ℓ(⌊ t ⌋^2)→ 1

Hence limΨ(t+1)/Ψ(t) = 1

Now, for every integer m ≥Ψ(0), choose t_m such that Ψ(t_m) = m. We have that (t_m) is non-decreasing and diverges. Since Ψ is increasing, it comes

    Ψ(t_m+1-1) ≤Ψ( ⌊ t_m+1⌋ - 1) ≤Ψ(t_m+1) = m+1

Hence, Ψ(⌊ t_m+1⌋ - 1) ∼ m.
Then, using Lemma <ref>, we get that

    K_⌊ t_m+1⌋ - 1∼Ψ( ⌊ t_m+1⌋ - 1) ∼ m    a.s.

Finally, let n ≥Ψ(0), let m_n = min{ m  |  n ∈{⌊ t_m⌋ ,...,⌊ t_m+1⌋ - 1 }},

    K_⌊ t_m_n⌋/Ψ( ⌊ t_m_n+1⌋ -1) ≤K_n/Ψ(n)≤K_⌊ t_m_n+1⌋ - 1/Ψ( ⌊ t_m_n⌋)

Since t_m_n→∞, both bounds converge to 1 almost surely, which gives the result.




As for Proposition <ref>, we only need to show that for d ≥ 1,

    Ψ_d(n) ∼σΓ(d-σ) /d! n^2 σμ^2 σℓ(n^2)

Therefore the proof is very similar to the one of Proposition <ref>. However, there are some technicalities we need to address since here v ↦v^d/d!e^-rv is neither convex nor decreasing. Like previously, we will lower bound and upper bound Ψ_d(n) by two quantities that are equivalent to σΓ(d-σ) /d! n^2 σμ^2 σℓ(n^2).

Let us first introduce some notations. Let (V_i)_i ∈ℕ i.i.d variables with distribution F. Let S_n = ∑_i=1^n V_i and φ_d(r,v) defined as

    φ_d(r,v) = r^d v^d/d! e^-r v.



Now, define A_d(r,n) for r > 0 by

    A_d(r,n) = 𝔼[φ_d(r,S_n^2) ]

Suppose that 0 < σ < 1, let ϵ > 0, recalling that 𝔼 V_i = μ, define B_ϵ = [nμ/ϵ + 1 , (1+ϵ)nμ ]. Let us notice that the law of large numbers gives us that ℙ (S_n ∈ B_ϵ) → 1.




  * Lower bound: We have that

    𝔼[1_S_n ∈ B_ϵφ_d(r, S_n^2) ] ≤ A_d(r,n)

Besides, for all v∈ B_ϵ and all r>0

    ϕ_d(r,v)≥r^d μ^2d n^2d/(1+ϵ)^2d d! e^-r n^2 (1+ϵ)^2 μ^2

hence

    ℙ (S_n ∈ B_ϵ) r^d μ^2d n^2d/(1+ϵ)^2d d! e^-r n^2 (1+ϵ)^2 μ^2≤𝔼[1_S_n ∈ B_ϵφ_d(r,S_n^2) ].

Therefore, using Lemma <ref> and since ℙ (S_n ∈ B_ϵ) → 1, we have for n large enough

    σΓ(d-σ) /d! (1+ϵ)^4d + 1 - 2σ n^2 σμ^2 σℓ(n^2) ≤∫ A_d(r,n) ρ(dr).




  * Upper bound: We have that

    A_d(r,n) = 𝔼[1_S_n ∈ B_ϵφ_d(r, S_n^2) ] + 𝔼[1_S_n ∉B_ϵφ_d(r, S_n^2)]

Like previously, since

    𝔼[1_S_n ∈ B_ϵφ_d(r,n S_n) ] ≤ℙ (S_n ∈ B_ϵ) r^d (1+ϵ)^2dμ^2d  n^2d/d! e^-r n^2 μ^2/(1+ϵ)^2

We find that for n large enough,

    ∫ A_d(r,n) ρ(dr) ≤σΓ(d-σ) (1+ϵ)^4d + 1 - 2σ/d! n^2 σμ^2 σℓ(n^2) + ∫𝔼[1_S_n ∉B_ϵφ_d(r, S_n^2)] ρ(dr)

Therefore, we only need to prove that

    ∫𝔼[1_S_n ∉B_ϵφ_d(r, S_n^2)] ρ(r)dr = o(n^2σℓ(n^2)).

In order to do so, we split the integral with respect to r in two parts, an integral over (0,1/n^2) and an integral over (1/n^2,∞) and show that both are o(n^2σℓ(n^2)). Since φ_d(r,v) ≤ 1,

    ∫_1/n^2^∞𝔼[1_S_n ∉B_ϵφ_d(r, S_n^2)] ρ(dr)    ≤ℙ (S_n ∉B_ϵ) ∫_1/n^2^∞ρ(dr) 
       = ℙ (S_n ∉B_ϵ) ρ(1/n^2)
       = o(n^2σℓ(n^2))


where the last line follows from the law of large numbers and Assumption (<ref>).
Besides,

    ∫_0^1/n^2𝔼[1_S_n ∉B_ϵφ_d(r, S_n^2)] ρ(dr)    =   ∫_0^1/n^2𝔼[1_S_n ∉B_ϵr S_n^2/dφ_d-1(r, S_n^2)] ρ(dr) 
       ≤   ∫_0^1/n^2𝔼[1_S_n ∉B_ϵ r S_n^2] ρ(dr) 
       =   𝔼[1_S_n ∉B_ϵS_n^2/n^2]∫_0^1/n^2 r n^2 ρ(dr) 
       ≤   𝔼[1_S_n ∉B_ϵS_n^2/n^2] e ∫_0^1/n^2 r n^2 e^-r n^2ρ(dr) 
       ≤    8𝔼[1_S_n ∉B_ϵS_n^2/n^2] σΓ(d-σ) /d! n^2 σℓ(n^2)

where the last inequality holds for n large enough by Assumption (<ref>) and  Lemma <ref>. Now, we have that

    S_n^2/n^2 < 1/n∑_i=1^n V_i^2

Since (V_i^2)_i are i.i.d random variables in ℒ_1, we know that (1/n∑_i=1^n V_i^2)_n≥ 1 is uniformly integrable. Therefore, (1_S_n ∉B_ϵS_n^2/n^2)_n≥ 1 is uniformly integrable. Besides, using the law of large numbers, the sequence converges almost surely, and hence in probability, to 0. Therefore, lim_n 𝔼[1_S_n ∉B_ϵS_n^2/n^2] = 0, which concludes the proof.



For σ = 0, the previous computations for the upper bound give that almost surely, Ψ_d(n) = o(ℓ(n^2)) = o(Ψ(n)). Now let D > 1,

    𝔼∑_d ≥ D K_n,d = Ψ(n) - ∑_d=1^D-1Ψ_d(n) ∼ℓ(n^2)

And since x ↦∑_d≥ Dφ_d(1,x) is non-decreasing, (𝔼∑_d ≥ D K_n,d)_n is non-decreasing, therefore, similarly to the proof for σ = 0 for (K_n)_n, we find that

    ∑_d ≥ D K_n,d∼𝔼∑_d ≥ D K_n,d∼ℓ(n^2)   a.s

Therefore, we finally find that

    K_n,D/K_n = ∑_d ≥ D K_n,d - ∑_d ≥ D+1 K_n,d/K_n→ 0    a.s
















§ GIBBS SAMPLER





As mentioned in the main text, the observed graph can be directed or undirected,  binary or count, and can have missing entries we would like to predict. Denote by B the observed graph. Here we describes the steps of a Gibbs algorithm with stationary distribution

    p(K_n,(r_k,v_1:n,k)_k=1,…,K_n,θ| B).

Notice that observing the full matrix B=A corresponds to a weighted and directed graph with no missing entry. Let ℐ denote the set of all possible edges. In the directed case, ℐ = { (i,j) | 1≤ i,j ≤ n } and on the undirected case ℐ = { (i,j) | 1≤ i ≤ j ≤ n }. We say that (i,j) is not observed if we don't know the value of A_i,j. Remark that (i,j) can be observed and still A_i,j = 0. Denote 𝒪 the set of all observed entries and 𝒪^c = ℐ∖𝒪, the set on unobserved entry.  For all unobserved entry (i,j) ∈𝒪^c, set B_i,j = -1


Additionally, to deal with the unknown number of active communities K_n, we use auxiliary slice variables s_i,j for all (i,j) ∈ℐ, details are given in the following paragraphs. Denote s the smallest non-zero slice variable s_i,j for (i,j) ∈ℐ. By definition of the slice variables, r_k≥ s for all k=1,…,K_n. Let

    G   =∑_k r_kδ_v_1:n,kr_k≥ s:=∑_k=1^K_nr_kδ_v_1:n,k

be the CRM corresponding to the set of active or inactive communities with weight r_k≥ s, of (almost surely finite) cardinality K_n≥ K_n. Denote Z_ijk≥ 0 the associated community interactions, and Z_k=(Z_ijk).





 §.§ Directed graph



For each observed pair (i,j)∈𝒪, we define the slice variable as

    s_ij|(r_k,Z_ijk)_k=1,…,K_n∼ (0, min_{k|Z_ijk≥ 1}r_k )

if A_ij > 0 and s_ij = 0 otherwise. For each non observed entry (i,j) ∈𝒪^c, we define s_i,j by (<ref>) if { k | Z_ijk≥ 1}≠∅ and

    s_ij|(r_k,Z_ijk)_k=1,…,K_n∼ (0, 1)

otherwise.



  §.§.§ Gibbs sampler step 1 for weighted graph on observed entries



Updating (Z_k)_k=1,...,K_n| (s, G), θ, B on observed entries indexes 


We sample (Z_l)_l=1,..,K_n associated to all atoms of G and keep only the non empty communities. For every (i,j) ∈𝒪 such that A_i,j > 0. define the random variable m_ij = min_{l|Z_ijl≥ 1}r_l. Then, writing the joint distribution it comes that independently for every such (i,j),

    ℙ((Z_ijl)_l=1,..,K_n | (s, G^+),θ, B_ij) ∝∏_i,j1/m_ij1_s_ij < m_ij((Z_ijl)_l; B_ij,(p_ijl)_l)

where  is the multinomial distribution and p_ijl = r_l v_il v_jl/∑_t=1^K_n r_t v_it v_jt. Let p_ijl = r_l v_il v_jl. To simplify the notations, let us suppose that the atoms of G are in decreasing order. Remark that the indexing of Z is different from the one of Z, the second corresponding to the one of the truncated random measure. For each observed edge (i,j) independently, we can proceed in 4 phases for this step.


  * Sample m_ij from the locations of G such that ℙ(m_ij = r_L) ∝(∑_l=1^L p_ijl)^B_ij - (∑_l=1^L-1 p_ijl)^B_ij/r_L1_s_ij < r_L.

  * For l > L, set Z_ijl = 0

  * Sample Z_ijL∼(B_ij,p_ijL/∑_l=1^L p_ijl), where  is the zero truncated binomial distribution

  * Sample (Z_ij1,..,Z_ijL-1) ∼(B_ij-Z_ijL,(p_ijl/∑_t=1^L-1 p_ijt)_l ≤ L-1   )





  §.§.§ Gibbs sampler step 1 for unweighted graph on observed entries


In this setting we observe a binary matrix B_ij = 1_A_ij > 0. Then the first step of the Gibbs sampler is modified and becomes:


Updating (Z_k)_k=1,...,K_n| (s, G), θ, B on observed entries indexes


For each observed edge (i,j)∈𝒪 independently do


  * Sample m_ij from the locations of G such that

    ℙ(m_ij = r_L) ∝e^∑_k=1^L p_ijl-e^∑_k=1^L-1 p_ijk/r_L1_s_ij < r_L

Suppose m_ij = r_L

  * For l > L, set Z_ijl = 0

  * Sample Z_ijL∼(p_ijL), where  is the zero truncated Poisson distribution

  * For l < L, sample Z_ijl∼ Poiss(p_ijl)





  §.§.§ Gibbs sampler step 1 on unobserved entries


For each unobserved entry (i,j)∈𝒪^c, knowing s_ij, we define L_0 = max{ k |  r_k > s_ij}.



  * Draw 1_A_ij=0, which is a Bernoulli with parameter

    p = 1/1 + ∑_L=1^L_0e^∑_k=1^L p_ijk-e^∑_k=1^L-1 p_ijk/r_L


  * If A_i,j≠ 0, then use subsection <ref>. Otherwise, set all counts of that entry to zero





 §.§ Undirected graph


In the undirected graph, we suppose that for i ≠ j, B_ij = A_ij + A_ji and B_ii = A_ii. Besides, in this setting we actually don't need to sample Z_ijk for all (i,j,k) but only Z_ijk + Z_jik.
For each observed pair (i,j)∈𝒪, we define the slice variable as

    s_ij|(r_k,Z_ijk+Z_jik)_k=1,…,K_n∼ (0, min_{k|Z_ijk+Z_jik≥ 1}r_k )

if B_ij > 0 and s_ij = 0 otherwise. For each non observed entry (i,j) ∈𝒪^c, we define s_ij by (<ref>) if { k | Z_ijk + Z_jik≥ 1}≠∅ and

    s_ij|(r_k,Z_ijk+Z_jik)_k=1,…,K_n∼ (0, 1)

otherwise. Then Step 2 and 3 remain unchanged. For step 1, simply replace p_ijk by 2 p_ijk for i ≠ j.



 §.§ Proofs for the Gibbs sampler step 1



  §.§.§ Weighted graph

Here will give the posterior distribution of the count matrices and show that

    (v_k, Z_k)_k=1,...,K_n| (s, G), θ, B dist= (v_k, Z_k)_k=1,...,K_n| s, G, θ, B


In order to do so, we derive the RHS posterior distribution. Let us first notice that given G, sampling the non zero counts and corresponding locations is equivalent to sampling (Z_k)_k for k ∈ℕ.
As stated previously, we can treat each edge (i,j) independently. Therefore, we sample the sequence (Z_ijk)_k. Here we suppose that the communities come with decreasing activity order. Let the random variable L = max{k  |  Z_ijk > 0} (supposing that the (r_k)_k are decreasing). And let p_ijk = r_k v_ik v_jk


    ℙ((Z_ijk)_k| s, G, θ, A_ij)
       ∝   ℙ((Z_ijk)_k|G, θ, A_ij) ×ℙ(s_ij|(Z_ijk)_k, G, θ, A_ij)
       ∝   1_∑_k Z_ijk = A_ijA_ij!/∏_k=1^L Z_ijk!∏_k=1^L p_ijk^Z_ijk×1/r_L1_s_ij < r_L 
       ∝   (∑_k=1^L p_ijk)^A_ij (1-(∑_k=1^L-1 p_ijk/∑_k=1^L p_ijk)^A_ij  )/r_L1_s_ij < r_L
       ×   1_1 ≤ Z_ijL≤ A_ijA_ij!/Z_ijL!(A_ij-Z_ijL)! (p_ijL/∑_k=1^L p_ijk)^Z_ijL (∑_k=1^L-1 p_ijk/∑_k=1^L p_ijk)^A_ij-Z_ijL/1-(∑_k=1^L-1 p_ijk/∑_k=1^L p_ijk)^A_ij
       ×   1_∑_k=1^L-1 Z_ijk = A_ij - Z_ijL(A_ij-Z_ijL)!/∏_k=1^L-1 Z_ijk!∏_k=1^L-1 (p_ijk/∑_k=1^L-1 p_ijk)^Z_ijk

This shows how we can sample in three steps these variables. Let us remark that the second part corresponds to the distribution of a zero truncated binomial and that the third part corresponds to the distribution of a multinomial. We also notice that only the elements of G are actually needed.



  §.§.§ Unweighted Graph

We proceed similarly for the unweighted graph

    ℙ((Z_ijk)_k| s, G, θ, B_ij = 1)
       ∝   ℙ((Z_ijk)_k|G, θ, B_ij = 1) ×ℙ(s_ij|(Z_ijk)_k, G, θ, B_ij = 1)
       ∝   1_Z_ijL≠ 0∏_k=1^L p_ijk^Z_ijk/Z_ijk ! e^-p_ijk× e^∑_k=1^L p_ijk×1/r_L1_s_ij < r_L 
       ∝   e^∑_k=1^L p_ijk-e^∑_k=1^L-1 p_ijk/r_L1_s_ij < r_L
       ×   1_Z_ijL≠ 01/1-e^-p_ijLp_ijk^Z_ijL/ Z_ijL!  e^-p_ijL
       ×   1_Z_ij(L+1),... = 0∏_k=1^L-1p_ijk^Z_ijk/Z_ijk ! e^-p_ijk






  §.§.§ Prediction

Here we show how to update the missing entries we try to predict. Let us recall that for a predicted count, if it is positive, we define the slice variable as previously. However, if the count is equal to zero, then the slice variable is simply uniform over [0,1]. Now let L_0 = max{ k |  r_k > s_ij}


    ℙ((Z_ijk)_k| s, G, θ)
       ∝   ℙ((Z_ijk)_k|G, θ) ×ℙ(s_ij|(Z_ijk)_k, G, θ)
       ∝   1_Z_ij(L_0+1),... = 0 e^∑_k=1^L_0 p_ijk∏_k=1^L_0p_ijk^Z_ijk/Z_ijk ! e^-p_ijk× (1_A_ij≠ 01/r_L1_s_ij < r_L  + 1_A_ij = 0)

Now let

    f((Z_ijk)_k) = 1_Z_ij(L_0+1),... = 0  e^∑_k=1^L_0 p_ijk∏_k=1^L_0p_ijk^Z_ijk/Z_ijk ! e^-p_ijk× (1_A_ij≠ 01/r_L1_s_ij < r_L  + 1_A_ij = 0)


Using <ref>, it comes that

    𝔼 ( f((Z_ijk)_k) | A_ij≠ 0) = e^∑_k=1^L_0 p_ijk∑_L=1^L_0e^∑_k=1^L p_ijk-e^∑_k=1^L-1 p_ijk/r_L


Besides,

    𝔼 ( f((Z_ijk)_k) | A_ij = 0) = f(0) = 1


Therefore, here we proceed in two steps, first we sample the binomial 1_A_ij = 0 with parameter

    p = 1/1 + e^∑_k=1^L_0 p_ijk∑_L=1^L_0e^∑_k=1^L p_ijk-e^∑_k=1^L-1 p_ijk/r_L


Then, conditioning on the event A_ij≠ 0, we use <ref> to proceed.



 §.§ Proof for the Gibbs step 2


Here we show how we can update the parameters θ = (κ,σ,τ,α,β) using a Metropolis-Hastings update. First, let us derive the posterior distribution of the hyperparameters.



We write G = G' + ∑_c = 1^K r̃_c δ_ṽ_c where G' is the non observed part. And we note G' the restriction of G' to the locations which intensity is larger than min s.


    p(θ | (s,G), v,Z)    ∝    p( θ | ∑_c = 1^K r̃_c δ_ṽ_c, s, G', Z) 
       ∝    p( θ , ∑_c = 1^K r̃_c δ_ṽ_c, s, G', Z) 
       ∝    p( θ , ∑_c = 1^K r̃_c δ_ṽ_c, Z)   p( s, G' |  θ , ∑_c = 1^K r̃_c δ_ṽ_c, Z)


Now let us derive consider the first part

    p( θ , ∑_c = 1^K r̃_c δ_ṽ_c, Z)    ∝    p(θ)  p(∑_c = 1^K r̃_c δ_ṽ_c | θ) p ( Z| ∑_c = 1^K r̃_c δ_ṽ_c) 
       ∝    p(θ) p(∑_c = 1^K r̃_c δ_ṽ_c | θ) 
       ∝    p(θ) e^-Ψ(n)∏_c=1^K ρ_κ,σ,τ(r̃_c) f_α,β(ṽ_c)


Now let us consider the second part

    p( s, G' |  θ , ∑_c = 1^K r̃_c δ_ṽ_c, Z)    ∝   
    p(s | ∑_c = 1^K r̃_c δ_ṽ_c, Z) p(G' | θ, s ) 
       ∝    p(G' | θ, s ) 
       ∝    e^-Ψ'(min s,n)∏_k ρ_κ,σ,τ(r'_k) f_α,β(v'_k)

where (r'_k) and (v'_k) are respectively the intensities and locations of G'




Let

    π(θ) = e^-Ψ (n) -Ψ'(min s,n)∏_t=1^T μ_θ(r_t,v_t),

where we are taking the product over the T atoms and jumps of G and

    Ψ' (s,n) = ∫_r > s,v e^-r |v|^2 ∏_i=1^n (f(v_i) dv_i) ρ(r) dr.

The posterior satisfies p(θ | (s,G)) ∝ p(θ) π(θ). With our particular choice of distribution of the CRM, the multivariate integrals are reduced to one dimensional integrals, which makes the algorithm tractable. Indeed, we find that


    Ψ (n) + Ψ'(min s,n) = κ/σ∫_0^+∞ (τ+ς^2)^σ [ σ/Γ(1-σ)Γ(-σ,(τ+ς^2)min s ) + 1] f_nα,β(ς)dς - κτ^σ/σ





    Ψ (n)=
    κ/σ∫_0^+∞ (τ+ς^2)^σ  f_nα,β(ς)dς - κτ^σ/σ   ,   if σ > 0  or σ < 0
    κ∫_0^+∞log(τ+ς^2) f_nα,β(ς)dς - κlogτ   ,   if σ =  0

and

    Ψ'(min s,n)=
    κ/Γ (1-σ)∫_0^+∞ (τ+ς^2)^σΓ(-σ,(τ+ς^2)min s ) f_nα,β(ς)dς,      if σ > 0  or σ < 0
    κ∫_0^+∞∫_(τ+ς^2)min s^+∞e^-r/r dr f_nα,β(ς)dς   ,   if σ =  0




We use the following priors:

    1-2σ   ∼    Gamma(a_σ,b_σ) 
    κ   ∼    Gamma(a_κ,b_κ) 
    τ   ∼    Gamma(a_τ,b_τ) 
    α   ∼    Gamma(a_α,b_α) 
    β   ∼    Gamma(a_β,b_β)

And proposals

    1-2σ̃ | σ   ∼     Lognormal(log (1-2σ), Σ_σ) 
    κ̃ | κ   ∼    Lognormal(logκ, Σ_κ) 
    τ̃|τ   ∼    Lognormal(logτ, Σ_τ) 
    α̃ | α   ∼    Lognormal(logα, Σ_α) 
    β̃|β   ∼    Lognormal(logβ, Σ_β)

We find that

    logp(θ̃) q(θ|θ̃)/p(θ) q(θ̃ | θ)   =    a_σlog1-2σ̃/1-2σ + 2b_σ (σ̃-σ) 
       +    a_κlogκ̃/κ - b_κ (κ̃-κ) 
       +    a_τlogτ̃/τ - b_τ (τ̃-τ) 
       +    a_αlogα̃/α - b_α (α̃-α) 
       +    a_βlogβ̃/β - b_β (β̃-β)

And

    logπ(θ̃)/π(θ)   =   Ψ_θ(n) + Ψ'_θ(min s, n) - Ψ_θ̃(n) - Ψ'_θ̃(min s, n) 
       +    T logκ̃/κ - T logΓ(1-σ̃)/Γ(1-σ) -n T logΓ(α̃)/Γ(α)
       -    (σ̃ - σ) ∑_t log r_t - (τ̃ - τ) ∑_t r_t 
       +    (α̃-α) ∑_t,ilog v_t,i - (β̃ - β) ∑_t,i v_t,i
       +    n T (α̃logβ̃ - αlogβ )






 §.§ Sampling from the inhomogeneous CRM

In this section we show how we can sample from the inhomogeneous CRM G' with measure:

    μ'(dr,dv) = e^-r (∑_i=1^n v_i)^2ρ_κ,τ,σ(r) [∏_i=1^n f_α,β(v_i)dv_i] dr

Let us recall that f_α,β is the gamma pdf and ρ_κ,τ,σ the GGP intensity. From Section 4, we know that if we make the following change of variables
(v_1,..,v_n) ↦ (ς =∑_i v_i, ν_1 = v_1/s,...,ν_n = v_n/s), we get

    μ'(dr,dς,dν) = e^-r ς^2ρ_κ,τ,σ(r) f_nα,β(ς)dς dr Dir(dν_1,...,dν_n; α)

Hence, we can sample independently (r,ς) and ν. From one hand, (ν_1,..,ν_n) is sampled from a Dirichlet distribution with parameter α = (α,..,α).
On the other hand, the total sum ς and the intensity r are sampled from

    μ(r,ς)    =     e^-r ς^2ρ_κ,τ,σ(r) f_nα,β(ς) 
       =   κ/Γ (1-σ) e^-r (ς^2+τ) r^-1-σ f_nα,β(ς)

Now, to reduce the problem to sampling from a homogeneous CRM, let us consider the change of variable (r,ς) ↦ (r= r[τ+ς^2], s) which determinant is τ+ς^2. We find finally that

    μ(r,ς) = κ/Γ (1-σ)  e^-r r^-1-σ (ς^2+τ)^σ f_nα,β(ς)

Besides, since r≥ r τ  ∀ (i,j), we only need to sample the points such that r≥τmin s. Therefore, since τ > 0, we sample a finite number of atoms. Then, we only keep the points such as r > min s.
Finally, let us notice that in our setting, even with σ≤ 0

    ∫_ς (ς^2+τ)^σf_nα,β(ς) dς = σ/κΨ(n) + τ^σ

Therefore, we first sample the jumps from the levy measure

    ρ(r) =  σΨ(n) + κτ^σ/Γ (1-σ)  e^-r r^-1-σ1_r≥τmin s

using adaptive thinning <cit.>. Then, we sample ς with pdf ∝ (ς^2+τ)^σf_nα,β(ς) using rejection sampling.



