# A Bayesian Inference Method Using Monte Carlo Sampling for Estimating the Number of Communities in Bipartite Networks

Authors: Guo-Zheng Wang; Li Xiong; Hu-Chen Liu
Journal: Scientific Programming
 (2019)

---

## Abstract

Community detection is an important analysis task for complex networks, including bipartite networks, which consist of nodes of two types and edges connecting only nodes of different types. Many community detection methods take the number of communities in the networks as a fixed known quantity; however, it is impossible to give such information in advance in real-world networks. In our paper, we propose a projection-free Bayesian inference method to determine the number of pure-type communities in bipartite networks. This paper makes the following contributions: (1) we present the first principle derivation of a practical method, using the degree-corrected bipartite stochastic block model that is able to deal with networks with broad degree distributions, for estimating the number of pure-type communities of bipartite networks; (2) a prior probability distribution is proposed over the partition of a bipartite network; (3) we design a Monte Carlo algorithm incorporated with our proposed method and prior probability distribution. We give a demonstration of our algorithm on synthetic bipartite networks including an easy case with a homogeneous degree distribution and a difficult case with a heterogeneous degree distribution. The results show that the algorithm gives the correct number of communities of synthetic networks in most cases and outperforms the projection method especially in the networks with heterogeneous degree distributions.


---

## Body

## 1. Introduction

A bipartite network is a network with nodes of two types and edges connecting only nodes of different types. The decomposition of bipartite networks into communities (clusters, modules, or groups), i.e., community detection, plays an important role in revealing the structure of large networked systems, providing new insights into how the network is organized [1–4].Many methods [5–9] have been developed for community detection in bipartite networks in recent years. A fundamental shortcoming of most community detection methods is that they partition networks into a fixed number of groups. However, this number is usually unknown in real-world networks, and we need to mine such information from the network data. A lot of research [10–13] for making such efforts to determine the number of communities in bipartite networks has been proposed recently. There are three main problems in these methods. One is that they performed estimation through maximizing the modularity proposed in [10, 14] that is proved to be NP-hard [15, 16]; the second is that they gave the number of communities of mixed-type, which is nearly always substantially less efficient [6]; the third is that the projection method [15] performed poorly due to information loss. The heuristic methods proposed in [17, 18] for community detection in bipartite networks does not need the number of communities to be given a priori.In this paper, we propose a projection-free Bayesian inference method for determining the number of pure-type communities in a bipartite network. Our method builds mainly on the work as below: (i) the degree-corrected bipartite stochastic block model, proposed by Larremore et al. [6], is used to find the community structure of empirical networks with broad degree distributions; and (ii) the prior probability distribution over divisions of a network into groups and a new prior probability distribution based on a queueing-type process, both proposed by Riolo et al. [4], are used for calculating the number of communities in a unipartite network.In Section2, first, we present the first-principle derivation of a practical method, using the degree-corrected bipartite stochastic block model that is able to deal with networks with broad degree distributions, for estimating the number of pure-type communities of bipartite networks. Second, we propose a prior probability distribution over the partition of a bipartite network, with the community-type parameter ensuring that each community is pure type. In Section 3, we design a Monte Carlo algorithm incorporated with our proposed method and prior probability distribution. In the following section, we demonstrate our method on synthetic bipartite networks including an easy case with homogeneous degree distributions and a difficult case with a heterogeneous degree distribution. The results show that the proposed algorithm can determine the correct number of communities and perform better than our projection method in every case.


## 2. Methods

### 2.1. Degree-Corrected Bipartite Stochastic Block Model

The stochastic block model is a generative model used to produce networks containing blocks, groups, or communities. This model is very important in network science and is used for recovering the community structure in network data [2, 3]. The classic stochastic block model can be described as follows: divide the number of vertices N into K disjoint communities; any two vertices i and j are connected by an edge with probability ωgigj, which is an entry of a symmetric K × K matrix, and gi is the community of vertex i. However, the block model described above finds the community structure merely due to the degree sequence and fails to detect the known communities in a real-world network that has heterogeneous degree distributions [19]. Karrer and Newman [20] extended the classic stochastic block model including heterogeneity in the degrees of vertices and proposed the degree-corrected stochastic block model, which is proved to overcome the problems of the classic block model.Most stochastic block model community detection methods can be naturally applied to bipartite networks [20, 21]. Unfortunately, the stochastic block model often overfits bipartite data by mixing nodes of different types within communities and it is nearly always substantially less efficient [6]. Built on the work of Karrer and Newman [20], Larremore et al. [6] proposed the degree-corrected bipartite stochastic block model, which is employed in our calculations. In the degree-corrected bipartite stochastic block model, a bipartite network G is given with an Na×Nb bipartite asymmetric adjacency matrix B, where Na is the number of nodes of type-a and Nb is the number of nodes of type-b. Let A=0BBT0 be the N×N symmetric adjacency matrix of the network G with N=Na+Nb. The type-a nodes are divided into some number ka of communities, labeled 1,…,ka, and the Nb nodes of type-b are divided into kb communities, labeled ka+1,…,ka+kb. We express the matrix of community interrelationships as a k×k matrix, where k=ka+kb. Let gi again encode the community node i belongs to. Let ti be the type of vertex i and Tr be the type of community r, imposing the constraint(1)ti=Tgi,which indicates that node types and community types must match and ensure that communities will be pure type. We write(2)ka=∑r=1kδTr,type‐a,(3)kb=∑r=1kδTr,type‐b,where δ is the Kronecker delta. Let θi control the expected degree of node i and ωrs be the k×k symmetric matrix of parameters to control the number of edges between communities r and s. Following [4], the normalization of θi can be fixed by imposing the constraint(4)1n∑iθiδgi,r=1,where nr=∑iδgi,r is the number of nodes in community r. Following [22], we let the numbers of edges between nodes i and j follow a Poisson distribution with mean θiθjωgigj. Enforcing the bipartite constraint of equation (1) produces a restriction on ω:(5)ωrs=0,whenTr=TS.Given parametersg, k, θ, ω, and T for the specification of the mode, the probability of observing a bipartite network G with adjacency matrix A can be written as(6)PAg,k,θ,ω,T=∏i<jti≠tjθiθjωgigjAijAij!exp−θiθjωgigj.Allowing for the constraint of equation (4), the probability PAg,k,θ,ω,T can be simplified to the more convenient form of(7)PAg,k,θ,ω,T=∏iθidi×∏r<sTr≠Tsωrsmrsexp−nrnsωrs,where di is the observed degree of vertexes i and mrs=∑i,jti≠tjAijδgi,rδgj,s is the number of edges between communities r and s. We have neglected an overall multiplicative constant in (7) since it cancels out in later calculations. Note that a similar probability given by Larremore et al. [6] has been modified in equation (7) as follows:(i)
The numberk of communities, the objective we will estimate, is incorporated as an unknown quantity(ii)
The exponential expression is−nrnsωrs rather than −ωrs, with the normalization of θi under a different constraint conditionThen, we integrate out the irrelevant parametersθ and ω. We assume maximum-entropy (i.e., least informative) prior probability distributions on the parameters θ and ω. For θ, this means a uniform prior probability distribution over the regular simplex of values specified by equation (4). Then, we let the expected value of the edge probability ωrs be equal to the observed average edge probability in the network as a whole: p=2m/N2, where m is the total number of edges in the bipartite network. Then, the maximum-entropy prior probability distribution is an exponential distribution Pω=1/pe−ω/p. We assume the priors to be independent (conditioned on g, k, and T) so that Pω,θ∣k,g,T=Pω∣k,TPθ∣k,g,T and(8)PAg,k,T=∬PAg,k,θ,ω,TPθg,k,TPωk,Tdθdω.With these choices of priors, integration is performed on equation (8). Then, we have(9)PAg,k,T=∏rnrκrnr−1!nr+κr−1!×∏r<sTr≠Tsmrs!pnrns+1mrs+1,where κr=∑idiδgi,r and an overall multiplying constant has been discarded.


