Journal of Machine Learning Research 1 (2000) 1-48

Submitted 4/00; Published 10/00

Collapsed Variational Bayes Inference of Infinite Relational Model
Katsuhiko Ishiguro

ishiguro.katsuhiko@lab.ntt.co.jp

NTT Communication Science Laboratories
NTT Corporation
Kyoto 619-0237, Japan

arXiv:1409.4757v1 [cs.LG] 16 Sep 2014

Issei Satosato@r.dl.itc.u-tokyo.ac.jp

Information Technology Center
The University of Tokyo
Tokyo 113-0033, Japan

Naonori Ueda

ueda.naonori@lab.ntt.co.jp

NTT Communication Science Laboratories
NTT Corporation
Kyoto 619-0237, Japan

Editor: Unknown Editor

Abstract
The Infinite Relational Model (IRM) is a probabilistic model for relational data clustering that partitions objects into clusters based on observed relationships. This paper presents Averaged CVB(ACVB) solutions for IRM, convergence-guaranteed and practically useful fast Collapsed Variational Bayes (CVB) inferences. We first derive ordinary CVB and CVB0 for IRM based on thelower bound maximization. CVB solutions yield deterministic iterative procedures for inferring
IRM given the truncated number of clusters. Our proposal includes CVB0 updates of hyperparameters including the concentration parameter of the Dirichlet Process, which has not been studied in
the literature. To make the CVB more practically useful, we further study the CVB inference in two
aspects. First, we study the convergence issues and develop a convergence-guaranteed algorithm
for any CVB-based inferences called ACVB, which enables automatic convergence detection andfrees non-expert practitioners from difficult and costly manual monitoring of inference processes.
Second, we present a few techniques for speeding up IRM inferences. In particular, we describe
the linear time inference of CVB0, allowing the IRM for larger relational data uses. The ACVB solutions of IRM showed comparable or better performance compared to existing inference methods
in experiments, and provide deterministic, faster, and easier convergence detection.
Keywords: nonparametric Bayes, infinite relational models, collapsed variational Bayes inference, averaged CVB, relational data analysis

1. Introduction
Analysis of pairwise relational data, such as friend-links on social network services (SNS), customer records of purchases in online shops, and bibliographic citations between scientific articles,
is useful in many ways. Many statistical models for relational data have been presented in the literature (Clauset et al., 2008; Erosheva et al., 2004; Liben-Nowell and Kleinberg, 2003; Zhu et al.,2009). Among them, the infinite relational model (IRM) proposed by Kemp et al. (2006) achieves
simultaneous bi-clustering on the row and column dimensions of a given pairwise relational data
c 2000 Katsuhiko Ishiguro, Issei Sato, and Nanori Ueda.

Ishiguro, Sato, and Ueda

matrix. For example, in the case of customer records, rows and columns correspond to users anditems. In such a case, the row and column clusters are interpreted as latent user groups and item
topics, respectively. IRM adopts nonparametric Bayes modeling and so can automatically estimate
the number of clusters. This makes IRM a convenient tool for relational data analysis without theneed for careful model selection.
Two Bayesian inference algorithms are frequently used for probabilistic generative models including IRM: the Gibbs sampler and variational Bayes. The former guarantees asymptotic convergence to the true posteriors of random variables given infinitely many stochastic samples. Variational Bayes (VB) solutions often enjoy faster convergence with deterministic iterative computations
and massively parallel computation thanks to the factorization. The VB approaches also allow easy
and automatic detection of convergence. Instead of these favorable properties, VB only yields local
optimal solutions due to factorized approximated posteriors.
We can improve these inference methods by developing collapsed estimators, which integrateout some parameters from inferences. Collapsed Gibbs samplers are one of the best inference
solutions as they achieve faster convergence and better estimation than the original Gibbs samplers.
Recently, collapsed variational Bayes (CVB) solutions have been intensively studied, especially for
topic models such as latent Dirichlet allocation (LDA) (Teh et al., 2007; Asuncion et al., 2009;
Sato and Nakagawa, 2012) and HDP-LDA (Sato et al., 2012). The original paper (Teh et al., 2007)
examined a 2nd-order Taylor approximation of the variational expectation. A simpler 0th orderapproximated CVB (CVB0) solution also has been developed; it is an optimal solution in the senseof minimized α-divergence (Sato and Nakagawa, 2012). These papers report that CVB and CVB0
yield better inference results than VB solutions, even slightly better than exact collapsed Gibbs, indata modeling (Kurihara et al., 2007; Teh et al., 2007; Asuncion et al., 2009), link predictions, andneighborhood search (Sato et al., 2012).
Most IRM papers to date (Kemp et al., 2006; Ishiguro et al., 2012; Mø rup et al., 2010; Albers
et al., 2013) rely on (collapsed) Gibbs samplers. However, the automatic convergence detectionof stochastic sampling-based Gibbs is difficult to achieve (Cowles and Carlin, 1996). This is notpreferable for non-expert users to employ IRM in practical uses. Further, (Albers et al., 2013)
reported that the naive implementation of (collapsed) Gibbs is very slow in mixing for IRM applications. However, interestingly, there has been no attempt to use VB for IRM to the best of ourknowledge, even though VB allows easy and automatic detection of convergence, plus fast deterministic computations. One reason is that VB may perform poorly for IRM because IRM solvesdifficult partitioning problems with many local optima. CVB and CVB0 are promising alternatives
to VB, but most CVB studies have focused on topic models, which well suits simple Bag-of-Wordstyle data sets. The only exceptions are CVB for Probabilistic Context-Free Grammars (Wang andBlunsom, 2013b) and Hidden Markov Models (Wang and Blunsom, 2013a).
In this paper, we first formulate and derive the CVB inference of IRM for relational data analysis
as fast, deterministic and precise inference algorithms, which replace naive VB. Furthermore, we
derive update rules of hyperparameters based on CVB0; thus we can automatically optimize all
hyperparameters. In particular, the update of the concentration parameter of DP has not been studied
in the literature, which plays an important role in nonparametric Bayes. In Table 1, we summarize
the existing inference algorithms used in the previous IRM studies and this paper.
Next, we study CVB inference of IRM in two aspects to make the CVB inference easier to use
for practitioners. The first aspect is the convergence. The convergence behavior of CVB inferenceis still difficult to analyze theoretically, and there is no convergence guarantee for the general CVB
2

CVB Inference of IRM

Table 1: Inference algorithms presented in existing IRM works and this paper.
Paper(collapsed) Gibbs VB CVB CVB0 Comments
Kemp et al. (2006)
X
The seminal paper
Mø rup et al. (2010)
X
Application to fMRI
Hansen et al. (2011)X
GPU impl.
X
OpenMP impl.
Albers et al. (2013)
Ishiguro et al. (2012)
X
Noise filtering extension
This paper
X
X
X
X
Fully covers

inferences. This problem, interestingly, has not been much discussed in the literature. The sole
exception is (Foulds et al., 2013), which uses online stochastic learning valid for LDA. However,this is a tricky and problematic issue for practitioners who are not familiar with but want to try stateof-the-art machine learning techniques. Users are required to determine the convergence of CVB
inference manually: this is not an easy task for non-expert users. In that sense, CVB is not as favorable as naive VB and EM algorithms. In this paper, we empirically study the convergence behaviorsof CVB inference for IRM. We first monitor the naive variational lower bound and the pseudo
leave-one-out training log likelihood, and empirically show that both may serve as convergence detectors. Then, we develop a simple and effective technique that assures convergence of CVB for
any probabilistic model including IRM. The proposed annealing technique called Averaged CVB(ACVB) guarantees the convergence of CVB and allows automatic convergence detection. ACVB
has two advantages. First, ACVB posterior update offers assured convergence thanks to its simple
annealing mechanism. Second, the stationary point of the CVB lower bound is equivalent to theconverged solution of ACVB, if the lower bound has a stationary point (an issue unresolved in theliterature). Our formulation is applicable to any model, and is equally valid for CVB and CVB0.
Convergence-guaranteed ACVB is the preferred choice for practitioners who want to apply stateof-the-art inference to their problems. In Table 2, we summarize the existing CVB works and thispaper, based on applied models and the convergence issue.
The second aspect is computational cost. In naive implementation, CVB inference of IRM
requires square time over the number of objects. For example, if the buy-product relations X areobserved between N1 users and N2 items, and we assume K1 and K2 latent clusters among users
and items respectively, then the inference costs O (K1 K2 N1 N2 ) per iteration. This makes the IRM
an impractical solution for large relational data. In this paper, we describe how to mitigate thiscomputational cost, especially for the CVB0 solution. As a result, we show that we can solve the
CVB0 solution by O (K1 K2 L(N1 + N2 )), linear to the number of users and items, where L denotes
the average degree. Combining these techniques, we propose a practically useful ACVB0 solution
for IRM, with easy detection of guaranteed convergence and linear-time computation.
We experimentally show that ACVB solutions for IRM offer comparable, or even better performance in terms of data modeling (test likelihood) than naive variational Bayes on multiple synthetic
and real-world relational data sets. ACVB0 convergences are magnitude faster than the VB and
ACVB in terms of CPU times, presented stable and nice convergence behaviors throughout thedatasets. In addition, we demonstrate the scalability of the proposed ACVB0 solution to larger relational data by employing the linear time inference algorithm. Based on these findings, we conclude
that the ACVB0 inference solution on IRM is convenient and appealing for practitioners who work3

Ishiguro, Sato, and Ueda

Table 2: CVB-related studies summary: in terms of applied models and convergence.
Paper
Applied model
Convergence Comments
Teh et al. (2007)
LDA-CVB
The seminal paper
Asuncion et al. (2009)
LDA-CVB, CVB0
Introduces CVB0
Sato and Nakagawa (2012)
LDA-CVB, CVB0
Optimality analysis by
α-divergence
Foulds et al. (2013)
LDA-CVB0X
Stochastic approx. valid
for LDA
Kurihara et al. (2007)
DPM-CVBFirst attempt at DPM
Teh et al. (2008)
HDP-CVBFirst attempt at HDPSato et al. (2012)
HDP-CVB, CVB0
Approx. solution
Wang and Blunsom (2013a)
HMM-CVB0First attempt at HMM
Wang and Blunsom (2013b)
PCFG-CVB0First attempt at PCFG
This paper
IRM-CVB, CVB0
XFirst attempt at IRM,
convergence assurance
for any modelwith relational data; it shows good modeling performance, assures automatic convergence, and is
fast by linear time inference.
The contributions of this paper are three-fold.
1. We first present Collapsed Bayes solutions (CVB and CVB0) for inference of IRM, which
is used for relational data analysis. The CVB solutions are fast, precise and deterministic
inference algorithms. We also present update rules of hyperparameters, including the concentration parameter of the Dirichlet Process.
2. We empirically study the convergence behaviors of CVB solutions. Along with that, we
propose a simple but effective annealing technique called Averaged CVB (ACVB) that assuresthe convergence of the CVB inference for any model.
3. We show techniques to speed up the (A)CVB inference. In particular, one of them allows us
to solve (A)CVB0 of IRM in linear time to the number of objects. This linear time algorithm
is more effective when the relational data is sparse, which is typical in real datasets.
The rest of this paper is organized as follows. In the next section, we introduce the IRM model.
In the third section, we briefly review collapsed Gibbs sampler solutions of IRM and related works.
We introduce a naive VB solution of IRM in the fourth section. As stated, the VB solution is
not good for IRM, but it serves as a stepping stone for the CVB. Section 5 presents our CVB
solutions. Section 6 presents convergence issues, including convergence-assured annealing technique (ACVB). Section 7 discusses the speeding-up technique, including linear time inference of(A)CVB0. The eighth section is devoted to experimental evaluations, and the final section concludes
the paper.
4CVB Inference of IRM

2. Infinite Relational Models and Related WorksFirst, we introduce the infinite relational model (IRM) (Kemp et al., 2006), which estimates theunknown number of hidden clusters within a set of relational data. Then, we review some related
works.
2.1 Dirichlet Process Mixture
In the infinite relational model (IRM) (Kemp et al., 2006), the Dirichlet process (DP) is used asa prior for clusters of an unknown number; intuitively, it is equivalent to an infinite-dimensional
Dirichlet distribution. The Dirichlet Process Mixture (DPM) is a probabilistic generative modelthat uses DP for the prior of mixture proportions. We can implement DPM by using either a StickBreaking Process (SBP) (Sethuraman, 1994) or a Chinese restaurant process (CRP) (Blackwell and
MacQueen, 1973), which is a marginalized form of SBP. CRP is employed for the (collapsed) Gibbs
sampler, and SBP is employed for (collapsed) variational Bayes solutions typically.
First, let us start by explaining CRP. CRP is introduced as a probability distribution over apartitioning of N objects. Let zi = k, i ∈ {1, . . . , N}, k ∈ {1, . . . , K} denote the assignment of ith
object to the kth partition (cluster) among the total of K partitions. Then, the CRP is defined by thefollowing equations:QK(mk − 1)!
K
CRP(z1:N |α) = α QNk=1,(1)
i=1 (α + i − 1)
 \i


 mk

m\i
N−1+α
k > 0,p(zi = k|z\i , α) = 
(2)

α

 N−1+α=0.
m\ik
α > 0 is a hyperparameter called a concentration parameter. Equation 1 shows the joint probabilityof K partitions. The equation is rewritten as Eq. (2), which gives the probability of object i being
allocated to partition k given the assignments of other objects. mk denotes the number of objects
assigned to partition k, and m\ik denotes the same number excluding object i. The first part of Eq. (2)indicates that object i will be assigned to existing cluster k with a probability proportional to its
membership. The second part of Eq. (2) indicates that object i will be assigned to a new cluster witha probability proportional to α > 0. Repeatedly applying Eq. (2) for each object, we can randomly
generate cluster partitions of N objects. For each run of CRP, we will have different clustering ofobjects, and will also have a different number of resulting clusters.
Next, we explain SBP. SBP is another construction of DPM. In SBP, we explicitly sample an
infinite-dimensional vector of mixture proportions, while CRP directly samples cluster assignment
z without such a vector. SBP construction of DPM is described as follows:vk ∼ Beta (1, α) ,πk = vkv−1
Y

k = 1, . . . , ∞(1 − vl ) ,(3)

k = 1, . . . , ∞(4)

i = 1, . . . , N

(5)

l=1

zi ∼ Multinomial (π) .