### 2.2. Prior on Community Partitions

Our goal is to estimate the correct values ofka and kb for a given bipartite network using this model as the basis for a Bayesian model selection procedure. We have(10)Pg,k,TA=Pg,k,TPAg,k,TPA,where PAg,k,T is given by equation (6), and the probability P(A), which in the denominator of equation (10) has no effect on our results, is unknown but cancels out in later calculations. In this paper, our primary focus is to get the posterior distribution on k through summing over g; then, we choose a value for k and calculate ka and kb using equations (2) and (3). Now, we start to choose the prior Pg,k,T, which is often the most important and difficult task of the calculation in the case of Bayesian methods.


#### 2.2.1. PriorPgk on Community Partitions

If we know the number of communitiesk in advanced, let us choose the prior Pgk on community partitions of one type of node. We first employ the most commonly used approach, which is described as follows. The prior on the community partition probabilities γ is uniform under the constraint ∑rγr=1, where γr∈0,1, with which nodes are assigned to communities independently at random. We can get a particular community partition with the probability(11)Pgγ,k=∏i=1Nγgi=∏r=1kγrnr.The valuesγr fall on a regular k−1-dimensional simplex with volume k−1!, so its probability density is Pγk=k−1!. We integrate equation (11) over the simplex and get the following equation [4, 23, 24]:(12)Pgk=∫Pgγ,kPγkdγ=k−1!N+k−1!∏rnr!.Since the process above generates a uniform distribution over possible community sizes, we then introduce an alternative and simpler way (used by Riolo et al. [4]) to derive the prior Pgk. We have N+k−1k−1 possible ways to choose k communities with ∑r=1knr=N and N!/∏rnr! possible ways to place the nodes in the k communities; thus, any partition g of nodes to communities is given with the probability(13)Pgk=1N−1k−1N!/∏rnr!=k−1!N+k−1!∏rnr!,the same to equation (12) without the need for parameters γr.However, these two methods may generate partitionsg with empty communities. As in [4], we have the binomial coefficient N−1k−1 possible choices of k communities with nonempty ones and(14)PgK=1N−1k−1N!/∏rnr!.Then, we allow for two different types of partitions and get(15)Pgk,T=1Na−1ka−1Na!/∏Tr=type−anr!1Nb−1kb−1Na!/∏Tr=type−bnr!=1Na−1ka−1Nb−1kb−1Na!Nb!/∏rnr!.


#### 2.2.2. Choice of the Number of Communities

Some previous work has been done for the choices of priorP(k) over the number of communities itself, such as letting P(k) equal 1/N [23, 24] or 1/k! [21]. We again follow [4] and take a different approach, in which community partitions g and the number of communities k can be generated synchronously. We use a queueing-type mechanism for processing community partitions g. For one type, such as type-a, we order the Na nodes uniformly at random and the first node is placed in community 1. Then, we place each following node either (a) with probability 1−q in the same community as the previous node or (b) with probability q in the next community; for another type, we repeat the process above. This process ensures that all communities generated are not empty.For type-a, there are Na! (Nb! for type-b) possible ways to order the nodes, so the probability of each one occurring is the same as 1/Na! (1/Nb! for type-b). If k=ka+kb communities are generated finally, we must create k − 2 new communities (ka−1 new type-a ones and kb−1 new type-b ones). Because for one type each node except the first starts a new community with equal probability q, k communities with sizes n1,…,nk are generated with the probability:(16)1−qn1−1q1−qn2−1q⋯1−qnka−11−qnka+1−1q1−qnka+2−1q⋯1−qnka+kb−1=qk−21−qN−k.For each community of the same partitiong, the nodes can be rearranged in ∏rnr! ways. Thus, any given partition is generated in the process with the probability:(17)Pg,k,T=1Na!Nb!qk−21−qn−k∏r=1knr!.Given thatPg,k,T=Pk,TPgk,T,(18)Pk,T=Pg,k,TPgk=Na−1ka−1Nb−1kb−1qk−21−qN−k.We letq=μ/N−1, where the expected number of new communities created is µ. Then,(19)Pg,k,T=1−qNq2Na!Nb!μkN−μ−1k∏r=1knr!.In equation (19), 1−qN/q2Na!Nb! has no effect on our result and cancels out in later calculations.As in [4], we let µ = 1 and neglecting constants(20)Pg,k,T=N−2−k∏r=1knr!.Now, equation (10) can be written as(21)Pg,k,TA=N−2−k∏rnrκrnr!nr−1!nr+κr−1!×∏r<sTr≠Tsmrs!pnrns+1mrs+1,and here we allow for equations (9) and (20).Unfortunately, it is hard to sum overg since the sum has kN terms [4]. Instead, we approximate the distribution over k (ka and kb according to different types) by Markov chain Monte Carlo sampling.


## 2.1. Degree-Corrected Bipartite Stochastic Block Model

The stochastic block model is a generative model used to produce networks containing blocks, groups, or communities. This model is very important in network science and is used for recovering the community structure in network data [2, 3]. The classic stochastic block model can be described as follows: divide the number of vertices N into K disjoint communities; any two vertices i and j are connected by an edge with probability ωgigj, which is an entry of a symmetric K × K matrix, and gi is the community of vertex i. However, the block model described above finds the community structure merely due to the degree sequence and fails to detect the known communities in a real-world network that has heterogeneous degree distributions [19]. Karrer and Newman [20] extended the classic stochastic block model including heterogeneity in the degrees of vertices and proposed the degree-corrected stochastic block model, which is proved to overcome the problems of the classic block model.Most stochastic block model community detection methods can be naturally applied to bipartite networks [20, 21]. Unfortunately, the stochastic block model often overfits bipartite data by mixing nodes of different types within communities and it is nearly always substantially less efficient [6]. Built on the work of Karrer and Newman [20], Larremore et al. [6] proposed the degree-corrected bipartite stochastic block model, which is employed in our calculations. In the degree-corrected bipartite stochastic block model, a bipartite network G is given with an Na×Nb bipartite asymmetric adjacency matrix B, where Na is the number of nodes of type-a and Nb is the number of nodes of type-b. Let A=0BBT0 be the N×N symmetric adjacency matrix of the network G with N=Na+Nb. The type-a nodes are divided into some number ka of communities, labeled 1,…,ka, and the Nb nodes of type-b are divided into kb communities, labeled ka+1,…,ka+kb. We express the matrix of community interrelationships as a k×k matrix, where k=ka+kb. Let gi again encode the community node i belongs to. Let ti be the type of vertex i and Tr be the type of community r, imposing the constraint(1)ti=Tgi,which indicates that node types and community types must match and ensure that communities will be pure type. We write(2)ka=∑r=1kδTr,type‐a,(3)kb=∑r=1kδTr,type‐b,where δ is the Kronecker delta. Let θi control the expected degree of node i and ωrs be the k×k symmetric matrix of parameters to control the number of edges between communities r and s. Following [4], the normalization of θi can be fixed by imposing the constraint(4)1n∑iθiδgi,r=1,where nr=∑iδgi,r is the number of nodes in community r. Following [22], we let the numbers of edges between nodes i and j follow a Poisson distribution with mean θiθjωgigj. Enforcing the bipartite constraint of equation (1) produces a restriction on ω:(5)ωrs=0,whenTr=TS.Given parametersg, k, θ, ω, and T for the specification of the mode, the probability of observing a bipartite network G with adjacency matrix A can be written as(6)PAg,k,θ,ω,T=∏i<jti≠tjθiθjωgigjAijAij!exp−θiθjωgigj.Allowing for the constraint of equation (4), the probability PAg,k,θ,ω,T can be simplified to the more convenient form of(7)PAg,k,θ,ω,T=∏iθidi×∏r<sTr≠Tsωrsmrsexp−nrnsωrs,where di is the observed degree of vertexes i and mrs=∑i,jti≠tjAijδgi,rδgj,s is the number of edges between communities r and s. We have neglected an overall multiplicative constant in (7) since it cancels out in later calculations. Note that a similar probability given by Larremore et al. [6] has been modified in equation (7) as follows:(i)
The numberk of communities, the objective we will estimate, is incorporated as an unknown quantity(ii)
The exponential expression is−nrnsωrs rather than −ωrs, with the normalization of θi under a different constraint conditionThen, we integrate out the irrelevant parametersθ and ω. We assume maximum-entropy (i.e., least informative) prior probability distributions on the parameters θ and ω. For θ, this means a uniform prior probability distribution over the regular simplex of values specified by equation (4). Then, we let the expected value of the edge probability ωrs be equal to the observed average edge probability in the network as a whole: p=2m/N2, where m is the total number of edges in the bipartite network. Then, the maximum-entropy prior probability distribution is an exponential distribution Pω=1/pe−ω/p. We assume the priors to be independent (conditioned on g, k, and T) so that Pω,θ∣k,g,T=Pω∣k,TPθ∣k,g,T and(8)PAg,k,T=∬PAg,k,θ,ω,TPθg,k,TPωk,Tdθdω.With these choices of priors, integration is performed on equation (8). Then, we have(9)PAg,k,T=∏rnrκrnr−1!nr+κr−1!×∏r<sTr≠Tsmrs!pnrns+1mrs+1,where κr=∑idiδgi,r and an overall multiplying constant has been discarded.