Computing π = {π1 , π2 , . . . }, a mixing proportion vector of the mixture, is explained by an analogyof ”stick breaking”. We first assume a stick of length 1, which is the full portion of clusters. Then,
5

Ishiguro, Sato, and Ueda

we break the stick into a stick of length v1 (Eq. (3)) and the remainder. The proportion of the firstcluster π1 is the length of this stick, that is, v1 (Eq. (4)). Then, we again break the remaining stick,which has a length of 1 − v1 , based on the ratio of v2 (Eq. (3)). The proportion of the second cluster
π2 is equivalent to the broken stick, that is, v2 (1 − v1 ) (Eq. (4)). Repeating this process, we obtain
the infinite-dimensional vector π, whose sum equals 1, implying an infinite number of mixture
components. Using π as a mixing proportion, we can sample cluster assignments as in Eq. (5).
2.2 Generative Models of IRMs
Next, we describe a probabilistic generative model of IRM. We assume two-place relations throughout the paper, but extension to cover higher-order relations is straightforward.
IRM is an application of DP for relational data. Let us first assume a binary two-place relationon the two sets (domains) of objects, namely D1 × D2 → {0, 1}, where D1 = {1, . . . , i, . . . , N1 } andD2 = {1, . . . , j, . . . , N2 }. IRM divides the set of objects into multiple clusters based on the observed
relational data matrix of X = {xi, j ∈ {0, 1}}. Data entry xi, j ∈ {0, 1} denotes the existence of a relation
between a row (the first domain) object i ∈ {1, 2, . . . , N1 } and a column (the second domain) objectj ∈ {1, 2, . . . , N2 }. In an online purchase record case, the first domain corresponds to a user list, andan object i denotes a specific user i. The second domain corresponds to a list of product items, andan object j denotes a specific item j. The data entry xi, j represents a relation between the user i andthe item j: namely, the purchase record.
For such data, we define an IRM as follows:θk,l |ak,l , bk,l ∼ Beta ak,l , bk,l ,(6)

z1,i |α1 ∼ CRP (α1 ) ,(7)

z2, j |α2 ∼ CRP (α2 ) ,
xi, j |Z 1 , Z 2 , {θ} ∼ Bernoulli θz1,i ,z2, j .

(8)
(9)

In Eq. (6), θk,l is the strength of the relation between cluster k of the first domain and cluster l in
the second domain. z1,i in Eq. (7) and z2, j in Eq. (8) denote the cluster assignments in the firstdomain and in the second domain, respectively. Throughout the paper, we interchangeably choose
the 1-of-K representation of Z: z1,i = k is equivalently represented by z1,i,k = 1, z1,i,l,k = 0. For
each domain, we have a CRP prior; this indicates that each domain may have a different numberof clusters. Generating the observed relational data xi, j follows Eq. (9) conditioned by the cluster
N1
2
assignments Z 1 = {z1,i }i=1
, Z 2 = {z2, j }Nj=1
and the strengths, θ. A typical example of IRM is shownin Fig. 1. IRM infers appropriate cluster assignments of objects Z 1 = {z1,i } and Z 2 = {z2, j } given theobservation relation matrix X = {xi, j }. We can interpret the clustering as the permutation of object
indices so as to discover the “block” structure as in Fig. 1 (b).
2.2.1 Single-domain Model
As a special case, we can build an IRM for a binary two-place relation between the same domain
objects D × D → {0, 1}.
6

CVB Inference of IRM

2nd domain object j (sorted)
2nd domain object j

l=1

l=2

l=3

k=1

1st domain object i (sorted)

1st domain object i

k=2
k=3

(a)

(b)

Figure 1: Example of Infinite Relational Models (IRM). (a) Input observation X. (b) A visualizationof inferred clusters Z.

The probabilistic generative model of the single-domain IRM is described as follows:θk,l |ak,l , bk,l ∼ Beta ak,l , bk,l ,(10)

zi |α ∼ CRP (α) ,

(11)




xi, j |Z, {θ} ∼ Bernoulli θzi ,z j .

(12)

The generative model clearly shows the difference of the multi-domain IRM (Eqs. (6-9)) and thesingle-domain IRM (Eqs. (10-12)). In the single-domain IRM, there are only N objects in thedomain D, and they serve as either from-nodes or to-nodes in the network. Object indices i and jpoint to the same domain. On the other hand, the multi-domain IRM distinguishes the first domain
object i and the second domain object j.
Let us explain the difference by a simple SNS example. Imagine we are given the SNS relation
data where N users are mutually interconnected to others. An observed relation is a binary value
xi, j , which indicates there is a link from user i to user j. In the case of (multi-domain) IRM, thefirst domain is a collection of N users who act as from-nodes. The second domain is a collectionof N users who act as to-nodes. We assume that each user has a different “role” in reaching a linkto others (the first domain) and accepting a link from others (the second domain). Contrarily, in thecase of the single-domain IRM, each user i is assigned with a single cluster assignment. The user
has her own “role” in the network, and this role is used in either reaching a link or accepting a link.
Obviously, the single-domain IRM is not applicable when the number of from-nodes and tonodes are different. Thus, the model has a specific and limited applicability compared to the multidomain IRM. Afterward, we focus on the multi-domain IRM in this paper, but whole discussions
are also valid for single-domain IRM.
2.3 Related Works
IRM (Kemp et al., 2006) is a rather old model for formulating general relationship observations. One
drawback of IRM is that it allows a node to have only one cluster assignment. Mixed Membership7

Ishiguro, Sato, and Ueda

Stochastic Blockmodel (MMSB) (Airoldi et al., 2008) is a finite-cluster model that allows the nodes
to have multiple cluster assignment, and change the clusters edge by edge. (Miller et al., 2009)
employs the Indian Buffet Process (IBP, see (Griffiths and Ghahramani, 2011) for a review) tohandle countably infinite binary factors for each node. The Infinite Latent Attribute Model (Palla
et al., 2012) further allows flexible modeling of networks where an infinitely many “views” havetheir own clustering of nodes. Compared to these models, IRM is inferior in the potential modelingcapability but is superior in easier interpretation of the clustering results: after all, “multiple clusters
for a single node” is somewhat counter-intuitive for non-expert users. Of course, IRM is a simpler
model; thus the inference scheme is also much simpler than these advanced models.
Recently, probabilistic models for pure networks have attracted much attention in machine learning. By “pure” networks, we mean that observed relations are limited to the single-domain case: that
is, T × T → {1, 0}. (Ho et al., 2011) introduces a nested Chinese Restaurant Process (nCRP) (Bleiet al., 2010) to incorporate multiscale membership for the MMSB. (Ho et al., 2012) proposed abag of triangular representations of a network. The representation is based on the triplet of nodes.
Possible connections among three nodes are (i) all three nodes are connected in a circuit, (ii) all
three nodes are connected in a line (no link between a specific pair of nodes), (iii) two nodes areconnected and one is isolated, and (iv) all nodes are separated. (Yin et al., 2013) combines thetriangular representations with a simpler probabilistic generative model to achieve a scalable algorithm for large networks, which limits the cardinality of the triangle “cluster assignments” variety
in the likelihood function. (Yang et al., 2013) employ the edge structure and node attributes to findcommunities within large networks. The model is called CESNA, consisting of a soft-max-based
binary node attribute model and an affiliated network model (Yang and Leskovec, 2013). These
recent works make the model scalable against very large networks consisting of millions of nodes.
However, none of these works consider the cross-domain T 1 × T 2 relational observations that are
mainly discussed in this IRM paper.
As briefly described in the introduction, collapsed variational Bayes (CVB) solutions have been
intensively studied, especially for topic models such as latent Dirichlet allocation (LDA) (Teh et al.,2007; Asuncion et al., 2009; Sato and Nakagawa, 2012) and HDP-LDA (Sato et al., 2012). (Hensman et al., 2012) introduced a different view of CVB in a wide scope of exponential families. Onlya few researchers have examined CVB in models for structured data such as Probabilistic ContextFree Grammars (PCFG) (Wang and Blunsom, 2013b) and Hidden Markov Models (HMM) (Wangand Blunsom, 2013a). As stated, CVB solutions of IRM are first introduced to the literature to thebest of our knowledge.

3. Collapsed Gibbs Sampler Solution
Before deriving CVB solutions, we review the collapsed Gibbs sampler here to facilitate the derivation of CVB solutions. Let us define the counting statistics that are maintained during sampling.
nk,l and Nk,l denote the number of positive (x = 1) and negative (x = 0) relation observations in the(k, l)-cluster pairs, respectively. mk is the same quantity as used in Eqs. (1, 2).
XX
XX


nk,l =
z1,i,k z2, j,l xi, j , Nk,l =
z1,i,k z2, j,l 1 − xi, j ,(13)im1,k =

Xijz1,i,k ,im2,l =Xz2, j,l .

j

8

j

(14)

CVB Inference of IRM

3.1 Sampling z1,iWe review the collapsed Gibbs solution for inference of Z 1 = {z1,i }: the solution is completely
symmetric for the second domain, Z 2 .
In the Gibbs sampler approach, we repeat the following process. First, we select one object(1, i) (or (2, j)) from the data and take the object out from the model. More specifically, a clustering
assignment of the object z1,i is temporarily set empty (undefined). Then, we reassign (sample)
z1,i = k based on the posterior p z1,i = k .
To start, we divide the observations into two parts. Let us denote X(1,i) = {xi,· } as the set of allobservations concerning object i of the first domain. The remaining observations, hidden variables
excluding z1,i , and statistics computed on these data are denoted by \(1, i). Our target posterior isformulated as follows:



 

p z1,i = k|X, Z \(1,i), Z 2 ∝ p z1,i |Z 1\(1,i) p X|z1,i = k, Z 1\(1,i) , Z 2
1
 


Z 
\(1,i), Z 2 , X\(1,i) dΘ . (15)= p z1,i = k|Z 1
p X(1,i) |z1,i = k, Z 1\(1,i) , Z 2 , Θ p Θ|Z \(1,i)
1
The first term of the right-hand side of Eq. (15) becomes: \(1,i) existing clusters,
m1,k\(1,i)
p z1,i = k|Z 1∝α1
a new cluster,(16)

as in Eq. (2).
Thanks to the conjugacy of Eq. (6) and Eq. (9), we can easily evaluate the second integral termof the r.h.s. of Eq. (15). First, Eq. (17) is the posterior of the strength parameters Θ given all
information excepting the object i in the first domain.

 YY


\(1,i)\(1,i)\(1,i)
p Θ|Z \(1,i)
,
Z
,
X
=
Beta θk,l ; ak,l + nk,l, bk,l + Nk,l
,
(17)
2
1
k

l

 Y n+(1,i,k)
Nk,l+(1,i,k)
k,l
,
Z
,
Θ
=
θ
1
−
θ
.
p X(1,i) |z1,i = k, Z \(1,i)
2
k,l
1
k,l


(18)

l

n+(1,i,k) and N +(1,i,k) denote the statistics computed solely on X(1,i) given z1,i = k.
Combining Eq. (17) and Eq. (18), we obtain the following.
Z 
 

\(1,i)\(1,i)p X(1,i) |z1,i = k, Z \(1,i)
,
Z
,
Θ
p
Θ|Z
,
Z
,
X
dΘ
2
2
1
1


+(1,i,k)\(1,i)+(1,i,k)Y B ak,l + n\(1,i)
+
n
,
b
+
N
+
N
k,l
k,l
k,l
k,l
k,l


∝
, (19)\(1,i)\(1,i)
B
a
+
n
,
b
+
N
k,l
k,l
l
k,l
k,l
where B(·, ·) denotes the beta function.
Plugging Eq. (16) and Eq. (19) into Eq. (15) yields the posterior probability for sampling thecluster assignment of object i in the first domain:
\(1,i)+(1,i,k)\(1,i)+(1,i,k),bk,l +Nk,l+Nk,l\(1,i) Q B ak,l +nk,l  +nk,lmk is an existing cluster ,\(1,i)l B ak,l +n\(1,i),bk,l +Nk,l
 1,k
k,l


p z1,i = k|X, Z \(1,i),
Z∝\(1,i)+(1,i,k)2
1+n+(1,i,k),bk,l +Nk,l+Nk,lQ B ak,l +n\(1,i)

k,l
k,l




k is a new cluster .
α1 l\(1,i)\(1,i)
B ak,l +nk,l ,bk,l +Nk,l(20)
9Ishiguro, Sato, and Ueda

We iteratively take out z1,i from the statistics, compute the posterior and sample z1,i stochastically, then put z1,i back into the statistics.
3.2 Sampling z2, jWe also present the final result for the second domain. The derivation is symmetric to Z 1 ; thus, weomit details. Let us denote X(2, j) = {x·, j } as the set of all observations concerning the object j of
the second domain. The remaining observations, hidden variables, and statistics computed on thesedata are denoted by \(2, j).
Our goal is to sample z2, j based on the following equation:\(2, j)p z2, j = l|X, Z 1 , Z 2

Z 
 

\(2, j)\(2, j)\(2, j)∝ p z2, j = l|Z 2
p X(2, j) |z2, j = l, Z 1 , Z 2 , Θ p Θ|Z 1 , Z 2 , X\(2, j) dΘ . (21)Each term in the right-hand side of Eq. (21) is calculated as follows: \(2, j)
existing clusters,
m2,l\(2, j)p z2, j = l|Z 2∝α2
a new cluster,
Z

(22)


 