## 2.2. Prior on Community Partitions

Our goal is to estimate the correct values ofka and kb for a given bipartite network using this model as the basis for a Bayesian model selection procedure. We have(10)Pg,k,TA=Pg,k,TPAg,k,TPA,where PAg,k,T is given by equation (6), and the probability P(A), which in the denominator of equation (10) has no effect on our results, is unknown but cancels out in later calculations. In this paper, our primary focus is to get the posterior distribution on k through summing over g; then, we choose a value for k and calculate ka and kb using equations (2) and (3). Now, we start to choose the prior Pg,k,T, which is often the most important and difficult task of the calculation in the case of Bayesian methods.


### 2.2.1. PriorPgk on Community Partitions

If we know the number of communitiesk in advanced, let us choose the prior Pgk on community partitions of one type of node. We first employ the most commonly used approach, which is described as follows. The prior on the community partition probabilities γ is uniform under the constraint ∑rγr=1, where γr∈0,1, with which nodes are assigned to communities independently at random. We can get a particular community partition with the probability(11)Pgγ,k=∏i=1Nγgi=∏r=1kγrnr.The valuesγr fall on a regular k−1-dimensional simplex with volume k−1!, so its probability density is Pγk=k−1!. We integrate equation (11) over the simplex and get the following equation [4, 23, 24]:(12)Pgk=∫Pgγ,kPγkdγ=k−1!N+k−1!∏rnr!.Since the process above generates a uniform distribution over possible community sizes, we then introduce an alternative and simpler way (used by Riolo et al. [4]) to derive the prior Pgk. We have N+k−1k−1 possible ways to choose k communities with ∑r=1knr=N and N!/∏rnr! possible ways to place the nodes in the k communities; thus, any partition g of nodes to communities is given with the probability(13)Pgk=1N−1k−1N!/∏rnr!=k−1!N+k−1!∏rnr!,the same to equation (12) without the need for parameters γr.However, these two methods may generate partitionsg with empty communities. As in [4], we have the binomial coefficient N−1k−1 possible choices of k communities with nonempty ones and(14)PgK=1N−1k−1N!/∏rnr!.Then, we allow for two different types of partitions and get(15)Pgk,T=1Na−1ka−1Na!/∏Tr=type−anr!1Nb−1kb−1Na!/∏Tr=type−bnr!=1Na−1ka−1Nb−1kb−1Na!Nb!/∏rnr!.


### 2.2.2. Choice of the Number of Communities

Some previous work has been done for the choices of priorP(k) over the number of communities itself, such as letting P(k) equal 1/N [23, 24] or 1/k! [21]. We again follow [4] and take a different approach, in which community partitions g and the number of communities k can be generated synchronously. We use a queueing-type mechanism for processing community partitions g. For one type, such as type-a, we order the Na nodes uniformly at random and the first node is placed in community 1. Then, we place each following node either (a) with probability 1−q in the same community as the previous node or (b) with probability q in the next community; for another type, we repeat the process above. This process ensures that all communities generated are not empty.For type-a, there are Na! (Nb! for type-b) possible ways to order the nodes, so the probability of each one occurring is the same as 1/Na! (1/Nb! for type-b). If k=ka+kb communities are generated finally, we must create k − 2 new communities (ka−1 new type-a ones and kb−1 new type-b ones). Because for one type each node except the first starts a new community with equal probability q, k communities with sizes n1,…,nk are generated with the probability:(16)1−qn1−1q1−qn2−1q⋯1−qnka−11−qnka+1−1q1−qnka+2−1q⋯1−qnka+kb−1=qk−21−qN−k.For each community of the same partitiong, the nodes can be rearranged in ∏rnr! ways. Thus, any given partition is generated in the process with the probability:(17)Pg,k,T=1Na!Nb!qk−21−qn−k∏r=1knr!.Given thatPg,k,T=Pk,TPgk,T,(18)Pk,T=Pg,k,TPgk=Na−1ka−1Nb−1kb−1qk−21−qN−k.We letq=μ/N−1, where the expected number of new communities created is µ. Then,(19)Pg,k,T=1−qNq2Na!Nb!μkN−μ−1k∏r=1knr!.In equation (19), 1−qN/q2Na!Nb! has no effect on our result and cancels out in later calculations.As in [4], we let µ = 1 and neglecting constants(20)Pg,k,T=N−2−k∏r=1knr!.Now, equation (10) can be written as(21)Pg,k,TA=N−2−k∏rnrκrnr!nr−1!nr+κr−1!×∏r<sTr≠Tsmrs!pnrns+1mrs+1,and here we allow for equations (9) and (20).Unfortunately, it is hard to sum overg since the sum has kN terms [4]. Instead, we approximate the distribution over k (ka and kb according to different types) by Markov chain Monte Carlo sampling.


## 2.2.1. PriorPgk on Community Partitions