\(2, j)\(2, j)p X(2, j) |z2, j = l, Z 1 , Z 2 , Θ p Θ|Z 1 , Z 2 , X\(2, j) dΘj)+(2, j,l)\(2, j)+(2, j,l)
Y B ak,l + n\(2,+ nk,l , bk,l + Nk,l + Nk,lk,l. (23)∝

\(2, j)\(2, j)B ak,l + nk,l , bk,l + Nk,ll

In the above equation, n+(2, j,l) and N +(2, j,l) denote the statistics computed solely on X(2, j) given
z2, j = l. Thus, we obtain the Gibbs posterior probability for z2, j as follows:+(2, j,l)+(2, j,l)\(1,i)\(1,i),bk,l +Nk,l+Nk,l\(2, j) Q B ak,l +nk,l +nk,l




m
\(1,i)k 
B ak,l +n\(1,i),bk,l +Nk,l 2,l\(2, j)k,l
p z2, j = l|X, Z 1 , Z 2
∝
+(2,
j,l)
\(1,i)\(1,i)
Q B ak,l +nk,l +nk,l ,bk,l +Nk,l +Nk,l+(2, j,l)





α2 k
\(1,i)\(1,i)
B ak,l +nk,l ,bk,l +Nk,ll is an existing cluster ,l is a new cluster .

(24)
In practice, it is better to sample new assignments of objects in a domain-interleaving manner (z1,i → z2, j → z1,i0 → z2, j0 → . . . ) for faster convergence.
We can also sample hyperparameters and parameters by putting in hyper priors, or by solving
marginal likelihood maximization. We omit these procedures because they are out of our scope.
One difficulty in employing a Gibbs sampler is detection of convergence. A Gibbs sampler
assures asymptotic convergence to the true posteriors as the number of samples is infinitely many.
In practice, however, we will never have an infinite number of samples, so it is difficult to detect
convergence in a theoretically valid manner (Cowles and Carlin, 1996).
Another difficulty is the very slow mixing nature of collapsed Gibbs sampler on IRM, which
has been recently reported by (Albers et al., 2013). They showed that the several million iterations(sweeps) are not enough to mix the sampler, on 1000-nodes, real-world network data. One possible
reason is that one observed relation xi, j requires two hidden variables z1,i and z2, j , unlike topic
10

CVB Inference of IRM

models. To alleviate the slow mixing, we need to introduce much more sophisticated samplers suchas (Williamson et al., 2013). But such techniques would make it difficult to implement the sampler.
These two reasons motivate us to develop deterministic and fast VB-based inference solutions,though most of the existing IRM works rely on collapsed Gibbs sampler.

4. Variational Bayes Solution
No report to date has described a variational Bayes (VB) solution for IRM. However, it is beneficial
to quickly derive a VB solution for comparison with the proposed collapsed VB inference. In ageneral VB inference, we maximize the VB lower bound, which is defined as:

L =Z

q (Z 1 , Z 2 , Φ) log

p (X, Z 1 , Z 2 , Φ)dZ 1 dZ 2 dΦ ,q (Z 1 , Z 2 , φ)(25)

where Z 1 and Z 2 denote hidden variables, Φ denotes all associated parameters, X denotes all observations, and q(·) are the variational posteriors that approximate the true posteriors; all variational posteriors are assumed independent from each other. Maximizing the above lower bound isequivalent with minimizing the Kullback-Leibler divergence between the true posteriors p and thevariational posteriors q.
Generally speaking, the VB solution is analogous to the iterative process of the EM algorithm.
First, we maximize the VB lower bound w.r.t. the variational posteriors of hidden variables. Then,
we maximize the VB lower bound w.r.t. remaining parameters. This iteration monotonically increases the VB lower bound in Eq. (25); therefore, the VB solution halts automatically when it
reaches a local optimal point.

4.1 Generative models
We alter the generative model of IRM in two points. First, we use a Stick-Breaking Process (SBP)construction (Sethuraman, 1994) of the DPM (Eqs. (3, 4, 5)), instead of CRP (Eqs. (1, 2)). Second,
we “truncate” the maximum number of clusters, K, beforehand. Therefore, the VB solution of IRM
is doubly approximated: the independence of variational posteriors and the finite number of clusters.
Fixing the number of clusters seems to destroy the virtue of nonparametric Bayes: automatic model
selection. In practice, the SBP prior makes the unrepresented (unnecessary) clusters very small(very small weights) after inference. Therefore, it is easy to infer the true number of clusters evenif we “truncate” the infinite cluster representation.
11

Ishiguro, Sato, and Ueda

Here is the generative model of IRM for VB:
v1,k |α1 ∼ Beta (1, α1 ) ,π1,k = v1,k

(26)
K −1

k−1
Y

1
X

π1,m ,
1 − v1,m , π1,K1 = 1 −

(27)

m=1

m=1

z1,i |π1 ∼ Multinomial (π1 ) ,(28)

v2,l |α2 ∼ Beta (1, α2 ) ,(29)

π2,l = v2,l

l−1
Y

K −1

2
X

1 − v2,l , π2,K2 = 1 −π2,m ,

m=1

(30)

m=1

z2, j |π2 ∼ Multinomial (π2 ) ,θk,l |ak,l , bk,l ∼ Beta ak,l , bk,l ,xi, j |Z 1 , Z 2 , {θ} ∼ Bernoulli θz1,i ,z2, j .

(31)
(32)
(33)

The stick-breaking process is described in Eqs. (26-31). There are mainly two different parts
compared to the original generative models (Eq. (7) and Eq. (8)) formalized by CRP. The first one
is the introduction of cluster mixing proportional vectors π1 and π2 . The second one is the truncatedcluster numbers: K1 and K2 indicate the maximum number of truncated clusters for the first domain
and the second domain, respectively. Because the numbers of clusters are finite, we simply sample
the cluster assignment variables Z 1 and Z 2 from the multinomial distributions as in Eq. (28) andEq. (31).
4.2 Variational posteriors and the lower boundGiven the generative models, we introduce variational posteriors that are assumed independent fromeach other. Thanks to the model conjugacy, we can specify the form of the variational posteriors,
which are denoted by q(·) below:
K1


 Y

Beta v1,k ; α̂1,k , β̂1,k ,q v1 ; α̂1 , β̂1 =q (Z 1 ; π̂1 ) =

k=1
N1
Y

(34)


Multinomial z1,i ; π̂1,i ,

(35)



Beta v2,l ; α̂2,l , β̂2,l ,

(36)



Multinomial z2, j ; π̂2, j ,

(37)

K1 Y
K2


 Y

q Θ; â, b̂ =Beta θk,l ; âk,l , b̂k,l .

(38)



q v2 ; α̂2 , β̂2 =

i=1
K2
Y
l=1

q (Z 2 ; π̂2 ) =

N2
Y
j=1

k=1 l=1

12

CVB Inference of IRM

Following the definition, we obtain the following VB lower bound:"
log p (X, Z 1 , Z 2 , v1 , v2 , Θ)dZ 1 dZ 2 dv1 dv2 dΘdαL =q (Z 1 , Z 2 , v1 , v2 , Θ) log
log q (Z 1 , Z 2 , v1 , v2 , Θ)
= E Z 1 ,Z 2 ,Θ log p (X|Z 1 , Z 2 , Θ)
+ E Z 1 ,v1 log p (Z 1 |v1 )+ E Z 2 ,v2 log p (Z 2 |v2 )+ Ev1 log p (v1 |α1 )+ Ev2 log p (v2 |α2 )+ EΘ log p (Θ)
− E Z 1 log q (Z 1 )− E Z 2 log q (Z 2 )− Ev1 log q (v1 )− Ev2 log q (v2 )− EΘ log q (Θ) .

(39)
(40)
(41)
(42)
(43)
(44)
(45)
(46)
(47)
(48)
(49)

In the above equations, E x indicates the expectation of the predicate computed over the variational
posterior of x.
4.3 Variational posteriors of Z
We obtain the VB solution of the IRM by taking the derivative of the lower bound with respect tothe variational posterior parameters in Eqs. (34-38). Since the naive VB solution is not the primaryinterest of this paper, we omit the derivations and simply present the final results.
For the VB E-step, we compute the variational posteriors of hidden cluster assignment variables
Z 1 and Z 2 . The variational posterior parameters in Eq. (35) and Eq. (37) are shown below.
k−1 
k
X
 X


log π̂1,i,k = ψ α̂1,k +ψ β̂1,m −ψ α̂1,m + β̂1,m



m=1

−

N2 X
K2
X

m=1



π̂2, j,l ψ âk,l + b̂k,l

j=1 l=1

+

N2 X
K2
X

h
  i 
π̂2, j,l xi, j ψ âk,l + 1 − xi, j ψ b̂k,l + Const.

(50)

j=1 l=1
l−1
l
 X


 X 
log π̂2, j,l = ψ α̂2,l +ψ β̂2,m −ψ α̂2,m + β̂2,m
m=1

−
+

N1 X
K1
X
i=1 k=1
N1 X
K1
X

m=1



π̂1,i,k ψ âk,l + b̂k,lh  i 
π̂1,i,k xi, j ψ âk,l + 1 − xi, j ψ b̂k,l + Const.

i=1 k=1

13(51)

Ishiguro, Sato, and Ueda

In the above equations, ψ indicates the digamma function. π1,i,· and π2, j,· are to be normalized so as
to make the sum over cluster indices equal to 1.
4.4 Variational posteriors of parameters
In the VB-M step, we compute the posteriors of remaining parameters, Eq. (34), Eq. (36), andEq. (38). We again omit the derivations and present the final results.
α̂1,k = 1 +N1Xπ̂1,i,k ,(52)

π̂2, j,l ,

(53)

i=1

α̂2,l = 1 +

N2
X
j=1

β̂1,k = α1 +

N1 X
K1
X

π̂1,i,m ,(54)

π̂2, j,m ,(55)

π̂1,i,k π̂2, j,l xi, j ,(56)

π̂1,i,k π̂2, j,l (1 − xi, j ) .

(57)

i=1 m=k+1

β̂2,l = α2 +

N2 X
K2
X
j=1 m=l+1

âk,l = ak,l +

N1 X
N2
X
i=1 j=1b̂k,l = bk,l +

N1 X
N2
X
i=1 j=1

4.5 Optimizing hyperparameters
Optimization of the hyperparameters are also formulated as the maximization of the VB lower
bound.
It is easy to obtain update rules for hyperparameters by taking derivatives of the lower bound.
Employing the fixed-point methods presented in (Iwata et al., 2012; Minka, 2000; Wallach, 2008),
we have the following update rules for hyperparameters.
K1 i ,
α1 = PK h 
1
ψ
α̂
+
β̂
−
ψ
β̂1,k
1,k
1,k
k=1
K2

 i ,
α2 = PK h 2
l=1 ψ α̂2,l + β̂2,l − ψ β̂2,lψ ak,l + bk,l − ψ ak,l
ãk,l = ak,l ,ψ âk,l + b̂k,l − ψ âk,lψ ak,l + bk,l − ψ bk,l .
b̃k,l = bk,l ψ âk,l + b̂k,l − ψ b̂k,l
14

(58)
(59)
(60)
(61)

CVB Inference of IRM

5. Collapsed Variational Bayes (CVB) Solution of IRM
5.1 General Idea
The general idea of CVB inferences for hierarchical probabilistic models (Kurihara et al., 2007;
Teh et al., 2007, 2008; Asuncion et al., 2009; Sato and Nakagawa, 2012; Sato et al., 2012) is to
assume variational posteriors of hidden variables of the model where parameters are marginalizedout beforehand.
In Eq. (25), parameters Φ are not marginalized (collapsed) out in ordinary VB inference (Attias,
2000; Bishop, 2006). Thus, we need to compute the variational posteriors of the parameters as well.
The variational posteriors of the parameters impact the inference results, and this may increase thedanger of being trapped at a bad local optimal point.
CVB inference first marginalizes out the parameters in an exact way (as in a collapsed Gibbs
sampler). After that, remaining hidden variables are assumed to be independent from each other.
We can avoid the effect of parameter estimations and can reduce the number of quantities to beinferred because parameters are already marginalized. Further, it is known that the lower bound of
CVB is always tighter than that of the original VB (Teh et al., 2007). The formal definition of CVB
lower bound for IRM is:
L (Z 1 , Z 2 ) =Z

q (Z 1 , Z 2 ) log

p (X, Z 1 , Z 2 )dZ.
q (Z 1 , Z 2 )(62)

As evident, this is the same formulation as Eq. (25) except for the marginalized parameters.
To the best of our knowledge, this is the first work to formulate and derive CVB solutions for
IRM.
5.2 Generative model
Similar to VB, CVB also employs the truncated version of the IRM generative model. Let us
denote the truncated number of clusters of the first domain and the second domain as K1 and K2 ,respectively. For readers’ convenience, we present the SBP presentation of IRM again:
v1,k |α1 ∼ Beta (1, α1 ) ,(63)

v2,l |α2 ∼ Beta (1, α2 ) ,π1,k = v1,kπ2,l = v2,l

k−1
Y
m=1
l−1
Y


1 − v1,m , π1,K1 = 1 −
1 − v2,m , π2,K2 = 1 −

m=1

(64)
KX
1 −1
m=1
KX
2 −1

v1,m ,(65)

v2,m ,(66)

m=1

z1,i |π1 ∼ Multinomial (π1 ) ,(67)

z2, j |π2 ∼ Multinomial (π2 ) ,θk,l |ak,l , bk,l ∼ Beta ak,l , bk,l ,xi, j |Z 1 , Z 2 {θ} ∼ Bernoulli θz1,i ,z2, j .

(68)

15

(69)
(70)

Ishiguro, Sato, and Ueda

5.3 Counting statistics
The statistics required by the CVB solutions are defined in the same way as for the Gibbs samplers,
except for Eq. (72), which represents a kind of negative membership.

m1,k =N1
X

I(z1,i = k) =i=1

M1,k =N1X

z1,i,k ,m2,l =i=1

I(z1,i > k) =

K1
X

N2
X

N1 X
N2
X

z1,i,k z2, j,l xi, j ,I(z2, j = l) =j=1

m1,k0 ,

M2,l =

k0 =k+1i=1

nk,l =

N1
X

N2
X

N2
X

I(z2, j > l) =j=1

Nk,l =N1 X
N2X

i=1 j=1

z2, j,l ,

(71)

j=1
K2
X

m2,l0 ,

(72)

l0 =l+1z1,i,k z2, j,l (1 − xi, j ) .

(73)

i=1 j=1

5.4 Variational posterior of Z
Before deriving the variational posteriors, there are two points to note concerning the differencewith the original VB.
First, we assume the form of the variational posterior distributions beforehand in the case of VB.
Then, we directly compute the variational posterior parameters. This is possible because the wholegenerative model retains the conjugacy. In the case of CVB inference, however, the conjugacy does
not hold because we marginalize out parameters (v1 , v2 , Θ) in the inference. Therefore, we cannot
assume specific forms of variational posteriors, q (Z 1 ) and q (Z 2 ).
Second, in the case of naive VB, we write down the actual lower bound Eq. (25) for VB(Eqs.
(39-49). In the case of CVB, we do not explicitly write down the lower bound. The reason is thatconvenient forms of the lower bound are different for the hidden variables Z and the hyperparameters α, a, b. Both forms are equivalent, but in practice it is easier to choose different representations
to derive inference algorithms.
In fact, the procedure of CVB inference resembles the collapsed Gibbs samplers more thanordinary VB inferences. We take one object out from the model, recompute the posterior of theobject cluster assignment, and put the object back in the model. A difference is that CVB computes
soft cluster assignments of Z = {Z 1 , Z 2 } while the Gibbs sampler samples hard assignments foreach process. We repeat this process on all objects, and one iteration of CVB inference is done.
Let us derive the update rule of the hidden cluster assignment of the first domain z1,i . First,
we modify the representation of the CVB lower bound Eq. (62). The integral is replaced by thesummation because Z is discrete.
X X 
 \(1,i) p X(1,i) , X\(1,i) , z1,i , Z 1\(1,i) Z 2 \(1,i)
q z1,i q Z 1
L (z1,i , Z 1 , Z 2 ) =q (Z 2 ) log
  \(1,i) 
(Z
)
q
z
q
Z
q
\(1,i)1,i
2
z1,i Z
1
,Z 2
1

X
n



=
E Z \(1,i) ,Z 2 q z1,i log p X(1,i) |X\(1,i) , z1,i , Z \(1,i)
, Z2
1
z1,i

1




+ log p z1,i |Z \(1,i)− log q z1,i1o+ (Terms that are not related to z1,i ) .

16

(74)

CVB Inference of IRM

As in the case of the Gibbs sampler solution, X(1,i) = {xi,· } denotes the set of all observations
concerning object i of the first domain. The remaining observations, hidden variables excluding z1,i ,and statistics computed on these data are denoted by \(1, i). E x [y] indicates the expectation of y on
the variational posterior of x.
The above rewriting yields a Gibbs-like likelihood and prior terms in Eq. (15), averaged overother posteriors. Taking the derivative of Eq. (74) w.r.t. q(z1,i ) and equating it to zero, we have thefollowing update rule for q(z1,i ):

h

i
h

i

\(1,i)\(1,i)(1,i) \(1,i)
q z1,i ∝ exp E Z \(1,i) ,Z 2 log p X |X
, z1,i , Z 1 , Z 2 + E Z \(1,i) ,Z 2 log p z1,i |Z 1
. (75)
1

1

Before evaluating the expectations, we need to derive two terms inside E. Let us start from theprior part. Using Eq. (63) and Eq. (67), we readily obtain the following:p (Z 1 , v1 |α1 ) = α1K1

K1
Y



Beta v1,k ; m1,k + 1, M1,k + α1 B m1,k + 1, M1,k + α1 .

k=1

p (Z 1 |α1 ) =Z

p (Z 1 , v1 |α1 ) dv1 = α1K1



K1
Y
Γ m1,k + 1 Γ M1,k + α1.
Γ m1,k + M1,k + α1 + 1
k=1

(76)



Now we are ready to evaluate p z1,i = k|Z \(1,i).
1p (Z 1 |α1 )p Z \(1,i)
|α1
1



p z1,i = k|Z \(1,i)
,
α
=
1
1

m\(1,i)1,k + 1=\(1,i)
M1,k0 + α1k−1Y\(1,i)\(1,i)\(1,i)
m\(1,i)1,k + M1,k + α1 + 1 k0 =1 m1,k0 + M1,k0 + α1 + 1

.

(77)

The likelihood term is also easily available thanks to conjugacy:\(1,i)p X(1,i) |z1,i = k, Z \(1,i)
,
Z
,
X
2
1
 


 
\(1,i)+(1,i,k)\(1,i)+(1,i,k)\(1,i)K2
Y Γ ak,l + bk,l + n\(1,i)
Γ
a
+
n
+
n
Γ
b
+
N
+
N
+
N
k,l
k,l
k,l
k,l
k,l
k,l
k,l
k,l

 
 
 . (78)=\(1,i)\(1,i)\(1,i)\(1,i)+(1,i,k)+(1,i,k)
Γ
a
+
n
Γ
b
+
N
Γ
a
+
b
+
n
+
N
+
n
+
N
k,l
k,l
k,l
k,l
l=1
k,l
k,l
k,l
k,l
k,l
k,l
n+(1,i,k) and N +(1,i,k) denotes the statistics computed solely on X(1,i) given z1,i = k.
Plugging Eq. (77) and Eq. (78) into Eq. (75), then we obtain the variational posterior q(z1,i ).
q(z1,i ) must be normalized so that the summation for K1 clusters equals one.
In the same manner, we obtain the update rule for q(z2, j ).






h



q z2, j ∝ exp E Z ,Z \(2, j) log p X
1(2, j)|X\(2, j)2\(2, j), Z 1 , z2, j , Z 2\(2, j)p\(2, j)z2, j = l|Z 2 , α2



=

m2,l
\(2, j)m2,l+1\(2, j)+ M2,lii
\(2, j)+ E Z ,Z \(2, j) log p z2, j |Z 2
.
1 2
(79)
h

\(2, j)

l−1
Y

M2,l0
\(2, j)\(2, j)+ α2 + 1 l0 =1 m2,l0 + M2,l0

17

+ α2
+ α2 + 1.

(80)Ishiguro, Sato, and Ueda\(2, j)p X(2, j) |z2, j = l, Z 1 , Z 2 , X\(2, j)

 
 

K1 Γ a + b + n\(2, j) + N \(2, j) Γ a + n\(2, j) + n+(2, j,l) Γ b + N \(2, j) + N +(2, j,l)
Y
k,l
k,l
k,l
k,l
k,l
k,l
k,l
k,l
k,l
k,l
=

 
 . (81) \(2, j)\(2, j)\(2, j)+(2, j,l)+2 j jl\(2, j)
Γ
a
+
n
Γ
a
+
b
+
n
+
N
+
n
+
N
Γ
b
+
N
k,l
k,l
k,l
k,l
k=1
k,l
k,l
k,l
k,l
k,l
k,l
5.5 Computing Variational Expectations
The posterior of z1,i in Eq. (75) requires an expectation computation over the cluster assignmentsof Z 1\(1,i) and Z 2 . However, this is an intractable discrete combinatorial computation: there are
K1N1 −1 × K2N2 possible combinations. CVB inference approximates these expectations by Taylor
expansion. Let us denote the expectation of predicate x as a = E[x]. Then we have:f (x) ≈ f (a) + f 0 (a) (x − a) +1 00
f (a) (x − a)2 .
2(82)

Taking the expectations of both sides of Eq. (82), we obtain the following equation:

1
E[ f (x)] ≈ E[ f (a)] + E[ f 0 (a) (x − a)] + E[ f 00 (a) (x − a)2 ]2
100
2= f (a) + E[ f (a) (x − a) ]2
1
= f (E[x]) + f 00 (E[x]) V[x] .
2(83)

The 0th order term is constant. The 1st order term is canceled because x − a becomes zero by taking
the expectation. V denotes the posterior variance.
There are two types of approximations in CVB studies. The original CVB such as (Teh et al.,2007) employs 2nd-order Taylor approximation, and considers the variance as in Eq. (83). (Asuncion et al., 2009) reveals that the 0th-order Taylor approximation performs quite well in practice for
LDA. This is called the CVB0 solution, and it approximates the posterior expectation by

E[ f (x)] ≈ f (E[x]).

(84)

Obviously, the CVB0 solution is much simpler than that of the 2nd-order approximation, and is often
superior in terms of the perplexity of the learnt model (Asuncion et al., 2009; Sato and Nakagawa,
2012; Sato et al., 2012).
18

CVB Inference of IRM

Now we apply the Taylor approximation to two terms in the r.h.s. of Eq. (75). We pick one termfrom Eq. (78) and show how we can approximate the expectation.

h

i
h
i
\(1,i)E Z \(1,i) ,Z 2 log Γ ak,l + n\(1,i)
≈
log
Γ
E
a
+
n
\(1,i)k,lk,lk,lZ 1 ,Z 2
1
"

h
i
h
i2 #
1
\(1,i)\(1,i)\(1,i)+ Ψ E Z \(1,i) ,Z 2 ak,l + nk,l
E Z \(1,i) ,Z 2 ak,l + nk,l − E Z \(1,i) ,Z 2 ak,l + nk,l
1
1
1
2

h \(1,i) i
= log Γ ak,l + E Z \(1,i) ,Z 2 nk,l1
"h \(1,i) i
h \(1,i) i2 #1 \(1,i)E Z \(1,i) ,Z 2 nk,l − E Z \(1,i) ,Z 2 nk,l,
+ Ψ ak,l + E Z \(1,i) ,Z 2 nk,l
1
1
1
2

h \(1,i) i
= log Γ ak,l + E Z \(1,i) ,Z 2 nk,l1hih \(1,i) i
1 + Ψ ak,l + E Z \(1,i) ,Z 2 n\(1,i)
Vnk,l, (CVB)(85)\(1,i)k,l
Z 1 ,Z 2
1
2

h \(1,i) i
≈ log Γ ak,l + E Z \(1,i) ,Z 2 nk,l, (CVB0)(86)
1where Ψ is the trigamma function. For Eq. (77), the approximation is much simpler:hih \(1,i) i
E Z \(1,i) ,Z 2 m\(1,i)+ 1 (CVB,CVB0).
1,k + 1 = E Z \(1,i) ,Z 2 m1,k
1

(87)

1

From the above examples, we see that the expectation computations in Eq. (75) are achieved byreplacing counting statistics n, N, m, M with their expectations and variances.
The expectations and variances of the counting statistics in Eqs. (71-73) are computed as follows:

E[m1,k ] =N1Xq(z1,i,k ) ,E[m2,l ] =i=1

E[M1,k ] =

K1
X

N2
X

E[m1,k0 ] ,E[M2,l ] =

N1 X
N2
X

K2
X

N1 X
N2
X

E[m2,l0 ] ,(89)

l0 =l+1q(z1,i,k )q(z2, j,l )xi, j ,E[Nk,l ] =i=1 j=1

V[nk,l ] =

(88)

j=1

k0 =k+1

E[nk,l ] =q(z2, j,l ) ,

N1 X
N2
X

q(z1,i,k )q(z2, j,l )(1 − xi, j ) ,(90)

i=1 j=1q(z1,i,k ) 1 − q(z1,i,k ) q(z2, j,l ) 1 − q(z2, j,l ) xi,2 j ,

(91)




q(z1,i,k ) 1 − q(z1,i,k ) q(z2, j,l ) 1 − q(z2, j,l ) xi,2 j (1 − xi, j )2 .

(92)

i=1 j=1

V[Nk,l ] =

N1 X
N2
X
i=1 j=1

All expectations and variances are computed on q (Z 1 ) and q (Z 2 ).
19

Ishiguro, Sato, and Ueda

Based on Eqs. (88-92), we can easily derive the expectations and variances for the first domain
updates:
E[m\(1,i)
1,k ] =]=
E[n+(1,i,k)
k,l

X
i0 ,i
N2
X

\(1,i)]=
E[M1,k
q z1,i0 ,k = E[m1,k ] − q z1,i,k ,
q z2, j,l xi, j ,+(1,i,k)]=E[Nk,l

j=1

KX
k0 =k+1

N2
X



q z2, j,l 1 − xi, j ,

(93)

(94)

j=1


+(1,i,k)],E[n\(1,i)k,l ] = E[nk,l ] − q z1,i,k E[nk,l+(1,i,k)\(1,i)],] = E[Nk,l ] − q z1,i,k E[Nk,lE[Nk,l]=V[n+(1,i,k)k,l\(1,i)
E[m1,k
0 ]

(95)
(96)

N2
N2
X
X







 
2
+(1,i,k)2]=
q z2, j,l 1 − q z2, j,l xi, j , V[Nk,lq z2, j,l 1 − q z2, j,l 1 − xi, j ,
j=1

j=1

(97)


+(1,i,k)V[n\(1,i)],k,l ] = V[nk,l ] − q z1,i,k 1 − q z1,i,k V[nk,l
\(1,i)+(1,i,k)V[Nk,l ] = V[Nk,l ] − q z1,i,k 1 − q z1,i,k V[Nk,l
].

(98)
(99)



All expectations and variances are computed on q Z \(1,i)
and q (Z 2 ). By plugging Eqs. (93-99)
1
into Eq. (77) and Eq. (78), we can evaluate q(z1,i = k) for CVB and CVB0 solutions.
For the second domain updates, we have a completely symmetric story. Approximated expectations are:\(2, j)E[m2,l ] =X q z2, j0 ,l ,j0 , j+(2, j,l)E[nk,l ] =N1
X

q z1,i,k xi, j ,\(2, j)E[M2,lK2X]=

l0 =l+1+(2, j,l)E[Nk,l]=i=1\(2, j)E[m2,l0 ] ,

N1
X



q z1,i,k 1 − xi, j ,+(2, j,l)

]=

N1
X

(101)

i=1



\(2, j)+(2, j,l)E[nk,l ] = E[nk,l ] − q z2, j,l E[nk,l ] ,+(2, j,l)\(2, j)E[Nk,l ] = E[Nk,l ] − q z2, j,l E[Nk,l],V[nk,l

(100)

(102)
(103)
N

1
X
2



 
+(2, j,l)q z1,i,k 1 − q z1,i,k xi,2 j , V[Nk,l]=q z1,i,k 1 − q z1,i,k 1 − xi, j ,

i=1

i=1

(104)




\(2, j)+(2, j,l)V[nk,l ] = V[nk,l ] − q z2, j,l 1 − q z2, j,l V[nk,l ] ,





\(2, j)+(2, j,l)V[Nk,l ] = V[Nk,l ] − q z2, j,l 1 − q z2, j,l V[Nk,l
].

(105)
(106)

 \(2, j) All expectations and variances are computed on q (Z 1 ) and q Z 2. Plugging the above equations
into Eq. (80) and Eq. (81), we can then evaluate q(z2, j = l) for CVB and CVB0 solutions.
By iterating the variational posterior updates for all i, j, k, l in an interleaving manner, we obtainthe local optimal solutions of q(Z 1 ) and q(Z 2 ). The overall procedure for CVB inference of IRM is
described as pseudo codes in Fig. 2.
20

CVB Inference of IRM

input

Observation X = { xi, j } Hyperparameters α1, α2, { ak, l }, { bk, l }Maximum numbers of clusters K1, K2
Condition to convergenceprocess
initialize
initialize variational posteriors q(Z1 = { z1, i }), q(Z2 = { z2, j }) for all i, j, k, linitialize expected counting statistics in Eqs. (88-92) for all k, l
until satisfy convergence condition, iterate:
randomly select an object index from two domains
if object i in the 1st domain is selected:
compute expected statistics in Eqs. (93-99) for all k, lcompute expectations of Eqs. (77, 78) for all k
compute Eq. (75) for all knormalize Eq. (75) (sum to one) to obtain updated q(z1, i)
else if object j in the 2nd domain is selected:
compute expected statistics in Eqs. (100-106) for all k, lcompute expectations of Eqs. (80, 81) for all lcompute Eq. (79) for all lnormalize Eq. (79) (sum to one) to obtain updated q(z2, j)
end if-else
if necessary, update hyperparameters by computing Eqs. (108-111)
sort cluster indices k and l so as to align E[m1, k] and E[m2, l] in descending order

outputvariational posteriors q(Z1), q(Z2)Figure 2: Pseudo code of CVB inference for IRM.

21

Ishiguro, Sato, and Ueda

5.6 Optimizing hyperparameters
We can also derive CVB0-based update rules of hyperparameters by taking derivatives of Eq. (62)w.r.t. the hyperparameters. In particular, optimizing the concentration parameters of DP, α1 and α2are important. They play an important role in nonparametric Bayes modeling, but their update rule
for CVB has never been studied.
For the hyperparameters, we use another representation of the lower bound.
#X "p (X, Z 1 , Z 2 , α1 , α2 , a, b)
L (α1 , α2 , a, b) =q (Z 1 ) q (Z 2 ) log
q (Z 1 ) q (Z 1 )
Z 1 ,Z 2= E Z 1 ,Z 2 log p (X|Z 1 , Z 2 , a, b) + log p (Z 1 |α1 ) + log p (Z 2 |α2 ) − log q (Z 1 ) − log q (Z 2 )
"


#
XX
Γ ak,l + bk,l Γ ak,l + nk,l Γ bk,l + Nk,l=E Z 1 ,Z 2 logΓ ak,l Γ bk,l Γ ak,l + bk,l + nk,l + Nk,l
k
l
X



E Z 1 ,Z 2 log Γ M1,k + α1 − log Γ m1,k + M1,k + α1 + 1+ K1 log α1 ++ K2 log α2 +

k
X




E Z 2 ,Z 2 log Γ M2,l + α2 − log Γ m2,l + M2,l + α2 + 1

l





− E Z 1 ,Z 2 log q (Z 1 ) − E Z 1 ,Z 2 log q (Z 2 ) .

(107)

The last line is irrelevant to the hyperparameters.
To derive the update rules, we use the fixed point iteration technique used in (Iwata et al., 2012;
Minka, 2000; Wallach, 2008). For the concentration parameter α1 , we have the following inequality
between the current values of α1 and the new value α̂1 :log Γ M1,k + α̂1 − log Γ m1,k + M1,k + α̂1 + 1≥ log Γ M1,k + α1 − log Γ m1,k + M1,k + α1 + 1+ (α1 − α̂1 ) Ψ m1,k + M1,k + α1 + 1 − Ψ M1,k + α . (108)
We insert the above inequality into the lower bound Eq. (107), take derivative w.r.t. α̂1 with CVB0
approximation and equate to zero. Then, we obtain the following update rule for α1 :K1α̂1 = P 







 .
Ψ
E
m
+
E
M
Z 1 ,Z 2
1,k
Z 1 ,Z 2
1,k + α1 + 1 − Ψ E Z 1 ,Z 2 M1,k + α1
k

(109)

For α2 ,
K2







 .
α̂2 = P l Ψ E Z 1 ,Z 2 m2,l + E Z 1 ,Z 2 M2,l + α2 + 1 − Ψ E Z 1 ,Z 2 M2,l + α2

(110)

In the same manner, we have the update rules for the observation hyperparameters ak,l and bk,l .
 Ψ ak,l + E Z 1 ,Z 2 nk,l − Ψ ak,lâk,l = ak,l
 
 
.
(111)
Ψ ak,l + bk,l + E Z 1 ,Z 2 nk,l + E Z 1 ,Z 2 Nk,l − Ψ ak,l + bk,l Ψ bk,l + E Z 1 ,Z 2 Nk,l − Ψ bk,lb̂k,l = bk,l
 
 
.
(112)
Ψ ak,l + bk,l + E Z 1 ,Z 2 nk,l + E Z 1 ,Z 2 Nk,l − Ψ ak,l + bk,lAs evident, the update rules for CVB are very similar to those of VB (Eqs. (58-61)). The only
difference is that the CVB update rules incorporate the (approximated) counts, while the VB update
rules use the VB posterior parameters instead.
22

CVB Inference of IRM

6. Convergence Detection of CVB
It is theoretically guaranteed that each iteration of VB monotonically increases the variational lower
bound (Attias, 2000; Bishop, 2006), and VB eventually converges to its local optimal solutions
in finite iterations. Thus, the VB inference yields easy detection of convergence; the algorithm
automatically detects (technically sound) convergence by computing and monitoring the variational
lower bound (Eq. (25)) for each iteration.
Unfortunately, no theoretical guarantee of CVB convergence has been provided so far. Thisis due to the fact that we cannot correctly evaluate the posterior expectations over Z as we haveseen. What we try to find in the CVB solutions is a stationary point of a Taylor-approximated CVB
lower bound; thus, we are not sure that the procedure actually monotonically improves the true
lower bound. Convergence analysis of CVB inferences remains an important open problem in themachine learning field. However, the problem of CVB convergence is not so much discussed in theliterature since CVB inference yields better posterior estimations in many cases.
Convergence analysis of CVB remains an important but difficult problem. Instead of tackling
this problem directly, we study two aspects of CVB convergence in this paper. First, we empiricallystudy the convergence behaviors of CVB (in IRM) by monitoring a couple of quantities: a naive
VB lower bound and the pseudo leave-one-out log likelihood. We found that these two quantities
are potentially useful for CVB convergence detection. Second, we propose a simple and effective
technique to assure the automatic detection of CVB inference convergence called Averaged CVB(ACVB). We also prove that the ACVB reaches the stationary point of the true CVB lower bound,
if it exists.
From the non-expert user’s viewpoint, it is highly preferable if we can devise an easy convergence detection algorithm for CVB in general (not restricted to IRM). For many practitioners, it is
not easy to manually determine the convergence of the inference algorithms. This might be a part ofreasons why Maximum Likelihood estimators and EM-based algorithms are preferred. Therefore,
convergence-guaranteed ACVB would allow users to use CVB inference, which is more precise
than naive VB in theory, with automatic computation termination at guaranteed convergence.
To the best of our knowledge, (Foulds et al., 2013) is the only work that proposes convergenceassured CVB inference so far. This model is based on the Robbins and Monro stochastic approximation (Robbins and Monro, 1951) and is only valid for LDA-CVB0. More precisely, the solution
presented in (Foulds et al., 2013) is a MAP solution, leveraging the fact that the MAP solution isvery similar to the CVB0 solution in the case of LDA. On the contrary, ACVB proposed in thispaper is valid for any probabilistic model, and for both CVB and CVB0. The (Foulds et al., 2013)
approach is no longer valid for IRM because the MAP solution and the CVB0 solution are different.
6.1 Assessing Candidate Quantities for CVB Convergence Detection
We cannot correctly compute the true CVB lower bound in Eq. (62). Therefore, our first approach
would be to find some quantities that serve as the proxies of the true CVB lower bound.
For that purpose, we examine two quantities. The first candidate is the naive VB lower boundin Eq. (25). The VB lower bound is a lower bound of the true CVB lower bound. The variational
posteriors of parameters are computed using q(Z), which is obtained by the CVB inference. Using
these variational posteriors, we may compute the VB lower bound as a proxy of the true CVB lower
bound.
23

Ishiguro, Sato, and Ueda

VB Lowerbound & Pseudo LOO Log Likelihood (CVB, synth1)

VB Lowerbound & Pseudo LOO Log Likelihood (CVB, synth2)-6000-400000-8000-600000

Evaluation Quantity Value

Evaluation Quantity Value-10000-12000-14000-16000-18000-20000-22000-2400002040-800000-1000000-1200000-1400000

VB LowerboundPseudo LOO Log Likelihood-1600000

60-1800000
080

VB LowerboundPseudo LOO Log Likelihood

100

Iterations

20

40

60

80

100

Iterations

Figure 3: Evolutions of two quantities over CVB iterations. Solid lines indicate the evolutions ofnaive VB lower bound. Dashed lines indicate the evolution of pseudo leave-one-out log
likelihood on training data. Error bars denote the standard deviations. Left: computed on
Synth1 dataset. Right: computed on Synth2 dataset. For details of the datasets, see theexperiment section.

The second one is the pseudo leave-one-out log likelihood of the training data set. The CVB
solutions (and the Gibbs sampler) compute the predictive distribution of an object, say, z1,i , in aleave-one-out manner. Therefore, we might be able to detect the convergence of the CVB inference
by watching these predictive distributions. More precisely, we normalize Eq. (75) sum to one, inPcomputing q(z1,i ). The normalize constant, namely C1,i = k q z1,i,k serves as a pseudo leave-oneout log likelihood of the object (1, i) given the model. Then, the whole sum of this log likelihood,PP
C = i C1,i + j C2, j , is the pseudo log likelihood of the training data set.
Figure 3 and 4 present the evolutions of these two quantities in CVB (Fig. 3) and CVB0 (Fig. 4),
respectively. Hyperparameters are set to the best setup found in our experimental validations (seethe Experiment section). As evident from the cites, both of the quantities converge within a few
iterations.
It is interesting to see that the naive VB lower bound (presented in solid lines) increases as the
CVB inference proceeds. The CVB inference does not increase the VB lower bound directly, which
is a looser bound than that of CVB. However, the learned model actually decreases the discrepancy
between the variational posteriors and the true posteriors, resulting in increasing VB lower bound.
Unfortunately, the computational load is considerably heavy compared to the CVB updates. In
contrast, pseudo leave-one-out training log likelihood (presented in dashed lines) costs no extra
computational loads from the original CVB updates. Statistically, the property of the quantity ismore or less similar to the model evidence; thus the quantity would be a good choice for convergencedetection.
24

CVB Inference of IRM

VB Lowerbound & Pseudo LOO Log Likelihood (CVB0, synth2)

VB Lowerbound & Pseudo LOO Log Likelihood (CVB0, synth1)-400000

-6000
-8000

-600000

Evaluation Quantity Value

Evaluation Quantity Value-10000-12000-14000-16000-18000-20000-24000
0

20

40

60

80

-1000000-1200000-1400000-1600000

VB LowerboundPseudo LOO Log Likelihood-22000-800000-1800000
0100

VB LowerboundPseudo LOO Log Likelihood
20

40

60

80

100

Iterations

Iterations

Figure 4: Evolutions of two quantities over CVB0 iterations. Solid lines indicate the evolutions ofnaive VB lower bound. Dashed lines indicate the evolution of pseudo leave-one-out log
likelihood on training data. Error bars denote the standard deviations. Left: computed on
Synth1 dataset. Right: computed on Synth2 dataset. For details of the datasets, see theexperiment section.

6.2 Averaged CVB: Convergence technique for general CVB
Next, we propose a more direct and convergence-guaranteed technique for general CVB inferences.
The technique is based on monitoring the changes of q(Z). The rationale is simple: it is reasonable
to monitor q(Z) since the CVB solutions try to obtain the stationary point of the Taylor approximatedlower bound with respect to q(Z).
Then, we develop a simple annealing technique, called Averaged CVB (ACVB) to assure theconvergence of CVB solutions. We would like to emphasize that the discussion of ACVB is not
limited to the IRM: this technique is applicable to CVB inference on any model. Also, ACVB isvalid for CVB (2nd order) and CVB0 equally.
After a certain number of iterations for “burn-in”, we gradually decrease the portion of variational posterior changes using the following equation:
(s+1)

q̄

!
S
1
1 (s+1)
1 X (s)
(s)
(S )
= 1−
q̄ +
q
, or q̄ =
q ,
s+1
s+1
S s=1

(113)

where s denotes the iterations after completion of the “burn-in” period, q̄(s) denotes the “annealed”
variational posterior at the sth iteration, q(s) denotes the variational posterior by CVB inference atthe sth iteration, and S is the total number of iterations. After the “burn-in” period, we monitor theratio of changes of q̄ and detect convergence when the ratio falls below a predefined threshold. As
the final result, we do not finally use q(s) but q̄(s) . During the burn-in period, we monitor the changesof q: in most cases, q quickly converges before entering the annealing process.
Concerning the convergence of ACVB, there are two points to note. The first point is rather evident but makes the ACVB useful for practical CVB inference. ACVB updates assure convergence,
and we can easily detect the convergence by taking the difference of q̄ in successive iterations.
25

Ishiguro, Sato, and Ueda

Theorem 1 The averaged variational posterior q̄(s) is convergence-assured: ∀ > 0, ∃S 0 , s.t. ∀S >PN (S )−1)< .
q̄i − q̄(SS 0 ⇒ N1 i=1i
Proof Since!
S
S −1
1 X (s)1
1 X (s) 1 (S )q = 1−
q + q ,S s=1
S S − 1 s=1S
we haveS
S −1
S −1
1 1 X (s) 1 (S )1 X (s)1 X (s)
q −
q = −
q + q
S s=1
S − 1 s=1
S S − 1 s=1
S
S −1
1
1 1 X (s)≤|q | + |q(S ) |S S − 1 s=1
S

≤

1 1
1
2
(S − 1) + = .
S S −1
S
S

Thus,
N
S
S −1
1 X 1 X (s)2
1 X (s)qi −qi ≤ .
N i=1 S s=1
S − 1 s=1
S

If we set S 0 = 2 , then ∀S > S 0 ,
N
S
S −1
2
1 X 1 X (s)2
1 X (s)= .
qi −qi ≤ <
N i=1 S s=1
S − 1 s=1S
S0This meansN

1 X (S )−1)q̄ − q̄(S< .
i
N i=1 iThus, we can automatically stop the ACVB inference by using a stopping rule based on the difference of ACVB posteriors.
The second point is much noteworthy and validates the use of ACVB in Bayesian inference: wecan prove that the converged q̄ is asymptotically equivalent to the stationary point of the CVB lower
bound, if it exists (note that it is not clear whether the true CVB lower bound has a stationary pointin theory).
Theorem 2 If the variational posterior q(s) converges to a stationary point in the CVB lower bound,then the averaged variational posterior q̄(s) also converges to a stationary point in the CVB lower
bound.
26

CVB Inference of IRM

Proof Let q∗ be a stationary point in the CVB lower bound. By using this assumption,
lim q(s) = q∗ ⇔ ∀ > 0, ∃s0 s.t. ∀s > s0 ⇒ |q(s) − q∗ | < /2.

s→∞

Here, we defines0X(q(s) − q∗ ) = M > 0,s=1

and thus,
lim

s→∞

M
M
= 0 ⇔ ∀ > 0, ∃s00 s.t. ∀s > s00 ⇒< /2.
ss

When S 0 = max{s0 , s00 }, we have
∀S > S 0 , |q̄(S ) − q∗ | =

S
X
1

S
s=1

(q(s) − q∗ )

S
X
1 (s)
M
+
(q − q∗ )<S
S
s=S +10

S − S0≤ /2 +
/2 ≤ /2 + /2 = .
STherefore,
lim q̄(s) = q∗ .

s→∞

We want to stress that it remains unknown in the literature as to whether CVB inference has astationary point. However, we can still safely use ACVB because it assures convergence of theinference process and ACVB will find the ”true” solution if CVB has a stationary point. Suchsolutions for the convergence of CVB have been never studied, to the best of our knowledge.
Hereafter, we denote the (naive) CVB solution and the CVB0 solution, both with ACVB, as the
ACVB solution and the ACVB0 solution, respectively.

7. Speeding up CVB inferences of IRM
Convergence assurance by ACVB is beneficial for users because it enables automatic and easy
detection of inference convergence. However, the computational speed is also an important factor
for practical uses. In this section, we introduce two possible speed-up techniques for IRM-(A)CVB
solutions.
7.1 Cluster shrinkage
One drawback of CVB (and VB) inference is the computational cost per iteration. The Gibbs sampler dynamically shrinks and expands the cardinality of hidden clusters during inference. Thus, the27Ishiguro, Sato, and Ueda

computational cost of the Gibbs sampler is a function of the complexity of the hidden clusters. However, (C)VB solutions maintain all K clusters throughout the inference. Thus, there is no shrinkage
of clusters in CVB, which results in heavy computational costs, larger than the intrinsic complexityof the data.
A simple solution is to ignore small clusters. As described in (Kurihara et al., 2007), (C)VB
inferences need to maintain the order of clusters so that the memberships of clusters are aligned indescending order for better performance, in accordance with the stick-breaking process. Therefore,
it is easy to implement a heuristic that i) excludes small clusters from variational posterior updates,
and ii) avoids evaluating the contributions of small clusters in the updates.
k
For example, we can evaluate how small the cluster k is by Pml m. In practice, the small clusters
lk
become truly negligible in the sense of memberships: for example, Pml m∼ 1.0 × 10−5 or smaller.
lSetting a very small value for the threshold, the size of effective clusters automatically shrinks tothe size of the intrinsic data complexity. Then, the computational cost is reduced proportional to
the size of effective clusters. This dramatically speeds up the inference while barely harming theinference performance. In all experiments, we implemented this heuristic to eliminate unnecessary
computational costs for VB, CVB and CVB0 solutions.
7.2 Linear time inference for (A)CVB0
A naive implementation of IRM-CVB inference requires O (N1 N2 K1 K2 ) for one full sweep of hidand
den variable Z. We can see this instantly from the expectations E and variances V of n\(1,i)
k,l\(1,i)
Nk,l in Eqs. (94, 97). For update of q z1,i,k on specific i and k, we need to evaluate these expectations and evaluations for K2 times. This requires O (N2 K2 ) computations; thus the full sweep forN1 objects on K1 clusters requires O (N1 N2 K1 K2 ). The same holds for q(Z 2 ).
This prohibits applying IRM to larger data. However, for the (A)CVB0 solution, we can reduce
to O (L(N1 + N2 )K1 K2 ) where L denotes the average degree “1” links of objects, without any approximation. This is remarkable: we can solve IRM linear to the number of objects. Further, many
real-world relational data are very sparse: L is small. This makes (A)CVB0 even more efficient.
This is almost just an implementational issue, but we believe it is very beneficial for readers that areinterested in IRM for the first time.
To obtain this, we rewrite the variational expectation terms in Eq. (94) in the following way:+(1,i,k)E[nk,l
]=