If we know the number of communitiesk in advanced, let us choose the prior Pgk on community partitions of one type of node. We first employ the most commonly used approach, which is described as follows. The prior on the community partition probabilities γ is uniform under the constraint ∑rγr=1, where γr∈0,1, with which nodes are assigned to communities independently at random. We can get a particular community partition with the probability(11)Pgγ,k=∏i=1Nγgi=∏r=1kγrnr.The valuesγr fall on a regular k−1-dimensional simplex with volume k−1!, so its probability density is Pγk=k−1!. We integrate equation (11) over the simplex and get the following equation [4, 23, 24]:(12)Pgk=∫Pgγ,kPγkdγ=k−1!N+k−1!∏rnr!.Since the process above generates a uniform distribution over possible community sizes, we then introduce an alternative and simpler way (used by Riolo et al. [4]) to derive the prior Pgk. We have N+k−1k−1 possible ways to choose k communities with ∑r=1knr=N and N!/∏rnr! possible ways to place the nodes in the k communities; thus, any partition g of nodes to communities is given with the probability(13)Pgk=1N−1k−1N!/∏rnr!=k−1!N+k−1!∏rnr!,the same to equation (12) without the need for parameters γr.However, these two methods may generate partitionsg with empty communities. As in [4], we have the binomial coefficient N−1k−1 possible choices of k communities with nonempty ones and(14)PgK=1N−1k−1N!/∏rnr!.Then, we allow for two different types of partitions and get(15)Pgk,T=1Na−1ka−1Na!/∏Tr=type−anr!1Nb−1kb−1Na!/∏Tr=type−bnr!=1Na−1ka−1Nb−1kb−1Na!Nb!/∏rnr!.


## 2.2.2. Choice of the Number of Communities

Some previous work has been done for the choices of priorP(k) over the number of communities itself, such as letting P(k) equal 1/N [23, 24] or 1/k! [21]. We again follow [4] and take a different approach, in which community partitions g and the number of communities k can be generated synchronously. We use a queueing-type mechanism for processing community partitions g. For one type, such as type-a, we order the Na nodes uniformly at random and the first node is placed in community 1. Then, we place each following node either (a) with probability 1−q in the same community as the previous node or (b) with probability q in the next community; for another type, we repeat the process above. This process ensures that all communities generated are not empty.For type-a, there are Na! (Nb! for type-b) possible ways to order the nodes, so the probability of each one occurring is the same as 1/Na! (1/Nb! for type-b). If k=ka+kb communities are generated finally, we must create k − 2 new communities (ka−1 new type-a ones and kb−1 new type-b ones). Because for one type each node except the first starts a new community with equal probability q, k communities with sizes n1,…,nk are generated with the probability:(16)1−qn1−1q1−qn2−1q⋯1−qnka−11−qnka+1−1q1−qnka+2−1q⋯1−qnka+kb−1=qk−21−qN−k.For each community of the same partitiong, the nodes can be rearranged in ∏rnr! ways. Thus, any given partition is generated in the process with the probability:(17)Pg,k,T=1Na!Nb!qk−21−qn−k∏r=1knr!.Given thatPg,k,T=Pk,TPgk,T,(18)Pk,T=Pg,k,TPgk=Na−1ka−1Nb−1kb−1qk−21−qN−k.We letq=μ/N−1, where the expected number of new communities created is µ. Then,(19)Pg,k,T=1−qNq2Na!Nb!μkN−μ−1k∏r=1knr!.In equation (19), 1−qN/q2Na!Nb! has no effect on our result and cancels out in later calculations.As in [4], we let µ = 1 and neglecting constants(20)Pg,k,T=N−2−k∏r=1knr!.Now, equation (10) can be written as(21)Pg,k,TA=N−2−k∏rnrκrnr!nr−1!nr+κr−1!×∏r<sTr≠Tsmrs!pnrns+1mrs+1,and here we allow for equations (9) and (20).Unfortunately, it is hard to sum overg since the sum has kN terms [4]. Instead, we approximate the distribution over k (ka and kb according to different types) by Markov chain Monte Carlo sampling.


## 3. Monte Carlo Algorithm for Bipartite Networks

### 3.1. Our Algorithm

We design a Monte Carlo algorithm incorporated with the bipartite block model and prior probability distribution discussed above to apply the bipartite networks. We call our algorithm the bipartite network Monte Carlo algorithm (BMCA), and it is built on the unipartite network analysis of Riolo et al. [4]. Our algorithm fulfills the requirements of ergodicity and detailed balance [25].There are two types of steps used by BMCA:Type 1: moving one node from its current community to a different existing community. There are again two types of processes in this type of rearrangement. In the processes of the first type, BMCA decreases the number of communities (ka or kb) of the same type as the community whose last node moved, thereby decreasing the value of k by one. In the processes of the second type, the community the node moved from contains more than one node and the move here does not change the number of communities.Type 2: moving one node to a community newly created. The number of communities (ka or kb) of the same type as the community the node moved from and the value of k increase by one.The two types of steps described above make BMCA meet the requirement of ergodicity.Detailed balance requires that the rateRg,k,T⟶g′,k′,T′ goes from a state g,k,T to another state g′,k′,T′ and the opposite meet(22)Rg,k,T⟶g′,k′,T′Rg′,k′,T′⟶g,k,T=Pg′,k′,T′APg,k,TA=Pg′,k′,T′Pg,k,T×PAg′,k′,T′PAg,k,T,where we allow for equation (10). From (20), we have(23)Pg′,k′,T′Pg,k,T=N−2k−k′∏r=1knr′!∏r=1knr!.We considerRg,k,T⟶g′,k′,T′ as πg,k,T⟶g′,k′,T′αg,k,T⟶g′,k′,T′, where the previous part of the product represents the probability of proposing a move and the latter represents the probability chosen to satisfy the detailed balance condition for accepting the move. Then,(24)Rg,k,T⟶g′,k′,T′Rg′,k′,T′⟶g,k,T=πg,k,T⟶g′,k′,T′πg′,k′,T′⟶g,k,T×αg,k,T⟶g′,k′,T′αg′,k′,T′⟶g,k,T.BMCA is described as follows:Input: the bipartite adjacency matrix B and the node-type vector ti.Initial communities partition: using the process described in Section 2.2.2.Monte Carlo Sampling:(1)
(a) In each step of BMCA, we carry out a rearrangement of type 1 with probability1−1/N−1. If k=2 (i.e., ka=1 and kb=1), we do nothing. Otherwise, when k>2, first, we randomly select a community label r in the range 1,…,k. If the number of communities of type Tr is more than one, then we randomly select a community of type Tr labels s; otherwise, we turn to communities of another type, from which we randomly reselect a pair of communities, respectively, labels r and s. Then, we randomly select one node from community r and move it to community s. The number k=ka+kb of total communities remains constant.(b) In the process, if communityr becomes empty as a result that its last node is removed, the number of communities k decreases by one. In practice, for the type-a communities labels 1,…,ka, we can efficiently change the community ka to have label r and then change the community k to have label ka; specially, if r=ka, we only perform the latter relabeling. In such a process, the number ka of communities of type-a decreases by 1 and the number kb of communities of type-b remains constant. For type-b community labels ka+1,…,k, we can efficiently change the community k to have label r; specifically, if r=k, no relabeling is necessary. In such a process, the number kb of communities of type-b decreases by 1 and the number ka of communities of type-a remains constant. The number ka=ka+kb of total communities decreases by 1.(2)
Otherwise, we carry out a rearrangement of type 2 with probability1/N−1. We randomly select a community label r in the range 1,…,k. If there is only one node in community r, we do nothing. Otherwise, if Tr=type−a, we change community label ka+1 to k+1 and create a new empty community ka+1. Then, we randomly select a node from community r and move it to the newly created community ka+1. During this process, the number ka of communities of type-a increases by 1 and the number kb of communities of type-b remains constant. If Tr=type−b, we simply create a new empty community k+1 and no relabeling is necessary. Then, we randomly select a node from community r and move it to the newly created community k+1; during this process, the number kb of communities of type-b increases by 1 and the number ka of communities of type-a remains constant. The number k=ka+kb of total communities increases by 1.(3)
We accept the rearrangement proposed above with acceptance probability [4]:(25)αg,k,T⟶g′,k′,T′=min1,PAg′,k′,T′PAg,k,T.(4)
Repeat steps 1–3.Output: the posterior probabilitiesPkaA and PkbA.Rearrangements of any type are performed, and we always have the following equation (see Table1):(26)πg,k,T⟶g′,k′,T′πg′,k′,T′⟶g,k,T=Pg′,k′,T′Pg,k,T.Table 1
Derivation of the relation ofπg,k,T⟶g′,k′,T′/πg′,k′,T′⟶g,k,T and Pg′,k′,T′/Pg,k,T.
Type of rearrangement
States
Probability of performing the move
Probability of selecting two distinct communities of the same type
Probability of selecting a specific node from or to distinct communities of the same type
Total probabilityπ
π
g
,
k
,
T
⟶
g
′
,
k
′
,
T
′
/
π
g
′
,
k
′
,
T
′
⟶
g
,
k
,
T
P
g
′
,
k
′
,
T
′
/
P
g
,
k
,
T (using equation (23))
Type 1 withk′=k
g
,
k
,
T
⟶
g
′
,
k
′
,
T
′
1 − 1/(N − 1)
1
/
k
a
k
a
−
1
+
k
b
k
b
−
1
1
/
n
r
1
−
1
/
N
−
1
×
1
/
k
a
k
a
−
1
+
k
b
k
b
−
1
×
1
/
n
r
n
s
′
/
n
r
n
s
′
/
n
r
g
′
,
k
′
,
T
′
⟶
g
,
k
,
T
1 − 1/(N − 1)
1
/
k
a
k
a
−
1
+
k
b
k
b
−
1
1
/
n
s
′
1
−
1
/
N
−
1
×
1
/
k
a
k
a
−
1
+
k
b
k
b
−
1
×
1
/
n
s
′
Type 1 withk′=k−1
g
,
k
,
T
⟶
g
′
,
k
′
,
T
′
1 − 1/(N − 1)
1
/
k
a
k
a
−
1
+
k
b
k
b
−
1
1 (nownr = 1)
1
−
1
/
N
−
1
×
1
/
k
a
k
a
−
1
+
k
b
k
b
−
1
N
−
2
n
s
′
N
−
2
n
s
′
g
′
,
k
′
,
T
′
⟶
g
,
k
,
T
1/(N − 1)
1
/
k
a
k
a
−
1
+
k
b
k
b
−
1
1
/
n
s
′
1
/
N
−
1
×
1
/
k
a
k
a
−
1
+
k
b
k
b
−
1
×
1
/
n
s
′
Type 2 withk′=k+1
g
,
k
,
T
⟶
g
′
,
k
′
,
T
′
1/(N − 1)
1
/
k
a
+
1
k
a
+
k
b
+
1
k
b
1
/
n
r
1
/
N
−
1
×
1
/
k
a
+
1
k
a
+
k
b
+
1
k
b
×
1
/
n
r
1
/
N
−
2
n
r
1
/
N
−
2
n
r
g
′
,
k
′
,
T
′
⟶
g
,
k
,
T
1 − 1/(N − 1)
1
/
k
a
+
1
k
a
+
k
b
+
1
k
b
1 (nowns′ = 1)
1
−
1
/
N
−
1
×
1
/
k
a
+
1
k
a
+
k
b
+
1
k
bTaking into consideration equations (22) and (24), the detailed balance condition can be written as(27)αg,k,T⟶g′,k′,T′αg′,k′,T′⟶g,k,T=PAg′,k′,T′PAg,k,T.Therefore, our algorithm satisfies the detailed balance with acceptance probability of equation (25) and will sample correctly from the distribution Pg,k,TA.


### 3.2. Output of BMCA

In our implementation, a given number of stepsti per node is performed in a Monte Carlo run on a bipartite network, and then we write approximate posterior probabilities:(28)Pk,TA=the timeskshowed out  in a Monte Carlo runS.For type-a, we can get approximately(29)PkaA=the timeskashowed out in a Monte Carlo runS,and for type-b,(30)PkbA=the timeskbshowed out in a Monte Carlo runS.The most likely number of type-a communities in a bipartite network is ka with the biggest value PkaA, and the number of type-b communities kb can be given in the same way.In order to avoid any bias in the results and improve the correctness of BMCA, instead of just usingPkaA, we performed a given number M of Monte Carlo runs for a bipartite network. Thus, we obtained the average value of each PikaA as the final posterior probabilities:(31)P′kaA=∑i=1MPikaAM,for type-a and(32)P′kbA=∑i=1MPikbAM,for type-b.We takeON+Na∗Nb to calculate the acceptance probability of each node move, so the time complexity of our algorithm is ONa2Nb+NaNb2.


## 3.1. Our Algorithm

We design a Monte Carlo algorithm incorporated with the bipartite block model and prior probability distribution discussed above to apply the bipartite networks. We call our algorithm the bipartite network Monte Carlo algorithm (BMCA), and it is built on the unipartite network analysis of Riolo et al. [4]. Our algorithm fulfills the requirements of ergodicity and detailed balance [25].There are two types of steps used by BMCA:Type 1: moving one node from its current community to a different existing community. There are again two types of processes in this type of rearrangement. In the processes of the first type, BMCA decreases the number of communities (ka or kb) of the same type as the community whose last node moved, thereby decreasing the value of k by one. In the processes of the second type, the community the node moved from contains more than one node and the move here does not change the number of communities.Type 2: moving one node to a community newly created. The number of communities (ka or kb) of the same type as the community the node moved from and the value of k increase by one.The two types of steps described above make BMCA meet the requirement of ergodicity.Detailed balance requires that the rateRg,k,T⟶g′,k′,T′ goes from a state g,k,T to another state g′,k′,T′ and the opposite meet(22)Rg,k,T⟶g′,k′,T′Rg′,k′,T′⟶g,k,T=Pg′,k′,T′APg,k,TA=Pg′,k′,T′Pg,k,T×PAg′,k′,T′PAg,k,T,where we allow for equation (10). From (20), we have(23)Pg′,k′,T′Pg,k,T=N−2k−k′∏r=1knr′!∏r=1knr!.We considerRg,k,T⟶g′,k′,T′ as πg,k,T⟶g′,k′,T′αg,k,T⟶g′,k′,T′, where the previous part of the product represents the probability of proposing a move and the latter represents the probability chosen to satisfy the detailed balance condition for accepting the move. Then,(24)Rg,k,T⟶g′,k′,T′Rg′,k′,T′⟶g,k,T=πg,k,T⟶g′,k′,T′πg′,k′,T′⟶g,k,T×αg,k,T⟶g′,k′,T′αg′,k′,T′⟶g,k,T.BMCA is described as follows:Input: the bipartite adjacency matrix B and the node-type vector ti.Initial communities partition: using the process described in Section 2.2.2.Monte Carlo Sampling:(1)
(a) In each step of BMCA, we carry out a rearrangement of type 1 with probability1−1/N−1. If k=2 (i.e., ka=1 and kb=1), we do nothing. Otherwise, when k>2, first, we randomly select a community label r in the range 1,…,k. If the number of communities of type Tr is more than one, then we randomly select a community of type Tr labels s; otherwise, we turn to communities of another type, from which we randomly reselect a pair of communities, respectively, labels r and s. Then, we randomly select one node from community r and move it to community s. The number k=ka+kb of total communities remains constant.(b) In the process, if communityr becomes empty as a result that its last node is removed, the number of communities k decreases by one. In practice, for the type-a communities labels 1,…,ka, we can efficiently change the community ka to have label r and then change the community k to have label ka; specially, if r=ka, we only perform the latter relabeling. In such a process, the number ka of communities of type-a decreases by 1 and the number kb of communities of type-b remains constant. For type-b community labels ka+1,…,k, we can efficiently change the community k to have label r; specifically, if r=k, no relabeling is necessary. In such a process, the number kb of communities of type-b decreases by 1 and the number ka of communities of type-a remains constant. The number ka=ka+kb of total communities decreases by 1.(2)
Otherwise, we carry out a rearrangement of type 2 with probability1/N−1. We randomly select a community label r in the range 1,…,k. If there is only one node in community r, we do nothing. Otherwise, if Tr=type−a, we change community label ka+1 to k+1 and create a new empty community ka+1. Then, we randomly select a node from community r and move it to the newly created community ka+1. During this process, the number ka of communities of type-a increases by 1 and the number kb of communities of type-b remains constant. If Tr=type−b, we simply create a new empty community k+1 and no relabeling is necessary. Then, we randomly select a node from community r and move it to the newly created community k+1; during this process, the number kb of communities of type-b increases by 1 and the number ka of communities of type-a remains constant. The number k=ka+kb of total communities increases by 1.(3)
We accept the rearrangement proposed above with acceptance probability [4]:(25)αg,k,T⟶g′,k′,T′=min1,PAg′,k′,T′PAg,k,T.(4)
Repeat steps 1–3.Output: the posterior probabilitiesPkaA and PkbA.Rearrangements of any type are performed, and we always have the following equation (see Table1):(26)πg,k,T⟶g′,k′,T′πg′,k′,T′⟶g,k,T=Pg′,k′,T′Pg,k,T.Table 1
Derivation of the relation ofπg,k,T⟶g′,k′,T′/πg′,k′,T′⟶g,k,T and Pg′,k′,T′/Pg,k,T.
Type of rearrangement
States
Probability of performing the move
Probability of selecting two distinct communities of the same type
Probability of selecting a specific node from or to distinct communities of the same type
Total probabilityπ
π
g
,
k
,
T
⟶
g
′
,
k
′
,
T
′
/
π
g
′
,
k
′
,
T
′
⟶
g
,
k
,
T
P
g
′
,
k
′
,
T
′
/
P
g
,
k
,
T (using equation (23))
Type 1 withk′=k
g
,
k
,
T
⟶
g
′
,
k
′
,
T
′
1 − 1/(N − 1)
1
/
k
a
k
a
−
1
+
k
b
k
b
−
1
1
/
n
r
1
−
1
/
N
−
1
×
1
/
k
a
k
a
−
1
+
k
b
k
b
−
1
×
1
/
n
r
n
s
′
/
n
r
n
s
′
/
n
r
g
′
,
k
′
,
T
′
⟶
g
,
k
,
T
1 − 1/(N − 1)
1
/
k
a
k
a
−
1
+
k
b
k
b
−
1
1
/
n
s
′
1
−
1
/
N
−
1
×
1
/
k
a
k
a
−
1
+
k
b
k
b
−
1
×
1
/
n
s
′
Type 1 withk′=k−1
g
,
k
,
T
⟶
g
′
,
k
′
,
T
′
1 − 1/(N − 1)
1
/
k
a
k
a
−
1
+
k
b
k
b
−
1
1 (nownr = 1)
1
−
1
/
N
−
1
×
1
/
k
a
k
a
−
1
+
k
b
k
b
−
1
N
−
2
n
s
′
N
−
2
n
s
′
g
′
,
k
′
,
T
′
⟶
g
,
k
,
T
1/(N − 1)
1
/
k
a
k
a
−
1
+
k
b
k
b
−
1
1
/
n
s
′
1
/
N
−
1
×
1
/
k
a
k
a
−
1
+
k
b
k
b
−
1
×
1
/
n
s
′
Type 2 withk′=k+1
g
,
k
,
T
⟶
g
′
,
k
′
,
T
′
1/(N − 1)
1
/
k
a
+
1
k
a
+
k
b
+
1
k
b
1
/
n
r
1
/
N
−
1
×
1
/
k
a
+
1
k
a
+
k
b
+
1
k
b
×
1
/
n
r
1
/
N
−
2
n
r
1
/
N
−
2
n
r
g
′
,
k
′
,
T
′
⟶
g
,
k
,
T
1 − 1/(N − 1)
1
/
k
a
+
1
k
a
+
k
b
+
1
k
b
1 (nowns′ = 1)
1
−
1
/
N
−
1
×
1
/
k
a
+
1
k
a
+
k
b
+
1
k
bTaking into consideration equations (22) and (24), the detailed balance condition can be written as(27)αg,k,T⟶g′,k′,T′αg′,k′,T′⟶g,k,T=PAg′,k′,T′PAg,k,T.Therefore, our algorithm satisfies the detailed balance with acceptance probability of equation (25) and will sample correctly from the distribution Pg,k,TA.