N2
X
X 



q z2, j,l xi, j =q z2, j,l ,+(1,i,k)E[Nk,l
]=

(114)

j∈J +

j=1

N2
X


 X 

q z2, j,l 1 − xi, j =q z2, j,l .

(115)

j∈J −

j=1

n
o
n
o
where J + = j : xi, j = 1 and J + = j : xi, j = 0 . It is evident that J + ∪ J − = {1, 2, . . . , N2 }. The key
observation is:+(1,i,k)E[n+(1,i,k)] + E[Nk,l
]=
k,l

N2
X 
 X 
 X


q z2, j,l +q z2, j,l =q z2, j,l = E[m2,l ] .
j∈J +

j∈J −

(116)

j=1

The right-most term is a membership count of clusters in the second domain, defined in Eq. (88),which can be cached during inference. To compute Eq. (114), we only need to take L (the average
28

CVB Inference of IRM

degree) computations, which is much smaller than N2 for many real-world relational data. Combined with Eq. (116), we can evaluate Eq. (94), namely Eqs. (114, 115) by O (L), instead of N2 .
Thus, one full sweep of q(Z) computation is reduced to O (L(N1 + N2 )K1 K2 ).
This linear-time inference algorithm has two limitations. One is that this does not allow missing entries within X. More precisely, we can run the algorithm, but the resulting inference is notaccurate. This is because the above linear-time inference cannot correctly evaluate the existence of
missing data. If we encounter the relational data with missing entries, we need some preprocessing
to impute the missing entries. Another limitation is that this algorithm is not applicable to the CVB
solution: the posterior variation V[n] requires square terms of q(z), for which we cannot use thetrick of Eq. (116).
In the case of collapsed Gibbs, we can implement it in a similar way by replacing q(z) with I(z)on the current sample of Z. For VB, the update algorithm is very different from the others, but ingeneral, VB allows massive parallelization on the updates of q(z).
(Hansen et al., 2011; Albers et al., 2013), which employed collapsed Gibbs, presented a different
representation to avoid direct computations of negative counts (N in this paper). However, they did
not analyze the order of linear inference computations, nor did not emphasize the usefulness of thesparsity.