## 3.2. Output of BMCA

In our implementation, a given number of stepsti per node is performed in a Monte Carlo run on a bipartite network, and then we write approximate posterior probabilities:(28)Pk,TA=the timeskshowed out  in a Monte Carlo runS.For type-a, we can get approximately(29)PkaA=the timeskashowed out in a Monte Carlo runS,and for type-b,(30)PkbA=the timeskbshowed out in a Monte Carlo runS.The most likely number of type-a communities in a bipartite network is ka with the biggest value PkaA, and the number of type-b communities kb can be given in the same way.In order to avoid any bias in the results and improve the correctness of BMCA, instead of just usingPkaA, we performed a given number M of Monte Carlo runs for a bipartite network. Thus, we obtained the average value of each PikaA as the final posterior probabilities:(31)P′kaA=∑i=1MPikaAM,for type-a and(32)P′kbA=∑i=1MPikbAM,for type-b.We takeON+Na∗Nb to calculate the acceptance probability of each node move, so the time complexity of our algorithm is ONa2Nb+NaNb2.


## 4. Example Application

Here, we demonstrate our algorithm on synthetic bipartite networks generated by the bipartite stochastic block model [6] and find that it works well in most cases.There is often some noise in empirically observed networks because of both errors in the measurements and missing data [26]. Therefore, we employ a mixing model to generate noisy synthetic networks for testing the robustness of our algorithm. In this model, we specify g, and the expected value of an edge is given by(33)ω=λωplanted+1−λωrandom,where the parameters ωplanted and ωrandom are used to generate a pure planted community structure and no community structure, respectively; and the mixing parameter λ∈0,1 is used to control various levels of uniformly random noise. Following Larremore et al. [6], we let ωrsplanted=mrs and ωrsrandom=κrκs/m, where θi=di/κgi. We employ this model to create synthetic networks of an easy case with a homogeneous degree distribution and a difficult case with a heterogeneous degree distribution.We performedM=10 Monte Carlo runs for each network of S=10000 steps per node and found that BMCA can determine the correct number of communities for synthetic networks and outperform the projection method in every case.


### 4.1. An Easy Case

In the easy case, we use the model above to create the synthetic networks including four type-a communities and four type-b communities, i.e., ka=kb=4, and each node has the same degree (i.e., the network with planted community structure has a homogeneous degree distribution). We let the number of each type of node equal to 1000, and all communities are equally sized as 250. Then, we let m1,5=m2,6=m3,7=m4,8=2500 (i.e., the total number of edges is 10000), and let the block structure matrix ωplanted be defined with ω1,5planted=ω2,6planted=ω3,7planted=ω4,8planted=2500 (the symmetric entry has the same value). Moreover, the random structure matrix ωrandom can be defined with ωi,jrandom=625i=1,2,3,4;j=5,6,7,8. Finally, with these specifications, we create networks of the easy case to test our algorithm (the code to create these networks can be downloaded from [27]).As the mixing parameterλ increases, i.e., the level of noise is decreased, BMCA begins to estimate the correct number of communities for the network in an easy case when λ=0.5; in addition, the fraction of correct communities number of the network calculated by BMCA increases as a whole (blue line in Figure 1). Then, we use our method to derive the synthetic mixing networks generated with λ=0.6 and λ=0.65, indicated by red circles in Figure 1(a) and Figure 1(b) and showing the posterior probabilities of the number of communities in the network in Figure 2(a) and Figure 2(b). As shown in Figures 1 and 2, we are able to estimate the correct number of different types of communities when λ≥0.65 with PBMCA′ka=4A=0.6 (PBMCA′kb=4A=0.6); then, when λ≥0.8, the proposal probability of the correct number of communities PBMCA′ka=4A and PBMCA′kb=4A is equal to 1. However, the projection method gives poorer results (see the red lines in Figure 1), and it begins to estimate the correct number of different types of communities at λ=0.6 with PPM′ka=4A=0.1 (PPM′kb=4A=0.1) and when λ=1 PPM′ka=4A=0.9 (PPM′kb=4A=0.8). When λ≥0.5, PBMCA′ka=4A (PBMCA′kb=4A) is always bigger than PPM′ka=4A(PPM′kb=4A), as shown in Figure 3.Figure 1
Test of BMCA against the projection method (PM) on synthetic networks generated using the bipartite stochastic block model in an easy case. Posterior probabilitiesP′ka=4A (a) and P′kb=4A (b) are functions of the mixing parameter λ.
(a)
(b)Figure 2
Posterior probabilitiesP′kaA of the number of communities calculated for the synthetic networks generated using the bipartite stochastic block model in an easy case. The synthetic networks were generated with λ=0.6 (a) and λ=0.65 (b) in the easy case.
(a)
(b)Figure 3
Test of BMCA against the projection method (PM) on synthetic networks generated using the bipartite stochastic block model in a difficult case. Posterior probabilitiesP′ka=2A (a) and P′kb=3A (b) are functions of the mixing parameter λ.
(a)
(b)


### 4.2. A Difficult Case

In the difficult case, the synthetic networks are created with two type-a communities ka=2 and three type-b communities ka=3, and the degree of each node is different (i.e., the network with planted community structure has a heterogeneous degree distribution). The communities are set with different sizes, and we divide 700 type-a nodes evenly into 2 communities {350, 350} and 300 type-b nodes into 3 communities {100, 150, 150}. Let m1,3=m2,4=2500 and m1,5=m2,5=1500; i.e., the total number of edges is 8000. Then, the block structure matrix ωplanted can be defined with ω1,3Planted=ω2,4Planted=2500 and ω1,5Planted=ω2,5Planted=1500 (the symmetric entry has the same value), and the random network matrix ωrandom can be defined with ω1,3random=ω1,4random=ω2,3random=ω2,4random=1250 and ω1,5random=ω2,5random=1500. The symmetric entry has the same value. Finally, with the specification above, we create networks of the difficult case to test our algorithm (the code to create these networks can be downloaded from [27]).As shown by the blue line in Figure3, when the level of noise is decreased, BMCA begins to estimate the correct number of communities for the network in a difficult case when λ=0.5, and the posterior probabilities P′ka=2A and P′kb=3A increase sharply after λ≥0.55. Especially, when λ≥0.75, the proposal probability PBMCA′ka=2A=1 PBMCA′kb=3A=1; i.e., we are always able to determine the correct number of communities. We used our method to calculate the probability of the number of communities in the synthetic mixing networks generated with λ=0.55 and λ=0.6, indicated as red circles in Figure 3, and show the result in Figure 4. As seen from Figures 3 and 4, our method can estimate the correct number of communities ka=2 and kb=3 in the synthetic network when λ≥0.6 for the difficult case. However, the projection method fails to estimate the correct number of communities even when λ=1 and there is no noise, as shown by the red line in Figure 3.Figure 4
Posterior probabilities of the number of communities calculated for the synthetic networks generated using the bipartite stochastic block model in a difficult case. Posterior probabilitiesP′kaA of the number of type-a communities in the synthetic networks generated with λ=0.55 (a) and λ=0.55 (b) and P′kbA of the number of type-b communities in the synthetic networks generated with λ=0.55 (c) and λ=0.6 (d).
(a)
(b)
(c)
(d)


### 4.3. Further Testing

We tested our method on synthetic networks of two different sizes of community as the number of network communities increases. The networks were generated using the bipartite stochastic block model withλ=1, and the other parameters are set as listed in Table 2, where θi=1,…,N is heterogeneous as real-world networks and the mean node degree of the network for the figures is 10. The number of communities ka or kb estimated using BMCA was correct until the actual number of communities increased to about 7, which is about 14 for k. The results are shown in Figure 5.Table 2
The parameters set for synthetic networks.
Figure no.
n
r
r
=
1
,
…
,
k
θ
i
i
=
1
,
…
,
N
ω
r
s
,
r
=
1
,
…
,
k
a
s
=
r
+
k
a and ka=1,…,10
Figure5(a)
250
Homogeneous
2500
Figure5(b)
250
Heterogeneous
2500
Figure5(c)
500
Homogeneous
5000
Figure5(d)
500
Heterogeneous
5000Figure 5
Tests of our algorithm on synthetic bipartite networks. In each subgraph, triangles represent results derived from Monte Carlo runs with the known correct partitions (ground truth initial conditions), while circles represent runs started with a random initial partition of nodes to community (random initial conditions).
(a)
(b)
(c)
(d)The results show a tendency to underestimate the number of communities for higher actual numbers of communities, especially when the size of the community changes from 250 (Figures5(a) and 5(b)) to 500 (Figures 5(c) and 5(d)) and that of the networks increase correspondingly. However, these calculations appear when BMCA is run with a random initialization partition of nodes to communities. When BMCA is started on the same network with community partitions corresponding exactly to the planted community structure (yellow triangles), we always find the accurate number of communities.Even so, the underestimation ofka (or kb) occurs not because the correct community partition fails to maximize the posterior probability Pka∣A but rather because BMCA has not run for long enough to find the maximum. The method is theoretically sound, but when the number of possible community partitions kN increases very rapidly with k, the numerical calculation becomes too demanding [4]. We can possibly design a more efficient Monte Carlo algorithm to solve this problem although it offers some useful information of a lower bound on the number of communities in the network given by BMCA.


## 4.1. An Easy Case