8. Experiments
In this section, we present the experimental validations. In summary, we confirmed the following
facts.
1. ACVB inferences achieved better modeling performances than the naive VB in large relational data sets. No significant differences in smaller relational data.
2. The ACVB0 solution is the fastest among deterministic inferences.
3. ACVB0 with linear time computation scales very well agasint large relational data.
8.1 Procedure
We compare the performance of proposed Averaged CVB solutions (ACVB, ACVB0) with a naive
variational Bayes (VB), which is a baseline deterministic inference. As a reference, we also include
comparisons with the collapsed Gibbs samplers (Gibbs) with very small number of iterations.
Initializations and hyperparameter choices are important for fair comparisons of inference methods. We employ hyperparameter updates for all solutions: fixed point iterations for VB, ACVB, and
ACVB0 and hyper-prior sampling for Gibbs. We test several initial hyperparameter values, andreport the results computed on the best hyperparameter setting. All hidden variables are initialized
in a completely random manner: we use the uniform distribution to assign soft values of p(zi = k).
In the case of Gibbs, we perform hard assignments of zi = k to the most weighted cluster. For VB,
ACVB, and ACVB0 solutions, we normalized the assigned weights sum to one over clusters. All
solutions without Gibbs require the number of truncated clusters a priori. To assess the effect oftruncation level, the experiments examined K1 = K2 = K ∈ {20, 40, 60}. In practice, we just need toprepare a sufficient number of K to explain data complexity.
Data modeling performance is evaluated by averaged test data marginal log likelihood. Given a
relational data matrix, we exclude a part of the relational observations (roughly 10% of matrix en29

Ishiguro, Sato, and Ueda

tries) from the inference as held-out test data. After the inference is finished, we compute marginal
log likelihoods of these test data. The number of test data, and chosen data entries are randomized
for every evaluation run. Thus, we average the log likelihoods by the number of test data entries. A
per-test data entry log likelihood is computed for 20 runs with different initializations and hyperparameter settings.
To compare the inference solutions in terms of the computational cost, we also monitored theconvergence behaviors and the computational time of the solutions. For the VB solution, we monitored the VB lower bound. For ACVB and ACVB0 solutions, we use the annealed posteriors. We
determined solutions’ convergence by the relative changes of the monitored quantity: if the changes
were smaller than 0.001% of the current value of the quantity, we assumed that the algorithm had
converged. As stated in the speeding-up section, we employed a cluster shrinkage technique for
VB, ACVB and ACVB0. We did not utilize the linear-time inference because of the existence oftest data that must be kept missing during inferences.
For the reference Gibbs sampler, we iterated the sampling procedure 3, 000 times; the first 1, 500
iterations were discarded as the burn-in period. As repeated explained, the collapsed Gibbs would
require tremendous amount of iterations (some millions) to obtain better modeling results (Albers
et al., 2013). However, we can not afford such computational resources for the collapsed Gibbs,
which has no easy way to detect convergence, not preferable for practitioners.
8.2 Datasets
We prepared several synthetic and real-world datasets for the experiments; they allow us to assess
the inferences in several scales and densities.
We generated two synthetic relation datasets. The size and true numbers of clusters of thesedatasets were: N1 = 100, N2 = 200, K1 = 4, K2 = 5 (synth 1), and N1 = 1000, N2 = 1500, K1 =7, K2 = 6 (synth 2).
The first real-world relational dataset is the Enron e-mail dataset (Klimt and Yang, 2004). Thisis a famous relational dataset used in many studies (Tang et al., 2008; Fu et al., 2009; Ishiguro
et al., 2010, 2012). We extracted monthly e-mail transactions for 2001. The dataset contained
N = N1 = N2 = 151 company members of Enron. xi, j = 1(0) if there is (not) an e-mail sent frommember i to member j. Out of twelve months, we selected the transactions of June (Enron Jun.),August (Enron Aug.), October (Enron Oct.), and December (Enron Dec.).
The second real-world relational dataset is the Lastfm dataset.1 This dataset contains several
records for the Last.fm music service, including lists of users’ most listened-to musicians, tag assignments for artists, and friend relations between users. We employ the friend relations between
N = N1 = N2 = 1892 users (Lastfm UserXUser). xi, j = 1(0) if there is (not) a friend relationfrom a user i to a user j. This dataset is 10 times larger than the Enron dataset in the number ofobjects, and 100 times larger in the number of matrix entries. We also employ the artist-tag relations between 17, 632 artists and 11, 946 tags. Since the relation matrix is too large for inferenceof Gibbs and naive VB, we truncate the number of artists and tags. The original observations are
the co-occurrence counts of (artist name, tag) pairs. We binarize the observations as to whether
the (artist name, tag) pair counts is greater than 1 or not: that is, we ignore one single occasional
co-occurrence of (artist name, tag). If the counts are greater than 1, then the observation entries are
set to 1; otherwise, set to 0. Then, all rows (artists) and columns (tags) that have no “1” entries are
1. provided by HetRec2011. http://ir.ii.uam.es/hetrec2011/30

CVB Inference of IRM

Table 3: Data sizes used in our experiments.
DatasetN1
N2# of “1” entries
synth1100
200
7608
synth21000 1500
578802Enron Jun.
151
151
257
Enron Aug.
151
151
439
Enron Oct.
151
151
707
Enron Dec.
151
151
377
Lastfm UserXUser
1892 1892
21512
Lastfm ArtistXTag
6099 1088
23253

Density
38.0%
38.6%
1.13%
1.93%
3.10%
1.66%
0.60%
0.35%

removed. The resulting binary matrix consists of N1 = 6099 artists and N2 = 1088 tags (Lastfm
ArtistXTag). xi, j = 1(0) if the artist i is (not) attached with the tag word j more than once.
The data sizes and densities are summarized in Table 3. The data sizes are rather small comparedto CVB research in LDA (Asuncion et al., 2009; Sato and Nakagawa, 2012; Sato et al., 2012). One
reason is that IRM deals with two hidden variables (z1,i , z2, j ) for one observation (xi, j ), while LDA
requires one hidden variable for one observation (word). This makes the inference difficult andhinders the scale up of the problem. In fact, existing IRM studies work in a somewhat similar
volume of datasets (Kemp et al., 2006; Ishiguro et al., 2010, 2012). Also, we cannot implement
linear time inference when we have missing data (test data) for model evaluations.
8.3 Results
8.3.1 Numerical Performance
The modeling performances of the solutions are presented in Table 4 (K = 20), Table 5 (K = 40),
and Table 6 (K = 60). They show the averages of test data marginal log likelihood after convergence.
Results of the best setup are presented for each solution. In addition, we conducted statisticalsignificance tests using t-tests.
These results reveal characteristics of the solutions in a few aspects.
First, ACVB inferences are significantly better than VB for larger datasets: synth2, and two
Lastfm datasets. Especially, we confirmed that ACVB0 always performed significantly better thanVB, and often recorded significantly better results than ACVB for those datasets. This indicates thatin potential ACVB inferences are superior to the naive VB inference as expected.
Second, we found no advantages of ACVB inferences over VB for smaller datasets: synth1
and Enron datasets. Specifically, the VB performed significantly better than ACVB solutions insynth1 data. The data is an artificial, dense and small cross-domain relation data. In such cases,
the VB still may finds good estimations of the true parameters θ. If so, VB may obtain better test
data log likelihood since ACVB marginalizes out all possibilities of the parameters, including ”bad”
estimations. Anyway, the synth1 data set is very small and dense. In general, we don’t face such datain our practical data analysis thus the results on larger and sparser data cases are more informative
for practical uses.
Third, the 3, 000 iterations of collapsed Gibbs sampler did not work well as we expected. Interestingly, in the case of Lastfm ArtixtXTag with K = 60, the Gibbs sampler performs significantly
31

Ishiguro, Sato, and Ueda