In the easy case, we use the model above to create the synthetic networks including four type-a communities and four type-b communities, i.e., ka=kb=4, and each node has the same degree (i.e., the network with planted community structure has a homogeneous degree distribution). We let the number of each type of node equal to 1000, and all communities are equally sized as 250. Then, we let m1,5=m2,6=m3,7=m4,8=2500 (i.e., the total number of edges is 10000), and let the block structure matrix ωplanted be defined with ω1,5planted=ω2,6planted=ω3,7planted=ω4,8planted=2500 (the symmetric entry has the same value). Moreover, the random structure matrix ωrandom can be defined with ωi,jrandom=625i=1,2,3,4;j=5,6,7,8. Finally, with these specifications, we create networks of the easy case to test our algorithm (the code to create these networks can be downloaded from [27]).As the mixing parameterλ increases, i.e., the level of noise is decreased, BMCA begins to estimate the correct number of communities for the network in an easy case when λ=0.5; in addition, the fraction of correct communities number of the network calculated by BMCA increases as a whole (blue line in Figure 1). Then, we use our method to derive the synthetic mixing networks generated with λ=0.6 and λ=0.65, indicated by red circles in Figure 1(a) and Figure 1(b) and showing the posterior probabilities of the number of communities in the network in Figure 2(a) and Figure 2(b). As shown in Figures 1 and 2, we are able to estimate the correct number of different types of communities when λ≥0.65 with PBMCA′ka=4A=0.6 (PBMCA′kb=4A=0.6); then, when λ≥0.8, the proposal probability of the correct number of communities PBMCA′ka=4A and PBMCA′kb=4A is equal to 1. However, the projection method gives poorer results (see the red lines in Figure 1), and it begins to estimate the correct number of different types of communities at λ=0.6 with PPM′ka=4A=0.1 (PPM′kb=4A=0.1) and when λ=1 PPM′ka=4A=0.9 (PPM′kb=4A=0.8). When λ≥0.5, PBMCA′ka=4A (PBMCA′kb=4A) is always bigger than PPM′ka=4A(PPM′kb=4A), as shown in Figure 3.Figure 1
Test of BMCA against the projection method (PM) on synthetic networks generated using the bipartite stochastic block model in an easy case. Posterior probabilitiesP′ka=4A (a) and P′kb=4A (b) are functions of the mixing parameter λ.
(a)
(b)Figure 2
Posterior probabilitiesP′kaA of the number of communities calculated for the synthetic networks generated using the bipartite stochastic block model in an easy case. The synthetic networks were generated with λ=0.6 (a) and λ=0.65 (b) in the easy case.
(a)
(b)Figure 3
Test of BMCA against the projection method (PM) on synthetic networks generated using the bipartite stochastic block model in a difficult case. Posterior probabilitiesP′ka=2A (a) and P′kb=3A (b) are functions of the mixing parameter λ.
(a)
(b)


## 4.2. A Difficult Case

In the difficult case, the synthetic networks are created with two type-a communities ka=2 and three type-b communities ka=3, and the degree of each node is different (i.e., the network with planted community structure has a heterogeneous degree distribution). The communities are set with different sizes, and we divide 700 type-a nodes evenly into 2 communities {350, 350} and 300 type-b nodes into 3 communities {100, 150, 150}. Let m1,3=m2,4=2500 and m1,5=m2,5=1500; i.e., the total number of edges is 8000. Then, the block structure matrix ωplanted can be defined with ω1,3Planted=ω2,4Planted=2500 and ω1,5Planted=ω2,5Planted=1500 (the symmetric entry has the same value), and the random network matrix ωrandom can be defined with ω1,3random=ω1,4random=ω2,3random=ω2,4random=1250 and ω1,5random=ω2,5random=1500. The symmetric entry has the same value. Finally, with the specification above, we create networks of the difficult case to test our algorithm (the code to create these networks can be downloaded from [27]).As shown by the blue line in Figure3, when the level of noise is decreased, BMCA begins to estimate the correct number of communities for the network in a difficult case when λ=0.5, and the posterior probabilities P′ka=2A and P′kb=3A increase sharply after λ≥0.55. Especially, when λ≥0.75, the proposal probability PBMCA′ka=2A=1 PBMCA′kb=3A=1; i.e., we are always able to determine the correct number of communities. We used our method to calculate the probability of the number of communities in the synthetic mixing networks generated with λ=0.55 and λ=0.6, indicated as red circles in Figure 3, and show the result in Figure 4. As seen from Figures 3 and 4, our method can estimate the correct number of communities ka=2 and kb=3 in the synthetic network when λ≥0.6 for the difficult case. However, the projection method fails to estimate the correct number of communities even when λ=1 and there is no noise, as shown by the red line in Figure 3.Figure 4
Posterior probabilities of the number of communities calculated for the synthetic networks generated using the bipartite stochastic block model in a difficult case. Posterior probabilitiesP′kaA of the number of type-a communities in the synthetic networks generated with λ=0.55 (a) and λ=0.55 (b) and P′kbA of the number of type-b communities in the synthetic networks generated with λ=0.55 (c) and λ=0.6 (d).
(a)
(b)
(c)
(d)


## 4.3. Further Testing

We tested our method on synthetic networks of two different sizes of community as the number of network communities increases. The networks were generated using the bipartite stochastic block model withλ=1, and the other parameters are set as listed in Table 2, where θi=1,…,N is heterogeneous as real-world networks and the mean node degree of the network for the figures is 10. The number of communities ka or kb estimated using BMCA was correct until the actual number of communities increased to about 7, which is about 14 for k. The results are shown in Figure 5.Table 2
The parameters set for synthetic networks.
Figure no.
n
r
r
=
1
,
…
,
k
θ
i
i
=
1
,
…
,
N
ω
r
s
,
r
=
1
,
…
,
k
a
s
=
r
+
k
a and ka=1,…,10
Figure5(a)
250
Homogeneous
2500
Figure5(b)
250
Heterogeneous
2500
Figure5(c)
500
Homogeneous
5000
Figure5(d)
500
Heterogeneous
5000Figure 5
Tests of our algorithm on synthetic bipartite networks. In each subgraph, triangles represent results derived from Monte Carlo runs with the known correct partitions (ground truth initial conditions), while circles represent runs started with a random initial partition of nodes to community (random initial conditions).
(a)
(b)
(c)
(d)The results show a tendency to underestimate the number of communities for higher actual numbers of communities, especially when the size of the community changes from 250 (Figures5(a) and 5(b)) to 500 (Figures 5(c) and 5(d)) and that of the networks increase correspondingly. However, these calculations appear when BMCA is run with a random initialization partition of nodes to communities. When BMCA is started on the same network with community partitions corresponding exactly to the planted community structure (yellow triangles), we always find the accurate number of communities.Even so, the underestimation ofka (or kb) occurs not because the correct community partition fails to maximize the posterior probability Pka∣A but rather because BMCA has not run for long enough to find the maximum. The method is theoretically sound, but when the number of possible community partitions kN increases very rapidly with k, the numerical calculation becomes too demanding [4]. We can possibly design a more efficient Monte Carlo algorithm to solve this problem although it offers some useful information of a lower bound on the number of communities in the network given by BMCA.


## 5. Conclusions

In our paper, a new projection-free Bayesian inference method for determining the number of pure-type communities in a bipartite network has been introduced. First, we present the first principle derivation of a practical method, using the degree-corrected bipartite stochastic block model that is able to deal with networks with broad degree distributions, for estimating the number of pure-type communities of bipartite networks. Second, we propose a prior probability distribution over the partition of a bipartite network, with type parameterT ensuring that each community is pure type. Third, we design a Monte Carlo algorithm incorporated with our proposed method and prior probability distribution. We have illustrated the performance of the method with applications to a wide range of synthetic bipartite networks, including an easy case with homogeneous degree distributions and a difficult case with heterogeneous degree distributions. The results show that the proposed algorithm can determine the correct number of communities and perform better than the projection method, especially in networks with heterogeneous degree distributions.However, our method underestimates the number of communities when the number of communities becomes large. The reason for this is due to the number of possible community partitions increasing very rapidly with the increase in the number of communities and because our algorithm has not run for long enough to find the posterior probability. Thus, our future work will focus on (i) finding a method that can more efficiently sample the posterior distribution over community partitions to correctly estimate a large number of communities in the network and (ii) extending applications on real-world data sets.



---

*Source: 9471201-2019-12-09.xml*