Table 4: Marginal test data log likelihood per test data entry (K = 20, 10% test data). Parenthesized
numbers indicate standard deviations. Larger values are better. Boldfaces indicate thebest method, which is significantly better than the method(s) marked with ∗ (by t-test,
p = 0.05).

DatasetSynth1Synth2Enron Jun.
Enron Aug.
Enron Oct.
Enron Dec.
Lastfm(UserXUser)
Lastfm(ArtistXTag)

Gibbs-0.3696* (0.0282)-0.3721* (0.0191)-0.0585 (0.0107)-0.0830* (0.0109)-0.1268* (0.0103)-0.0740 (0.0119)-0.0283* (0.0006)

VB-0.3260 (0.0155)-0.3737* (0.0090)-0.0547 (0.0087)-0.0789 (0.0076)-0.1164* (0.0091)-0.0693 (0.0068)-0.0287* (0.0005)

ACVB-0.3337 (0.0120)-0.3348* (0.0107)-0.0559 (0.0086)-0.0766 (0.0095)-0.1098 (0.0095)-0.0686 (0.0107)-0.0271* (0.0005)

ACVB0-0.3372* (0.0163)-0.3261 (0.0016)-0.0540 (0.0064)-0.0763 (0.0081)-0.1099 (0.0107)-0.0685 (0.0088)-0.0267 (0.0005)-0.0160* (0.0003)-0.0165* (0.0003)-0.0161* (0.0002)-0.0158 (0.0003)

Table 5: Marginal test data log likelihood per test data entry (K = 40, 10% test data). Parenthesized
numbers indicate standard deviations. Larger values are better. Boldfaces indicate the best
method which is significantly better than the method(s) marked with ∗ (by t-test, p = 0.05).

DatasetSynth1Synth2Enron Jun.
Enron Aug.
Enron Oct.
Enron Dec.
Lastfm(UserXUser)
Lastfm(ArtistXTag)

Gibbs-0.3657* (0.0224)-0.3712* (0.0122)-0.0595* (0.0110)-0.0838* (0.0088)-0.1256 (0.0111)-0.0750* (0.0095)-0.0280 (0.0008)*

VB-0.3246 (0.0141)-0.3743 (0.0108)*-0.0541 (0.0107)-0.0795 (0.0088)-0.1143 (0.0112)-0.0672 (0.0062)-0.0289* (0.0004)

ACVB-0.3430* (0.0197)-0.3285* (0.0050)-0.0531 (0.0068)-0.0770 (0.0084)-0.1141 (0.0125)-0.0688 (0.0099)-0.0272* (0.0005)

ACVB0-0.3380* (0.0146)-0.3254 (0.0015)-0.0558 (0.0070)-0.0766 (0.0072)-0.1145 (0.0115)-0.0678 (0.0101)-0.0267 (0.0004)-0.0161 (0.0003)-0.0167* (0.0004)-0.0162 (0.0003)-0.0162 (0.0003)

better than others. To explain this, we focus on the fact that the ACVB0 with K = 20 is significantly
better than the Gibbs with K = 60. It indicates that the data has much smaller complexity than we
expected. With greater K, the (C)VB inference algorithms may trapped at bad local optimum. Contrary, the collapsed Gibbs sampler obtained stable but not good solutions regardless of initial K. As
32

CVB Inference of IRM

Table 6: Marginal test data log likelihood per test data entry (K = 60, 10% test data). Parenthesized
numbers indicate standard deviations. Larger values are better. Boldfaces indicate thebest method, which is significantly better than the method(s) marked with ∗ (by t-test,
p = 0.05).

DatasetSynth1Synth2Enron Jun.
Enron Aug.
Enron Oct.
Enron Dec.
Lastfm(UserXUser)
Lastfm(ArtistXTag)

Gibbs-0.3668* (0.0155)-0.3624* (0.0190)-0.0573 (0.0087)-0.0862* (0.0102)-0.1281 (0.0148)-0.0794* (0.0120)-0.0283 (0.0006)*

VB-0.3281 (0.0127)-0.3736* (0.0090)-0.0569 (0.0069)-0.0772 (0.0098)-0.1162 (0.0100)-0.0682 (0.0098)-0.0287* (0.0005)

ACVB-0.3379* (0.0154)-0.3261 (0.0015)-0.0572 (0.0103)-0.0781 (0.0107)-0.1145 (0.0115)-0.0686 (0.0109)-0.0272 (0.0006)*

ACVB0-0.3452* (0.0113)-0.3258 (0.0015)-0.0563 (0.0064)-0.0754 (0.0105)-0.1139 (0.0096)-0.0690 (0.0124)-0.0267 (0.0006)-0.0160 (0.0003)-0.0167* (0.0003)-0.0163* (0.0002)-0.0163* (0.0003)

reported in (Albers et al., 2013), the collapsed Gibbs for IRM would require millions of iterations to
obtain better results. Thus it is perfectly possible that the collapsed Gibbs outperforms all VB-based
techniques provided the sophisticated sampling techniques and much more iterations.
Figure 5 and 6 present examples of obtained clustering for Synth2 and Lastfm UserXUser data
in K = 60. All object indices in the cites are sorted so that the objects are grouped into blocks inthe cites. Horizontal and vertical color lines indicates the borderlines of object clusters for the firstdomain i and the second domain j, respectively. We show the MAP assignments: we assign an
object into the cluster with the highest posterior probability.
8.3.2 Computational load and convergence behaviors
To assess the computational loads of four solutions, we have monitored CPU times for convergence.
First, we report the overall trends in convergence CPU times based on the average convergencetime presented in Table 7 (K = 20), Table 8 (K = 40), and Table 9 (K = 60). Aside from thecollapsed Gibbs, which has no no definite way to detect convergence of the inference, ACVB0
was magnitude-faster than the naive VB and the ACVB (2nd order) for almost all datasets. There
are several possible reasons. First, the update equations of (A)CVB0 is much simpler than thatof 2nd-order ACVB. Second, ACVB inference has fewer unknown variables to estimate than the
VB. Third, count statistics maintained in ACVB0 are able to efficiently cache and compute thanks
to their simplicity. Fourth, the landscape of ACVB0 posteriors may have smoother charactersitics
than those of the ACVB and the VB. Concerning the ACVB and the VB, the VB was faster whentruncated K is small. Also, the VB was faster than ACVB for dense synthetic data. In other cases,
the ACVB was faster than the VB.
33

Ishiguro, Sato, and Ueda

Figure 5: MAP clustering assignments of Synth2 dataset. All object indices are sorted.

34

CVB Inference of IRM

Figure 6: MAP clustering assignments of Lastfm UserXUser dataset. All object indices are sorted.

35

Ishiguro, Sato, and Ueda

Table 7: Average CPU time for convergence of the solutions (K = 20, 10% test data). Parenthesized
numbers indicate standard deviations. Computational time for convergence detection is
excluded.
DatasetSynth1Synth2Enron Jun.
Enron Aug.
Enron Oct.
Enron Dec.
Lastfm (User x User)
Lastfm (Artist X Tag)

Gibbs
72.90 (6.617)
3593.45 (22.17)
62.70 (17.33)
66.85 (5.480)
94.50 (7.710)
46.10 (1.136)
15238 (811.4)
35074 (1975)

VB
1.750 (1.894)
42.10 (4.571)
48.95 (37.14)
26.25 (32.64)
80.15 (38.70)
70.70 (43.65)
4627 (8006)
11237 (12183)

CVB
33.45 (23.10)
153.3 (69.34)
159.7 (92.97)
39.55 (5.408)
113.4 (54.49)
58.80 (26.45)
17668 (10427)
73186 (44579)

CVB0
1.550 (1.071)
13.30 (11.20)
3.400 (2.130)
1.250 (0.433)
3.400 (2.437)
3.800 (2.482)
400.5 (213.4)
1024 (358.9)

We also need to note that the CPU times for convergence are deeply affected by the convergencethreshold. In our experiments, we choose the threshold of 1.0 × 10−5 relative changes of the monitored quantities for VB, ACVB and ACVB0 solutions. If we change the threshold to 1.0 × 10−4 ,convergence times of these solutions become 10 times or more faster.
Finally, we show a few plots of test data likelihood evolutions over CPU times. Figure 7,
8, 9, and 10 respectively illustrate the time evolution of test data likelihood versus CPU time on
different datasets. End points of the plots indicate the average convergence time of (AC)VBs, or
3,000 iterations of collapsed Gibbs sampler. For all cases, we observe fast convergence evolutionsof ACVB0.
From these experimental results, we conclude that ACVB0 solutions are good for practitioners
who require good enough clustering results (possibly not a global optimum) with very fast computations and assured convergence.
8.4 Large data clustering experiment
To further demonstrate the usefulness of IRM with ACVB0, we conduct further experiments onclustering of a larger dataset.
Our scenario is a typical situation of practical relational data analysis. Our goal is to perform
clustering of relational data, hoping to extract some knowledge from the data. We do not need toevaluate the test generalization performance; thus we assume no missing entries within the relational
matrix X (or impute missing entries in preprocessing). Therefore, we can use a linear time ACVB0
inference. We cannot evaluate the data modeling performance by test data, thus we only show thecomputational time until convergence, with different N and different truncation level K.
We employ five relational data for clustering experiments. First, we borrow the two largest
datasets from the previous experiments: Lastfm UserXUser with the size of N1 = N2 = 1892
users and Lastfm ArtistXTag with N1 = 6, 099 artists by N2 = 1, 088 tags. The third data is theMovielens-10M dataset. The data consists of ratings on N1 = 10, 681 unique movies by N2 =69, 878 unique users. As the name indicates, there are about 10 million ratings. We treat all rated
36

CVB Inference of IRM

betterEnron Jun. K=40

Averaged marginal test data log likelihood

Gibbs
VB
ACVB
ACVB0

CPU time [sec]Figure 7: Averaged test data marginal log likelihoods vs. inference CPU time on Enron Jun. data,
K=40. The horizontal axis denotes CPU time [sec], and the vertical axis denotes average
test data marginal log likelihoods per relation entry. Presented Gibbs results are those ofsampled assignments, not of averaged posteriors.

betterLastfm UserXUser K=20

Averaged marginal test data log likelihood

Gibbs
VB
ACVB
ACVB0

CPU time [sec]Figure 8: Averaged test data marginal log likelihoods vs. inference CPU time on Lastfm UserXUser
data, K=20. The horizontal axis denotes CPU time [sec], and the vertical axis denotes
average test data marginal log likelihoods per relation entry. Presented Gibbs results arethose of sampled assignments, not of averaged posteriors.

37

Ishiguro, Sato, and Ueda

betterLastfm ArtistXTag K=20

Averaged marginal test data log likelihood

Gibbs
VB
ACVB
ACVB0

number of iterations

Figure 9: Averaged test data marginal log likelihoods vs. inference CPU time on Lastfm ArtistXTagdata, K=20. The horizontal axis denotes CPU time [sec], and the vertical axis denotes
average test data marginal log likelihoods per relation entry. Presented Gibbs results arethose of sampled assignments, not of averaged posteriors.

betterLastfm ArtistXTag K=60

Averaged marginal test data log likelihood

Gibbs
VB
ACVB
ACVB0

number of iterations

Figure 10: Averaged test data marginal log likelihoods vs. inference CPU time on Lastfm ArtistXTag data, K=60. The horizontal axis denotes CPU time [sec], and the vertical axis
denotes average test data marginal log likelihoods per relation entry. Presented Gibbs
results are those of sampled assignments, not of averaged posteriors.

38

CVB Inference of IRM

Table 8: Average CPU time for convergence of the solutions (K = 40, 10% test data). Parenthesized
numbers indicate standard deviations. Computational time for convergence detection is
excluded.
DatasetSynth1Synth2Enron Jun.
Enron Aug.
Enron Oct.
Enron Dec.
Lastfm (User x User)
Lastfm (Artist X Tag)

Gibbs
48.45 (1.396)
3970 (58.21)
73.50 (7.533)
66.20 (3.789)
69.05 (4.685)
46.85 (1.424)
15172 (592.7)
21509 (1131)

VB
3.750 (2.118)
140.2 (12.51)
75.15 (69.41)
85.45 (92.13)
134.2 (79.49)
103.8 (76.04)
19802 (18438)
14406 (12184)

CVB
74.60 (32.67)
474.0 (111.5)
39.40 (8.581)
30.05 (5.172)
36.60 (7.158)
25.10 (6.196)
13809 (5017)
69993 (13541)

CVB02.700 (1.269)
9.450 (9.227)
1.700 (0.557)
2.200 (0.400)
4.200 (1.470)
4.600 (1.685)
492.3 (238.7)
2271 (287.7)

Table 9: Average CPU time for convergence of the solutions (K = 60, 10% test data). Parenthesized
numbers indicate standard deviations. Computational time for convergence detection is
excluded.
DatasetSynth1Synth2Enron Jun.
Enron Aug.
Enron Oct.
Enron Dec.
Lastfm (User x User)
Lastfm (Artist X Tag)

Gibbs
72.20 (4.976)
4147 (65.96)
79.50 (24.52)
75.25 (7.203)
65.05 (6.087)
59.95 (5.643)
14730 (785.9)33704 (2391)

VB
9.800 (4.045)
318.7 (6.034)
202.1 (151.0)
287.4 (166.7)
291.3 (126.6)
225.6 (150.2)
21341 (26698)
33567 (32485)

CVB
62.75 (28.98)
908.9 (64.66)
58.45 (12.61)
42.45 (9.516)
89.60 (34.02)
27.40 (5.054)
17357 (6121)
64431 (8146)

CVB03.850 (1.236)
11.55 (3.694)
1.000 (0.000)
2.100 (0.740)
2.200 (0.400)
3.550 (0.9206)
450.0 (182.8)
1870 (178.4)

entries (regardless of the rating points) as positive relations. Also, we prepare the fourth and fifthlargest datasets from the Netflix data. The data consists of ratings by N1 = 480, 189 unique users on
N2 = 17, 770 unique movies. We can use the full dataset, but we prepare two subsets of the Netflix
data to measure the impact of the number of non-zero elements, which would affect the average
degrees. Netflix-rate1 data consists of rating entries that have “1” (worst) values. There are about
4 million “1” entries, and we treat them as the positive relation between users and movies. Netflixrate5 data consists of those with “5” (best) values. There are about 23 million “5” entries that are
assumed as positive relations.
Table 10 presents the CPU times for convergence. We tested on several hyperparameter setups,
and report the average CPU times of the setup of the best training data log likelihood. As evident
39

Ishiguro, Sato, and Ueda

Table 10: Data sizes and ACVB0 CPU times [sec] in large data clustering experiments.
DatasetLastfm U.XU.
Lastfm A.XT.
Movielens-10M
Netflix-rate1
Netflix-rate5

N1

N2

Density

1892
6099
10681
480189480189

1892
1088
69878
17770
17770

0.60%
0.35%
1.34%
0.05%
0.27%

CPU times (linear)
K=20
K=40
K=60
183.2
236.6
194.6
414.3
621.9
826.6
2087
9569
22466
24973 112116 126482
13440 104397 208604

CPU times (naive)
K=20 K=40 K=60
400.5 492.3 450.0
1024
2271
1870
NA
NA
NA
NA
NA
NA
NA
NA
NA

from the table, the linear inference of ACVB0 enables clustering computations on large relationaldata.
We can observe that the computational times are affected by several factors, but not as predictedfrom the theory. The computational order of IRM-ACVB0 is O(L × (N1 + N2 ) × K1 × K2 ). Thus,
computational time would grows linear to the number of objects and the density, and square to K.
In general, datasets with large N1 , N2 took more CPU time for convergence but CPU time is notproportional to N. For the effect of data density, please consult the rows of Netflix data. We seethat the density does not necessarily governs the CPU times for convergence. The CPU time does
not grow larger as expected from the model complexity as well. We expect to have four and ninetimes larger CPU times at K = 40 models and K = 60 models, compared to the K = 20 models.
However, this does not hold for all datasets excepting Movielens-10M. This is because the ACVB0
model shrinks as the model discovers lesser numbers of latent clusters from the given data, thanks
to the cluster shrinkage technique introduced in the Speeding-Up section.
Also note the convergence CPU time is deeply affected by the threshold of convergence detection: if we loosen the threshold from 0.001% relative changes to 0.01%, convergence typically
becomes 10 times (or more) faster.
We argue that the convergence-guaranteed ACVB0 is especially beneficial for large data analysis. We can solve collapsed Gibbs in linear time as well, but several millions of iterations are notenough to obtain good posterior estimations (Albers et al., 2013). Also the collapsed Gibbs requires
to monitor the inference process because we have no measure to detect convergence. In large data
analysis, this is costly and painful. In contrast, ACVB0 does not require such elaboration because
it can detect the assured convergence easily. Combining the test data modeling results and very fastcomputation times, the proposed ACVB0 solution is a good practical choice for IRM, even for large
relational data.

9. Conclusion
In this paper, we proposed Averaged collapsed variational Bayes (ACVB) inference of the Infinite
Relational Model (IRM), which is a convergence-guaranteed and practically useful deterministic
inference algorithm to replace naive VB.
First, we formulated a CVB lower bound for IRM based on the standard procedure, which is
intractable to evaluate exactly. For this problem, we used Taylor approximations as in CVB researchon topic models, and derived the full formulations and the inference procedure for two types of40

CVB Inference of IRM

CVB solutions. We also provided the CVB0-based update rules of hyperparameters, including theconcentration parameter of the Dirichlet Process, which has been never reported in the literature.
To make the CVB inference more practically useful, we studied the CVB inference in two
aspects. First is the convergence issue, which is an open problem for the CVB inference. We started
by examining two possible quantities to assess the convergence of CVB solutions of IRM. After
that, we proposed a simple and effective annealing technique, Averaged CVB (ACVB), to assure
the convergence of CVB solutions. ACVB posterior update offers assured convergence thanks to its
simple annealing mechanism. Moreover, the stationary point of the CVB lower bound is equivalent
to the converged solution of ACVB, if the lower bound has a stationary point (an issue unresolvedin the literature). ACVB is applicable to any model, and is equally valid for CVB and CVB0.
The second aspect is the computational speed of CVB. We proposed a cluster shrinkage technique and a linear-time inference implementation. These techniques make the IRM inference morescalable against the data size, and open the door to larger and more complex relational data analysis
applications.
The resulting CVB solutions offer more precise inference than naive VB in experiments. At thesame time, the annealing ACVB technique allows us to automatically detect convergence and yields
short computational time. We also confirmed that the linear time inference of (A)CVB0 allows us
to analyze large two-place relational data.
As future work, we will further enhance inference speed. One possible solution is to stochastically approximate the sample size as in SGD. Recently, (Foulds et al., 2013) proposed such approximation for LDA. Another way is to parallelize the inference procedure, as (Hansen et al., 2011;
Albers et al., 2013) have examined the parallelization of collapsed Gibbs samplers on IRM. It isalso important to explore efficient CVB algorithms for more advanced models such as MMSB andits followers (Airoldi et al., 2008; Miller et al., 2009; Griffiths and Ghahramani, 2011). Aside fromthe representation of multiple cluster assignments, a few studies have headed toward to other issues.
For example, (Fu et al., 2009; Ishiguro et al., 2010) focused on dynamics of network evolution in thecontext of stochastic blockmodels (MMSB and IRM). Subset IRM (Ishiguro et al., 2012) is anotherextension of IRM that automatically ”filters out” nodes from the clustering that are not so informative to group. Applying CVB for these models may make it easier for practitioners to examine thedepth of various relational data.

References
Edoardo M Airoldi, David M Blei, Stephen E Fienberg, and Eric P Xing. Mixed Membership
Stochastic Blockmodels. Journal of Machine Learning Research, 9:1981–2014, 2008.
Kristoffer Jon Albers, Andreas Leon Aagard Moth, Morten Mø rup, and Mikkel N . Schmidt. Large
Scale Inference in the Infinite Relational Model: Gibbs Sampling is not Enough. In Proceedingsof the IEEE International Workshop on Machine Learning for Signal Processing (MLSP), 2013.
Arthur Asuncion, Max Welling, Padhraic Smyth, and Yee Whye Teh. On Smoothing and Inference
for Topic Models. In Proceedings of the 25th Conference on Uncertainty in Artificial Intelligence(UAI), 2009.
Hagai Attias. A Variational Bayesian Framework for Graphical Models. In Advances in Neural
Information Processing Systems 12 (Proceedings of NIPS), pages 209–215, 2000.
41

Ishiguro, Sato, and Ueda

Christopher M. Bishop. Pattern Recognition and Machine Learning. Springer-Verlag New York,2006.
D. Blackwell and J. B. MacQueen. Ferguson Distributions via Polya urn schemes. The Annals of
Statistics, 1(2):353–355, 1973.
David M. Blei, Thomas L. Griffiths, and Michael I. Jordan. The nested chinese restaurant process
and bayesian nonparametric inference of topic hierarchies. Journal of the ACM, 57(2):7:1–30,
2010.
A. Clauset, C. Moore, and M. E. J. Newman. Hierarchical Structure and the Prediction of Missing
Links in Networks. Nature, 453:98–101, 2008.
Mary Kathryn Cowles and Bradley P. Carlin. Markov Chain Monte Carlo Convergence Diagnostics: A Comparative Review. Journal of the American Statistical Association (JASA), 91(434):883–
904, 1996.
E. Erosheva, S. Fienberg, and J. Lafferty. Mixed-membership Models of Scientific Publications.
Proceedings of the National Academy of Sciences of the United States of America (PNAS), 101(Suppl 1):5220–5227, 2004.
James Foulds, Levi Boyles, Christopher DuBois, Padhraic Symyth, and Max Welling. Stochastic
collapsed variational Bayesian inference for latent Dirichlet allocation. In Proceedings of the
19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD),
2013.
Wenjie Fu, Le Song, and Eric P. Xing. Dynamic mixed membership blockmodel for evolving
networks. In Proceedings of the 26th Annual International Conference on Machine Learning(ICML), 2009. ISBN 9781605585161.
Thomas L. Griffiths and Zoubin Ghahramani. The Indian Buffet Process : An Introduction andReview. Journal of Machine Learning Research, 12:1185–1224, 2011.
Toke Jansen Hansen, Morten Mø rup, and Lars Kai Hanse. Non-parametric Co-clustering of Large
Scale Sparse Bipartite Networks on the GPU. In Proceedings of the IEEE International Workshop
on Machine Learning for Signal Processing (MLSP), 2011.
James Hensman, Magnus Rattray, and Neil D Lawrence. Fast Variational Inference in the Conjugate
Exponential Family. arXive, page 1206.5162v2, 2012.
Qirong Ho, Ankur P Parikh, and Eric P Xing. Multiscale Community Blockmodel for Network
Exploration. In Proceedings of the 14th International Conference on Artificial Intelligence andStatistics (AISTATS), 2011.
Qirong Ho, Junming Yin, and Eric P. Xing. On Triangular versus Edge Representations - Towards
Scalable Modeling of Networks. In Advances in Neural Information Processing Systems 25(Proceedings of NIPS), 2012.
Katsuhiko Ishiguro, Tomoharu Iwata, Naonori Ueda, and Joshua Tenenbaum. Dynamic Infinite
Relational Model for Time-varying Relational Data Analysis. In Advances in Neural Information
Processing Systems 23 (Proceedings of NIPS), 2010.
42

CVB Inference of IRM

Katushiko Ishiguro, Naonori Ueda, and Hiroshi Sawada. Subset Infinite Relational Models. In Proceedings of the 15th International Conference on Artificial Intelligence and Statistics (AISTATS),2012.
Tomoharu Iwata, Takeshi Yamada, Yasushi Sakurai, and Naonori Ueda. Sequential Modeling ofTopic Dynamics with Multiple Timescales. ACM Transactions on Knowledge Discovery from
Data, 5(4):19:1–19:27, 2012.
Charles Kemp, Joshua B Tenenbaum, Thomas L Griffiths, Takeshi Yamada, and Naonori Ueda.
Learning Systems of Concepts with an Infinite Relational Model. In Proceedings of the 21st
National Conference on Artificial Intelligence (AAAI), 2006.
Bryan Klimt and Yiming Yang. The Enron Corpus : A New Dataset for Email Classification Research. In Proceedings of the European Conference on Machine Learning (ECML), 2004.
Kenichi Kurihara, Max Welling, and Yee Whye Teh. Collapsed Variational Dirichlet Process Mixture Models. In Proceedings of the 20th International Joint Conference on Artificial Intelligence(IJCAI), 2007.
David Liben-Nowell and Jon Kleinberg. The Link Prediction Problem for Social Networks. In
Proceedings of the Twelfth Annual ACM International Conference on Information and Knowledge
Management (CIKM), 2003.
Kurt T Miller, Thomas L Griffiths, and Michael I Jordan. Nonparametric Latent Feature Models
for Link Prediction. In Advances in Neural Information Processing Systems 22 (Proceedings of
NIPS), 2009.
Thomas P. Minka.
Estimating a Dirichlet distribution.
us/um/people/minka/papers/dirichlet/, 2000.

http://research.microsoft.com/en-

Morten Mø rup, Kristoffer Hougaard Madsen, Anne Marie Dogonowski, Hartwig Siebner, andLars Kai Hansen. Infinite Relational Modeling of Functional Connectivity in Resting State fMRI.
In Advances in Neural Information Processing Systems 23 (Proceedings of NIPS), 2010.
Konstantina Palla, David A Knowles, and Zoubin Ghahramani. An Infinite Latent Attribute Modelfor Network Data. In Proceedings of the 29th International Conference on Machine Learning(ICML), 2012.
H. Robbins and S. Monro. A stochastic approximation method. The Annals of Mathematical Statistics, pages 400–407, 1951.
Issei Sato and Hiroshi Nakagawa. Rethinking Collapsed Variational Bayes Inference for LDA. In
Proceedings of the 29th International Conference on Machine Learning (ICML), 2012.
Issei Sato, Kenichi Kurihara, and Hiroshi Nakagawa. Practical collapsed variational bayes inference for hierarchical dirichlet process. In Proceedings of the 18th ACM SIGKDD international
conference on Knowledge discovery and data mining (KDD), 2012.
J. Sethuraman. A Constructive Definition of Dirichlet Process. Statistica Sinica, 4:639–650, 1994.
43

Ishiguro, Sato, and Ueda

Lei Tang, Huan Liu, Jianping Zhang, and Zohreh Nazeri. Community evolution in dynamic multimode networks. In Proceeding of the 14th ACM SIGKDD international conference on Knowledge
discovery and data mining (KDD), 2008.
Yee Whye Teh, David Newman, Max Welling, and D Neaman. A Collapsed Variational BayesianInference Algorithm for Latent Dirichlet Allocation. In Advances in Neural Information Processing Systems 19 (Proceedings of NIPS), 2007.
Yee Whye Teh, Kenichi Kurihara, and Max Welling. Collapsed Variational Inference for HDP. In
Advances in Neural Information Processing Systems 20 (Proceedings of NIPS), 2008.
Hanna M Wallach. Structured Topic Models for Language. PhD thesis, University of Cambridge,
2008.
Pengyu Wang and Phil Blunsom. Collapsed Variational Bayesian Inference for Hidden Markov
Models. In Proceedings of the 16th International Conference on Artificial Intelligence and Statistics (AISTATS), 2013a.
Pengyu Wang and Phil Blunsom. Collapsed Variational Bayesian Inference for PCFGs. In Proceedings of the Seventeenth Conference on Computational Natural Language Learning (ACL),
2013b.
Sinead A Williamson, Avinava Dubey, and Eric P Xing. Parallel Markov Chain Monte Carlo for
Nonparametric Mixture Models. In Proceedings of the 30th International Conference on Machine
Learning (ICML), volume 28, 2013.
Jaewon Yang and Jure Leskovec. Overlapping Community Detection at Scale : A Nonnegative
Matrix Factorization Approach. In Proceedings of the 6th ACM International Conference on Web
Search and Data Mining (WSDM), 2013.
Jaewon Yang, Julian Mcauley, and Jure Leskovec. Community Detection in Networks with Node
Attributes. In Proceedings of the IEEE Conference on Data Mining (ICDM), 2013.
Junming Yin, Qirong Ho, and Eric P Xing. A Scalable Approach to Probabilistic Latent Space
Inference of Large-Scale Networks. In Advances in Neural Information Processing Systems 26(Proceedings of NIPS), 2013.
Shenghuo Zhu, Kai Yu, and Yihong Gong. Stochastic Relational Models for Large-scale Dyadic
Data using MCMC. In Advances in Neural Information Processing Systems 21 (Proceedings of
NIPS), 2009.

44