Bayesian Mean-parameterized Nonnegative Binary Matrix Factorization

TL;DR

This paper introduces a Bayesian mean-parameterized nonnegative binary matrix factorization framework that avoids link functions, includes novel inference algorithms, and automatically determines the number of latent factors, improving interpretability and performance.

Contribution

It proposes a unified Bayesian mean-parameterized NMF model for binary data, with new inference algorithms and a nonparametric extension for automatic latent dimension selection.

Findings

Achieves comparable or better results than existing methods.

Automatically detects the number of relevant components.

Effective in tasks like dictionary learning and missing data prediction.

Abstract

Binary data matrices can represent many types of data such as social networks, votes, or gene expression. In some cases, the analysis of binary matrices can be tackled with nonnegative matrix factorization (NMF), where the observed data matrix is approximated by the product of two smaller nonnegative matrices. In this context, probabilistic NMF assumes a generative model where the data is usually Bernoulli-distributed. Often, a link function is used to map the factorization to the $[0,1]$ range, ensuring a valid Bernoulli mean parameter. However, link functions have the potential disadvantage to lead to uninterpretable models. Mean-parameterized NMF, on the contrary, overcomes this problem. We propose a unified framework for Bayesian mean-parameterized nonnegative binary matrix factorization models (NBMF). We analyze three models which correspond to three possible constraints that respect the mean-parametrization without the need for link functions. Furthermore, we derive a novel collapsed Gibbs sampler and a collapsed variational algorithm to infer the posterior distribution of the factors. Next, we extend the proposed models to a nonparametric setting where the number of used latent dimensions is automatically driven by the observed data. We analyze the performance of our NBMF methods in multiple datasets for different tasks such as dictionary learning and prediction of missing data. Experiments show that our methods provide similar or superior results than the state of the art, while automatically detecting the number of relevant components.

Tables

Table 1. Table 1 : Bernoulli matrix factorization methods considered in the literature. Bernoulli, Gamma, Dirichlet distributions are denoted as Ber, Ga, and Dir, respectively. Gradient refers to Gradient-based optimization.

Reference	Likelihood	Prior / Constr.	Estimation
Sammel et al. (1997)	$\text{Ber}\left(\sigma \left({\left[\mathrm{𝐖𝐇}\right]}_{fn}\right)\right)$	$w_{fk}\sim \text{Normal}(\cdot )$	MC-EM
		$h_{kn}\in ℝ$
Tipping (1999)	$\text{Ber}\left(\sigma \left({\left[\mathrm{𝐖𝐇}\right]}_{fn}\right)\right)$	$w_{fk}\sim \text{Normal}(\cdot )$	vEM
		$h_{kn}\in ℝ$
Collins et al. (2002)	$\text{Ber}\left(\sigma \left({\left[\mathrm{𝐖𝐇}\right]}_{fn}\right)\right)$	$w_{fk}\in ℝ$	Gradient
		$h_{kn}\in ℝ$
Schein et al. (2003)	$\text{Ber}\left(\sigma \left({\left[\mathrm{𝐖𝐇}\right]}_{fn}\right)\right)$	$w_{fk}\in ℝ$	ALS
		$h_{kn}\in ℝ$
Meeds et al. (2007)	$\text{Ber}\left(\sigma \left({\left[\mathrm{𝐖𝐗𝐇}\right]}_{fn}\right)\right)$	$w_{fk}\sim \text{Ber}(\cdot )$	Gibbs
		$h_{kn}\sim \text{Ber}(\cdot )$
Kabán and Bingham (2008)	$\text{Ber}\left({\left[\mathrm{𝐖𝐇}\right]}_{fn}\right)$	$w_{kn}\sim \text{Beta}(\cdot )$	VB
		${𝐡}_n\sim \text{Dir}(\cdot )$
Bingham et al. (2009)	$\text{Ber}\left({\left[\mathrm{𝐖𝐇}\right]}_{fn}\right)$	${\sum }_k{w}_{fk}=1$	EM
		$h_{kn}\in [0,1]$
Tomé et al. (2013)	$\text{Ber}\left(\sigma \left({\left[\mathrm{𝐖𝐇}\right]}_{fn}\right)\right)$	$w_{fk}\in {ℝ}_{+}$	Gradient
		$h_{kn}\in ℝ$
Larsen and Clemmensen (2015)	$\text{Ber}\left(\sigma \left({\left[\mathrm{𝐖𝐇}\right]}_{fn}\right)\right)$	$w_{fk}\in {ℝ}_{+}$	Gradient
		$h_{kn}\in {ℝ}_{+}$
Zhou (2015)	$\text{Ber}\left(f\left({\left[\mathrm{𝐖𝐇}\right]}_{fn}\right)\right)$	$w_{fk}\sim \text{Ga}(\cdot )$	Gibbs
		$h_{kn}\sim \text{Ga}(\cdot )$

Equations

$$\mathbf{V}\approx \mathbf{WH},$$

$$\mathbf{W},\mathbf{H}=\underset{\mathbf{W},\mathbf{H}}{\mathrm{argmin}}D(\mathbf{V}|\mathbf{WH})\text{s.t}\mathbf{W}\ge 0,\mathbf{H}\ge 0,$$

$$\underset{\mathbf{W},\mathbf{H}}{\mathrm{argmin}}D(\mathbf{V}|\mathbf{WH})=\underset{\mathbf{W},\mathbf{H}}{\mathrm{argmax}}p(\mathbf{V}|\mathbf{W},\mathbf{H}),$$

$$p(\mathbf{W},\mathbf{H}|\mathbf{V})=p(\mathbf{V}|\mathbf{W},\mathbf{H})p(\mathbf{W},\mathbf{H})/p(\mathbf{V}).$$

$$p(\mathbf{V}|\mathbf{W},\mathbf{H})=\prod _{f,n}\text{Bernoulli}\left(v_{fn}|\varphi ([\mathbf{WH}{]}_{fn})\right),$$

$$p(\mathbf{V}|\mathbf{W},\mathbf{H})=\prod _{f,n}\text{Bernoulli}(v_{fn}|[\mathbf{WH}{]}_{fn}),$$

$$v_{fn}\sim \text{Bernoulli}\left(\sigma \left([\mathbf{W}\mathbf{H}{]}_{fn}\right)\right),$$

$$v_{fn}\sim \text{Poisson}([\mathbf{W}\mathbf{H}{]}_{fn}).$$

$$v_{fn}\sim \text{Bernoulli}\left(f\left([\mathbf{W}\mathbf{H}{]}_{fn}\right)\right).$$

$$c_{fn}$$

$$v_{fn}$$

$$v_{fn}\sim \text{Bernoulli}([\mathbf{WH}{]}_{fn}).$$

$$\begin{array}{rl}& \text{Beta-Dir}\text{ (c1)}\\ & h_{kn}\sim \text{Beta}({\alpha }_k,{\beta }_k)\\ & {\mathbf{w}}_f\sim \text{Dirichlet}(\mathbf{\gamma })\end{array}\begin{array}{rl}& \text{Dir-Beta}\text{ (c2)}\\ & {\mathbf{h}}_n\sim \text{Dirichlet}(\eta )\\ & w_{fk}\sim \text{Beta}({\alpha }_k,{\beta }_k)\end{array}\begin{array}{rl}& \text{Dir-Dir}\text{ (c3)}\\ & {\mathbf{h}}_n\sim \text{Dirichlet}(\mathbf{\eta })\\ & {\mathbf{w}}_f\sim \text{Dirichlet}(\mathbf{\gamma })\end{array}$$

$${\mathbf{h}}_n$$

$${\mathbf{v}}_n|{\mathbf{h}}_n$$

$$\underset{\mathbf{W},\mathbf{H}}{\min }-\log p(\mathbf{V}|\mathbf{WH})=\sum _{f,n}{v}_{fn}\log ([\mathbf{W}\mathbf{H}{]}_{fn})+(1-v_{fn})\log (1-[\mathbf{W}\mathbf{H}{]}_{fn})$$

$$\text{subject to}(1)w_{fk},h_{kn}\ge 0,(2)\sum _k{w}_{fk}=1,(3)h_{kn}\in [0,1].$$

$$h_{kn}$$

$${\mathbf{w}}_f$$

$${\mathbf{z}}_{fn}|{\mathbf{w}}_f$$

$$v_{fn}|{\mathbf{h}}_n,{\mathbf{z}}_{fn}$$

$$p(v_{fn}|{\mathbf{w}}_f,{\mathbf{h}}_n)$$

$$=\sum _k{w}_{fk}{h}_{kn}^{v_{fn}}(1-h_{kn}{)}^{1-v_{fn}},$$

$$p(\mathbf{V},\mathbf{Z},\mathbf{W},\mathbf{H})=p(\mathbf{W})p(\mathbf{Z}|\mathbf{W})p(\mathbf{H})p(\mathbf{V}|\mathbf{H},\mathbf{Z})=$$

$$\left[\prod _{f=1}^F\left(p({\mathbf{w}}_f)\prod _{n=1}^{N}p({\mathbf{z}}_{fn}|{\mathbf{w}}_f)\right)\right]\left[\prod _{n=1}^N\left(\prod _{k=1}^{K}p(h_{kn})\prod _{f=1}^{F}p(v_{fn}|{\mathbf{h}}_n,{\mathbf{z}}_{fn})\right)\right].$$

$$p(\mathbf{V},\mathbf{Z})=$$

$$\prod _f\stackrel{p({\underline{\mathbf{Z}}}_f)}{\overbrace{\int p({\mathbf{w}}_f)\prod _{n}p({\mathbf{z}}_{fn}|{\mathbf{w}}_f)d{\mathbf{w}}_f}}\prod _n\stackrel{p({\mathbf{v}}_n|{\mathbf{Z}}_n)}{\overbrace{\int \prod _{k}p(h_{kn})\prod _{f}p(v_{fn}|{\mathbf{h}}_n,{\mathbf{z}}_{fn})d{\mathbf{h}}_n}},$$

$$L_{fk}$$

$$A_{kn}$$

$$p({\underline{\mathbf{Z}}}_f)$$

$$p({\mathbf{v}}_n|{\mathbf{Z}}_n)$$

$$L_{fk}^{\neg fn}$$

$$A_{kn}^{\neg fn}$$

$$p({\mathbf{z}}_{fn}|{\mathbf{Z}}_{\neg fn},\mathbf{V})\propto \prod _k{\left[({\gamma }_k+L_{fk}^{\neg fn})\frac{({\alpha }_k+A_{kn}^{\neg fn}{)}^{v_{fn}}({\beta }_k+B_{kn}^{\neg fn}{)}^{{\bar{v}}_{fn}}}{{\alpha }_k+{\beta }_k+M_{kn}^{\neg fn}}\right]}^{z_{fkn}}.$$

$${\mathbf{w}}_f|{\underline{\mathbf{Z}}}_f$$

$$h_{kn}|{\mathbf{Z}}_n,{\mathbf{v}}_n$$

$${\mathbb{E}}_{{\mathbf{w}}_f}[{\mathbf{w}}_f|{\underline{\mathbf{Z}}}_f]$$

Full text

1 Introduction

Nonnegative matrix factorization (NMF) is a family of methods that approximate a nonnegative matrix $\mathbf{V}$ of size ${F\times N}$ as the product of two nonnegative matrices,

[TABLE]

where $\mathbf{W}$ has size $F\times K$, and $\mathbf{H}$ has size ${K\times N}$, often referred to as the dictionary and the activation matrix, respectively. $K$ is usually chosen such that $FK+KN\ll~{}FN$, hence reducing the data dimension.

Such an approximation is often sought after by minimizing a measure of fit between the observed data $\mathbf{V}$ and its factorized approximation $\mathbf{WH}$, i.e.,

[TABLE]

where $D$ denotes the cost function, and where the notation $\boldsymbol{\mathbf{A}}\geq 0$ denotes nonnegativity of the entries of $\boldsymbol{\mathbf{A}}$. Typical cost functions include the squared Euclidean distance and the generalized Kullback-Leiber divergence (Lee and Seung, 2001), the $\alpha$-divergence (Cichocki et al., 2008) or the $\beta$-divergence (Févotte and Idier, 2011). Most of these cost functions underlie a probabilistic model for the data, such that minimization of the cost function is equivalent to joint maximum likelihood estimation of the factors (Singh and Gordon, 2008), i.e.,

[TABLE]

where $p$ is a probability distribution. As such, so-called Bayesian NMF can be considered, where the factors $\mathbf{W}$ and $\mathbf{H}$ are assumed to be random variables with prior distributions, and inference is based on their posterior distribution, i.e.,

[TABLE]

This has notably been addressed for different models such as Poisson (Cemgil, 2009), additive Gaussian (Schmidt et al., 2009; Alquier and Guedj, 2017), or multiplicative Exponential (Hoffman et al., 2010).

In this paper, we are interested in Bayesian NMF for binary data matrices. Binary matrices may represent a large variety of data such as social networks, voting data, gene expression data, or binary images. As we shall see in Section 2, a common practice is to consider the model

[TABLE]

where $\phi$ is a link function that maps the factorization $\mathbf{WH}$ to the $[0,1]$ range.111Distributions used throughout the article are formally defined in Appendix A. Although link functions are convenient since they allow the factors to be unconstrained, and sometimes result in tractable problems, they sacrifice the mean-parameterization of the Bernoulli likelihood (i.e. $\mathbb{E}[\mathbf{V}|\mathbf{WH}]=\phi(\mathbf{WH})$ instead of $\mathbb{E}[\mathbf{V}|\mathbf{WH}]=\mathbf{WH}$).

Mean-parameterized nonnegative binary matrix factorization (NBMF), however, does not rely on a link function —or equivalently, considers $\phi(\boldsymbol{\mathbf{WH}})=\boldsymbol{\mathbf{WH}}$— and assumes the likelihood of the data to be

[TABLE]

which implies $\mathbb{E}[\mathbf{V}|\mathbf{WH}]=\mathbf{WH}$. Mean-parameterization is an interesting property of a model because it makes the decomposition easy to interpret. Besides, in a Bernoulli likelihood, the product $\mathbf{WH}$ —and, in this paper, the individual factors as well— can be interpreted as probabilities. An additional advantage of dealing with probabilities is that they lay in a linear, continuous space, where we can apply off-the-shelf clustering methods over the latent factors. For instance, in recommender systems, we may want to cluster users by latent musical preferences, or by the latent type of product they buy. Or we may want a user to see only those categories with a probability higher than some threshold.

Our contributions in this paper are the following:

1. (a)

we present a unified framework for three Bayesian mean-parameterized NBMF models that place three possible constraints on the factors; 2. (b)

we derive a collapsed Gibbs sampler as well as collapsed variational inference algorithms, which have never been considered for these models; 3. (c)

we discuss the extension of the models to a nonparametric setting —where the number of latent components does not need to be fixed a priori— and propose an approximation that shows excellent results with real data.

We test the performance of the models for different tasks in multiple datasets and show that our models give similar or superior results to the state of the art, while automatically detecting the number of relevant components. The datasets, the algorithms, and scripts to replicate all the reported results are available through an R package. 222https://github.com/alumbreras/NBMF

2 Related work

2.1 Logistic PCA family

One of the earliest probabilistic approaches to model binary data matrices comes from PCA-related methods. The reformulation of PCA as a probabilistic generative model with a Gaussian likelihood (Sammel et al., 1997; Tipping and Bishop, 1999) opened the door to considering other likelihoods such as Bernoulli models, which are more appropriate for binary observations. We refer to it as logistic PCA. The maximum likelihood estimator in the model is given by

[TABLE]

where $\sigma$ is the logistic function $\sigma(x)=1/(1+e^{-x})$. Note that in this model the expectation is a non-linear transformation of the factors, such that $\mathbb{E}[\boldsymbol{\mathbf{V}}|\boldsymbol{\mathbf{WH}}]=\sigma(\boldsymbol{\mathbf{WH}})$.

There are multiple maximum likelihood estimation algorithms for logistic PCA. For instance, while Sammel et al. (1997) use a Monte-Carlo Expectation Minimization (MC-EM) algorithm, Tipping (1999) derives a faster variational EM (vEM) algorithm. Collins et al. (2002) generalize probabilistic PCA to the exponential family and propose a general algorithm that exploits the duality between likelihoods in the exponential family and Bregman divergences. Later, Schein et al. (2003) improved the algorithm of Collins et al. (2002), thanks to the optimization of a tight upper bound by Alternate Least Squares (ALS). Finally, note that inference could also be tackled with Polya-Gamma data augmentation schemes Polson et al. (2013).

Other models similar to logistic PCA have been proposed with various priors or constraints over the factors. Some examples are Hernandez-Lobato et al. (2014), where the factors are given Gaussian priors, Tomé et al. (2013), which allows one factor to have negative values, and Larsen and Clemmensen (2015), where both factors are nonnegative. Meeds et al. (2007) consider the same logistic link function but a three factor decomposition $\sigma(\mathbf{WXH}$), where $\mathbf{W}$ and $\mathbf{H}$ are binary factors that represent cluster assignments, and $\mathbf{X}$ is a real-valued matrix that encodes the relations between the clusters. The expectation in these models is always $\mathbb{E}[\boldsymbol{\mathbf{V}}|\boldsymbol{\mathbf{WH}}]=\sigma(\boldsymbol{\mathbf{WH}})$ or $\mathbb{E}[\boldsymbol{\mathbf{V}}|\boldsymbol{\mathbf{WXH}}]=\sigma(\boldsymbol{\mathbf{WXH}})$.

2.2 Poisson matrix factorization

For practical reasons, some works have considered Poisson matrix factorization (PMF) techniques for binary data. In this case the binary nature of the data is ignored and a Poisson likelihood is considered:

[TABLE]

Different flavors of PMF have been proposed, in frequentists or Bayesian settings, and can be found, for example, in Lee and Seung (2001); Canny (2004); Cemgil (2009); Zhou et al. (2012); Gopalan et al. (2014, 2015). An advantage of PMF is that it is mean-parameterized, i.e., $\mathbb{E}[\boldsymbol{\mathbf{V}}|\boldsymbol{\mathbf{WH}}]=\boldsymbol{\mathbf{WH}}$. Another useful advantage is that inference algorithms need only iterate over non-zero values, which makes them very efficient for sparse matrices. In our case, a significant disadvantage is that it assigns non-zero probabilities to impossible observations ($v_{fn}>1)$.

As we discussed above, a more reasonable choice consists in replacing the Poisson distribution with a Bernoulli distribution, possibly using some link function that maps the parameter into a $[0,1]$ range, ensuring a valid Bernoulli parameter. Unfortunately, unlike Poisson models, zeroes and ones under Bernoulli likelihoods do not represent counts but classes —a zero can be considered as a no, while a one can be considered as a yes— and the algorithms need to iterate over all the elements of the observation matrix. To bypass this, Zhou (2015); Zhou et al. (2016) proposed using the alternative link function $f(x)=1-e^{-x}$, coined Bernoulli-Poisson, such that

[TABLE]

Thanks to the new link function, the model can be “augmented” to a Poisson model by introducing latent variables $c_{fn}$ such that

[TABLE]

where $\mathbbm{1}$ is the set indicator function. By placing conjugate Gamma priors over the factors, posteriors can be obtained by Gibbs sampling. Expectation in this model is $\mathbb{E}[\boldsymbol{\mathbf{V}}|\boldsymbol{\mathbf{WH}}]=f(\boldsymbol{\mathbf{WH}})$. Figure 1 shows the logistic and Bernoulli-Poisson functions. Note that each link function has a different input domain, leading to different priors or constraints over the factors.

2.3 Bernoulli mean-parameterized matrix factorization

None of the above methods offers mean-parameterization and assumes a Bernoulli distribution over the data. The family of mean-parameterized Bernoulli models is the basis of the models presented in this paper. These models have first been introduced in Kabán and Bingham (2008) (binary ICA) and Bingham et al. (2009) (Aspect Bernoulli model) where the constraints $\sum_{k}w_{fk}=1,h_{kn}\in[0,1]$, or vice versa, are imposed on the factors. The trick is that these constraints induce convex combinations of binary elements such that $\sum_{k}w_{fk}h_{kn}\in[0,1]$, which gives a valid Bernoulli parameter. In Bingham et al. (2009) the constraints are imposed explicitly, and the maximum likelihood estimator is computed by using EM with an augmented version of the model. In Kabán and Bingham (2008) the constraint is imposed through Dirichlet and Beta priors over the factors, and Variational Bayes (VB) estimation of their posteriors is derived exploiting a similar augmentation scheme.

Table 1 presents a summary of the methods presented in Sections 2.1-2.3 that use a Bernoulli likelihood.

2.4 Others

Some models have also been proposed to find binary decompositions, that is, matrix factorizations where $\boldsymbol{\mathbf{W}}$ and $\boldsymbol{\mathbf{H}}$ contain binary elements. For instance, Zhang et al. (2009) aim to minimize a Euclidean distance or, equivalently, maximize a Gaussian likelihood under the binary constraint. Similarly, the “discrete basis problem” of Miettinen et al. (2008) aims to minimize a $L1$-norm under the same constraint. In Slawski et al. (2013), an algorithm is proposed to retrieve the exact factorization when one of the factors is constrained to be binary, and the other one to be stochastic, i.e., $\sum_{k}h_{kn}=1$. More recently Rukat et al. (2017) proposed a Bayesian model for the Boolean matrix factorization problem, where $\boldsymbol{\mathbf{W}}\boldsymbol{\mathbf{H}}$ is a Boolean product. Çapan et al. (2018) proposed sum-conditioned Poisson factorization models that apply to binary data.

In the general model defined by Eq. (5), we are essentially assuming that ${\mathbf{V}}\approx\phi(\boldsymbol{\mathbf{W}}\boldsymbol{\mathbf{H}})$. This can be seen as a one-layer generative network with input $\boldsymbol{\mathbf{H}}$, weight $\boldsymbol{\mathbf{W}}$ and non-linearity $\phi(\cdot)$. As such it is possible to conceive more general models by stacking various layers or modeling $\boldsymbol{\mathbf{W}}$ and $\boldsymbol{\mathbf{H}}$ as the outputs of deep networks themselves. In the context of recommendation systems, where binary matrices can represent either binary ratings or implicit feedback, this has been for example considered in He et al. (2017); Xue et al. (2017).

3 Mean-parameterized Bernoulli models

3.1 Models

Let us consider a mean-parameterized Bernoulli model for an observed binary matrix $\mathbf{V}$ and two latent factors $\mathbf{W}$, $\mathbf{H}$:

[TABLE]

To guarantee valid Bernoulli parameters, we can impose three possible sets of constraints on the factors such that $\sum_{k}w_{fk}h_{kn}\in[0,1]$:

[TABLE]

In a Bayesian setting, we may place Beta and Dirichlet priors over the factors to respect these constraints:

[TABLE]

where $\mathbf{h}_{n}$ denotes the $n$-th column of the matrix $\mathbf{H}$, and $\mathbf{w}_{f}$ denotes the $f$-th row of the matrix $\mathbf{W}$. The Beta parameters are positive real numbers $\alpha_{k},\beta_{k}\in\mathbb{R}_{++}$ and the Dirichlet parameters are $K$-dimensional vectors of positive real numbers $\boldsymbol{\mathbf{\gamma}},\boldsymbol{\mathbf{\eta}}\in\mathbb{R}_{++}^{K}$.

Note that each element $w_{fk}$ and $h_{kn}$ can be interpreted as a probability. We can either merely impose that the elements of a row $\mathbf{w}_{f}$ or a column $\mathbf{h}_{n}$ lie between [math] and $1$, or, more strongly, that they sum up to one. This implies a difference in modeling. On the one hand, imposing that the elements lie between [math] and $1$ induce non-exclusive components. On the other hand, the sum-to-one constraint induces exclusive components, i.e., the more likely a component is, the less likely are the others. Fig. 2 displays simulated matrices generated from each of the models.

The first two models, Beta-Dir and Dir-Beta, are symmetric. Indeed, estimating $\mathbf{W}$ (resp., $\mathbf{H}$) in one model is equivalent to estimating $\mathbf{H}$ (resp., $\mathbf{W}$) in the other model after transposing the matrix $\mathbf{V}$. As such, in the rest of the paper, we will only consider the Beta-Dir model and the Dir-Dir model.

The Aspect Bernoulli model of Bingham et al. (2009) is built over Eq. (11), and considers that the factors $\mathbf{W}$ and $\mathbf{H}$ are deterministic parameters which satisfy the constraint (c2). The factors are estimated by maximum likelihood with EM. The binary ICA of Kabán and Bingham (2008) corresponds to the Beta-Dir model, and inference is performed with VB. In the following sections, we present inference methods for the posterior distributions of $\mathbf{W}$ and $\mathbf{H}$ in the Beta-Dir and Dir-Dir models. Because these distributions are intractable, we propose novel collapsed Gibbs sampling and collapsed variational inference strategies. We also derive a nonparametric approximation where the number of latent dimensions $K$ does not need to be fixed a priori.

3.2 Connections with latent Dirichlet allocation (LDA).

The proposed models have connection with topic models, in particular with LDA (Blei et al., 2003), described by

[TABLE]

Here, $\mathbf{V}$ is the so-called “bag-of-words” representation of a corpus of documents, i.e., $v_{fn}$ represents the number of occurrences of word $f$ in document $n$. $L_{n}$ is the total number of words in document $n$, i.e., $L_{n}=\sum_{f}v_{fn}$. The columns of $\mathbf{W}$ are assumed to sum to 1, such that $w_{fk}$ represents the probability of word $f$ in topic $k$, and $h_{kn}$ represents the probability of topic $k$ in document $n$. In standard LDA, no prior distribution is assumed on $\mathbf{W}$, though it is common practice to assume Dirichlet distributions column-wise. The observation model is multinomial, with a total “budget” of $L_{n}$ words to distribute in each document $n$. In contrast, the binary data models presented in this article are based on independent Bernoulli observations, and cannot be recast as multinomial. An LDA-like model may arise if we assume that each column of $\mathbf{V}$ contains only one 1, and the rest are zeros (i.e. $L_{n}=1$ for all $n$). This setting is however not considered in the article.

3.3 Bayesian inference

In this paper, we opt for Bayesian inference of the latent factors ${\mathbf{W}}$ and ${\mathbf{H}}$. This means that we seek to characterize the posterior distribution $p({\mathbf{W}},{\mathbf{H}}|{\mathbf{V}})$. The posterior in the considered models is not available in closed form and we will resort to numerical approximations (Markov Chain Monte Carlo sampling, variational inference). Another possible route is to seek point estimates through either constrained maximum likelihood (ML) or maximum a posteriori (MAP) estimation. For example, given Eq. (11) and the set of constraints (c1), maximum likelihood estimation writes

[TABLE]

Such optimization problems may be tackled with, e.g., majorization-minimization (Bingham et al., 2009) or proximal gradient descent (Udell et al., 2016). Though it can be computationally more costly, we favored Bayesian inference for the following reasons. Characterizing the full posterior provides a measure of uncertainty over the latent parameters and over predicted values of ${\mathbf{V}}$. As we will show later, it allows in turn to infer the rank $K$ of the factorization, a very desirable property that is more difficult to obtain with point estimation methods. Finally, Bayesian inference is rather customary in topic models such as LDA or Discrete Component Analysis (Buntine and Jakulin, 2006), and our work intends to follow similar principles.

4 Inference in the Beta-Dir model

In this section, we derive a collapsed Gibbs sampler (Liu, 1994) for the Beta-Dir model. First, we will augment the model with latent indicator variables $\boldsymbol{\mathbf{Z}}$ so that it becomes conjugate. Then the collapsed Gibbs sampler consists in marginalizing out the factors $\boldsymbol{\mathbf{W}}$ and $\boldsymbol{\mathbf{H}}$, thus running a Gibbs sampler over the indicator variables $\boldsymbol{\mathbf{Z}}$ only. The interest of collapsed Gibbs sampling is that it offers improved mixing properties of the Markov chains, i.e., better exploration of the parameter space, thanks to the reduced dimensionality. We use a superscript, as in $x^{(j)}$, to indicate the $j$-th sample of a chain (after burn-in). After sampling, given a collection of samples $\boldsymbol{\mathbf{Z}}^{(j)}$ from the posterior, we will be able to directly sample from the posteriors of interest $p(\boldsymbol{\mathbf{W}}|\boldsymbol{\mathbf{V}})$ and $p(\boldsymbol{\mathbf{H}}|\boldsymbol{\mathbf{V}})$.

4.1 Augmented model

We can augment the Beta-Dir model with indicator variables $\boldsymbol{\mathbf{z}}_{fn}$, that contain component assignments from the Dirichlet factor, as shown in Bingham et al. (2009). More precisely, $\boldsymbol{\mathbf{z}}_{fn}$ is a vector of dimension $K$ with elements $z_{fkn}\in\{0,1\}$ such that only one element equals to one and all the others equal to zero. In other words, $\boldsymbol{\mathbf{z}}_{fn}\in\{\boldsymbol{\mathbf{e}}_{1},\ldots,\boldsymbol{\mathbf{e}}_{K}\}$, where $\boldsymbol{\mathbf{e}}_{k}$ is the $k$-th canonical vector of $\mathbb{R}^{K}$. The augmented model is a mixture model described by:333Some readers may be more accustomed to the alternative notation where the “one-hot” variable $\boldsymbol{\mathbf{z}}_{fn}$ is replaced by an integer-valued index $z_{fn}\in\{1,\ldots,K\}$. In this case, the Bernoulli parameter in Eq. (18) becomes $h_{z_{fn}n}$.

[TABLE]

Indeed, marginalizing $\boldsymbol{\mathbf{z}}_{fn}$ from Eqs. (17)-(18) leads to Eq. (11) as shown next. From Eqs. (17)-(18) we have

[TABLE]

and thus $p(v_{fn}=1|{\mathbf{w}}_{f},\boldsymbol{\mathbf{h}}_{n})=[\boldsymbol{\mathbf{W}}\boldsymbol{\mathbf{H}}]_{fn}$ and $p(v_{fn}=0|{\mathbf{w}}_{f},\boldsymbol{\mathbf{h}}_{n})=1-[\boldsymbol{\mathbf{W}}\boldsymbol{\mathbf{H}}]_{fn}$, i.e., $v_{fn}$ has the marginal distribution given by Eq. (11). A graphical representation of the augmented model is shown in Fig. 3-(a). Let us think of a recommender system application, where columns of $\boldsymbol{\mathbf{V}}$ are users, and rows are items. An interpretation of the above model is the following. Each item $f$ is characterized by a probability over topics, $\boldsymbol{\mathbf{w}}_{f}$. Then, for each user-item pair, a topic $k$ (indicated by $\boldsymbol{\mathbf{z}}_{fn}$) is activated with probability $w_{fk}$, and the probability that the user $n$ consumes this item ($v_{fn}=1$) is $h_{kn}$.

Denoting by $\mathbf{Z}$ the $F\times K\times N$ tensor with entries $z_{fkn}$, the joint probability in the augmented model is given by:

[TABLE]

4.2 Collapsed Gibbs sampler

Thanks to the previous augmentation, and exploiting conjugacy, we can now marginalize out $\boldsymbol{\mathbf{W}}$ and $\boldsymbol{\mathbf{H}}$ from Eq. (21). The marginalized distribution has the following structure:

[TABLE]

where $\underline{\mathbf{Z}}_{f}$ denotes the $K\times N$ matrix with entries $\{z_{fkn}\}_{kn}$, $\mathbf{Z}_{n}$ denotes the $F\times K$ matrix with entries $\{z_{fkn}\}_{fk}$. Let us define the following four variables that act as counters:

[TABLE]

where $\bar{v}_{fn}=1-v_{fn}$. $A_{kn}$ and $B_{kn}$ count how many times the component $k$ is associated to a “positive” observation ($v_{fn}=1$) and to a “negative” observation ($\bar{v}_{fn}=1$) in the $n$-th sample. Then we have the following expressions (derivations are available in Appendix B.1):

[TABLE]

The posterior of the indicator variables $p(\boldsymbol{\mathbf{Z}}|\boldsymbol{\mathbf{V}})$ is not available in closed form and the proposed collapsed Gibbs sampler consists in iteratively sampling each vector $\boldsymbol{\mathbf{z}}_{fn}$ given the current value of the other indicator vectors. Let $L_{fk}^{\neg fn}$, $M_{kn}^{\neg fn}$, $A_{kn}^{\neg fn}$, $B_{kn}^{\neg fn}$ be the state of the counters when the tube $(f,n)$ of the tensor $\boldsymbol{\mathbf{Z}}$ is left out of the sums:

[TABLE]

In Appendix B.2 we show that the conditional posterior of $\boldsymbol{\mathbf{z}}_{fn}$ given the remaining variables $\boldsymbol{\mathbf{Z}}_{\neg fn}$ is given by:

[TABLE]

This expression needs to be normalized to ensure a valid probability distribution. This can be easily done by computing the right-hand side of Eq. (25) for every of the $K$ possible values of $\boldsymbol{\mathbf{z}}_{fn}$ and normalizing by the sum. Eq. (25) shows that the probability of choosing a component $k$ depends on the number of elements already assigned to that component. More precisely, it depends on the one hand on the number of elements assigned to component $k$ in column $n$. On the other hand, it also depends on the proportion of elements in row $f$ assigned to component $k$ that explain ones (if $v_{fn}=1$) or zeros (if $v_{fn}=0$) in $\boldsymbol{\mathbf{V}}$ in the total number of elements associated to $k$ in that row (see Figure 4). The parameters $\gamma_{k},\alpha_{k},\beta_{k}$ act as pseudo-counts: they give a priori belief about how many elements are assigned to each component.

Our collapsed Gibbs sampling is summarized in Alg. (1). Note that, although Alg. (1) does not explicitly include it, we must draw samples during an initial burn-in phase (as required by any MCMC method) before collecting the last $J$ samples, after the chain has converged to the stationary distribution. Note also that the algorithm can readily deal with incomplete matrices by simply skipping missing entries (i.e., the loop over $f$ and $n$ only runs over available entries).

Latent factors posteriors. Thanks to conjugacy, the conditional posteriors of $\mathbf{W}$ and $\mathbf{H}$ given $\boldsymbol{\mathbf{Z}}$ and $\boldsymbol{\mathbf{V}}$ are given by:

[TABLE]

The conditional posterior expectations are given by:

[TABLE]

where we used the equalities $\sum_{kn}z_{fkn}=N$ and $A_{kn}+B_{kn}=\sum_{f}z_{fkn}$. Using the law of total expectation, i.e., $\mathbb{E}[X]~{}=~{}\mathbb{E}_{Y}[\mathbb{E}_{X}[X|Y]]$, and given a set a samples $\boldsymbol{\mathbf{Z}}^{(j)}$, it follows that the marginal posterior expectations of the latent factors can be computed as:

[TABLE]

Prediction. The predictive posterior distribution of an unseen data sample $v_{fn}^{*}$ given the available data $\boldsymbol{\mathbf{V}}$ is given by

[TABLE]

Because the predictive posterior is a Bernoulli distribution, its expectation is given by $\mathbb{E}[\boldsymbol{\mathbf{w}}_{f}\boldsymbol{\mathbf{h}}_{n}|\boldsymbol{\mathbf{V}}]$, which can be approximated using samples $\boldsymbol{\mathbf{W}}^{(j)}$, $\boldsymbol{\mathbf{H}}^{(j)}$ from the distributions given by Eqs. (26)-(27) given $\boldsymbol{\mathbf{Z}}=\boldsymbol{\mathbf{Z}}^{(j)}$:

[TABLE]

4.3 Collapsed variational inference

Given the collapsed model of Eqs. (22)-(24) we may derive a mean-field Collapsed Variational Bayes algorithm (CVB) (Teh et al., 2007) by assuming that the posterior factorizes as $q(\boldsymbol{\mathbf{Z}})=\prod_{fn}q(\boldsymbol{\mathbf{z}}_{fn})$. The key of CVB is that its free energy is a strictly better bound on the evidence than the free energy of the standard, i.e., uncollapsed, VB. We compute the CVB updates by applying the mean-field VB updates to the collapsed model:

[TABLE]

where the expectations are taken over the variational posterior. This leads us to

[TABLE]

The expectations of the form $\mathbb{E}_{q(z)}[\log(x+z)]$ are expensive to compute. A simpler alternative is CVB0 (Asuncion et al., 2009), which uses a zero-order Taylor approximation $\mathbb{E}_{q}(z)[\log(x+z)]\approx\log(x+\mathbb{E}_{q(z)}[z])]$ and has been shown to give, in some cases, better inference results than CVB. Under the CVB0 version our update becomes

[TABLE]

which has a similar structure to the collapsed Gibbs sampler in Eq. (25). Overall, the collapsed VB algorithm has the same structure as the Gibbs sampler summarized in Alg. (1). Note that when the data matrix is too large for batch processing, one can routinely resort to stochastic variational inference (Hoffman et al., 2013).

Latent factors posteriors. The variational distributions of the factors can be obtained from the uncollapsed version:

[TABLE]

The Taylor approximation breaks the theoretical guarantees of their superiority over the ones given by uncollapsed VB. Still, they have been reported to work better than VB in practice. Another drawback of the approximation is the loss of convergence guarantees. Although we do not address this issue here, this has been recently addressed by Ishiguro et al. (2017), where an annealing strategy is used to gradually decrease the portion of the variational posterior changes.

Prediction. The predictive posterior can be computed as in Eq. (32), and its expectation is computed using the variational approximations of the factors, i.e.,

[TABLE]

4.4 Approximating infinite components

Recall that in the augmented model, the component assignments $\boldsymbol{\mathbf{z}}_{fn}$ have a Discrete distribution such that (Eqs. (16)-(17))

[TABLE]

The variable $\boldsymbol{\mathbf{w}}_{f}$ may be integrated out leading to the expression of $p(\underline{\boldsymbol{\mathbf{Z}}}_{f})$ given by Eq.(23). In Appendix B.2, we show that the prior conditionals are given by:

[TABLE]

Let us assume from now that the Dirichlet prior parameters are such that $\gamma_{k}=\gamma/K$, where $\gamma$ is a fixed nonnegative scalar, so that:

[TABLE]

The conditional prior given by Eq. (41) is reminiscent of the Chinese Restaurant process (CRP) (Aldous, 1985; Anderson, 1991; Pitman, 2002). In the limit when $K\rightarrow\infty$, the probability of assigning $\boldsymbol{\mathbf{z}}_{fn}$ to component $k$ is proportional to the number $L_{fk}^{\neg fn}$ of current assignments to that component. Let $K^{+}$ denote the current number of non-empty components (i.e., such that $L_{fk}^{\neg fn}>0$). Then the probability of choosing an empty component is

[TABLE]

Note that the latter probability does not depend on $K$. In practice, we set $K$ to a large value and observed self-pruning of the number of components, hence achieving to automatic order selection, similar to Hoffman et al. (2010). Implementing exact inference in the truly nonparametric model

[TABLE]

is more challenging. This is because there is a CRP for each feature $f$, and some empty components may become unidentifiable in the limit. This is a known issue that could be addressed using for example a Chinese Restaurant Franchise process (Teh et al., 2006) but is beyond the scope of this article.

5 Inference in the Dirichlet-Dirichlet model

The methodology to obtain a collapsed Gibbs sampler for the Dir-Dir model is similar to the approach followed for the Beta-Dir model.

It is possible to augment the model with the same auxiliary variable ${\mathbf{z}}_{fn}$ and compute the expression of $p(\boldsymbol{\mathbf{z}}_{fn}|\boldsymbol{\mathbf{Z}}_{\neg fn},\boldsymbol{\mathbf{V}})$ in closed form. However, the expression of the conditional posterior, given in Appendix C, involves combinatorial computations and is infeasible in practice. As such, we propose an alternative Gibbs sampler that relies on a double augmentation, presented next. Obtaining a variational collapsed algorithm for the Dir-Dir model is not straightforward, even using the double augmentation, and is left for future work.

5.1 Fully augmented model

Unlike the Beta-Dir model, the Dir-Dir model is not fully conjugate after a first augmentation. We propose a second augmentation with a new indicator variable $\boldsymbol{\mathbf{c}}_{fn}\in\{\boldsymbol{\mathbf{e}}_{1},\ldots,\boldsymbol{\mathbf{e}}_{K}\}$, that plays a similar role to $\boldsymbol{\mathbf{z}}_{fn}$ The fully augmented version is:

[TABLE]

To show that this is a valid augmentation, note that $v_{fn}$ can only be nonzero (and equal to 1) if $\boldsymbol{\mathbf{c}}_{fn}=\boldsymbol{\mathbf{z}}_{fn}$. Then, the marginal probability of $v_{fn}=1$ is given by

[TABLE]

and we thus recover the Bernoulli model of Eq. (11) as announced. Compared to the Beta-Dir model and using our recommender system analogy, this means that, in each user-item pair, the user also activates one topic, and then consumes the item if the user active topic is equal to the item active topic. The Dir-Dir model makes a stronger assumption than the Beta-Dir since the user can only activate one topic per item. A graphical representation of the fully augmented model is given in Fig. 3-(b). In the following, we denote by $\boldsymbol{\mathbf{C}}$ the $F\times K\times N$ tensor with entries $c_{fkn}$, and by $\boldsymbol{\mathbf{C}}_{n}$ the $F\times K$ matrix with entries $\{c_{fkn}\}_{fk}$.

5.2 Collapsed Gibbs sampling

In this section we show that $\boldsymbol{\mathbf{W}}$ and $\boldsymbol{\mathbf{H}}$ can be marginalized from the joint probability of the fully augmented model and then propose a collapsed Gibbs sampler for $p(\boldsymbol{\mathbf{Z}},\boldsymbol{\mathbf{C}}|\boldsymbol{\mathbf{V}})$. The joint probability is given by:

[TABLE]

where

[TABLE]

and where $\delta$ denotes the Dirac delta function. Following Section 4.2 and Appendix B.2, the prior conditional are given by:

[TABLE]

where $Q_{kn}^{\neg fn}=\sum_{f^{\prime}\not=f}c_{f^{\prime}kn}$ and $L_{fk}^{\neg fn}=\sum_{n^{\prime}\not=n}z_{fkn^{\prime}}$ is as before. When $v_{fn}=1$, $\boldsymbol{\mathbf{c}}_{fn}$ and $\boldsymbol{\mathbf{z}}_{fn}$ must be assigned to the same component ($\boldsymbol{\mathbf{c}}_{fn}=\boldsymbol{\mathbf{z}}_{fn}$). To respect this constraint, we may sample them together from the posterior. Introducing the vector $\boldsymbol{\mathbf{x}}_{fn}$ such that $\boldsymbol{\mathbf{x}}_{fn}=\boldsymbol{\mathbf{z}}_{fn}=\boldsymbol{\mathbf{c}}_{fn}$, the conditional posterior is given by:

[TABLE]

When $v_{fn}=0$, we can assign to one of the two auxiliary variables any component not currently assigned to the other auxiliary variable. The respective conditional posteriors are given by:

[TABLE]

A pseudo-code of the resulting Gibbs sampler is given in Alg. (2). As with the Beta-Dir model, we set $\gamma_{k}=\gamma/K$ with $K$ large to emulate a nonparametric setting (note that $\boldsymbol{\mathbf{\eta}}$ does not need to depend on $K$ itself).

Latent factors posteriors and prediction. The conditional posteriors of the latent factors given $\boldsymbol{\mathbf{Z}}$ and $\boldsymbol{\mathbf{C}}$ are given by:

[TABLE]

As done with the Beta-Dir model, we may use the law of total expectation and the samples $\boldsymbol{\mathbf{z}}_{fn}^{(j)}$ to obtain Monte-Carlo estimates of the posterior expectations:

[TABLE]

As in the Beta-Dir model, $\boldsymbol{\mathbf{W}}$ and $\boldsymbol{\mathbf{H}}$ can be sampled in a second step given a collection of samples of $\boldsymbol{\mathbf{Z}}$ and $\boldsymbol{\mathbf{C}}$. The predictive posterior and its expectation can be computed as in Eqs. (32), (33).

6 Experiments

We show the performance of the proposed NBMF methods for different tasks in multiple datasets. The datasets, the algorithms, and scripts to replicate all the reported results are available through our R package.

6.1 Datasets

We consider five different public datasets, described next and displayed in Fig. 5.

Animals (animals). The animals dataset (Kemp et al., 2006) contains 50 animals and 85 binary attributes such as nocturnal, hibernates, small or fast. The matrix takes $v_{fn}=1$ if animal $n$ has attribute $f$.

Last.fm (lastfm). We use a binarized subset of the Last.fm dataset (Celma, 2010) where rows correspond to users and columns correspond to musical artists. The matrix takes $v_{fn}=1$ if user $n$ has listened to artist $f$ at least once. The matrix has $F=285$ rows and $N=1226$ columns.

Paleontological data (paleo). The NOW (New and Old Worlds) fossil mammal database contains information of fossils found in specific paleontological sites (NOW, 2018). From the original paleontological data, we build a matrix where each row is a genus, each column is a location, and $v_{fn}=1$ if genus $f$ has been found at location $n$. We used the same pre-processing as in Bingham et al. (2009) (i.e., we discarded small and infrequent genus, locations with only one genus and kept locations with longitude between 0 and 60 degrees East) and obtained a matrix with $F=253$ rows and $N=902$ columns.

Catalan parliament (parliament). We created a list of the current members of the Catalan parliament and collected the information of who follows whom on Twitter (March 2018). With this data, we created a square adjacency matrix where $v_{fn}=1$ if member $f$ follows member $n$. There are seven political groups represented. The resulting matrix has 135 rows and columns ($F=N$).

UN votes (unvotes). The United Nations General Assembly Voting Data is a dataset that contains the roll-call votes in the UN General Assembly between 1946-2017 (Voeten, 2013). Votes can be yes, no, or abstention. From this data, we created a matrix where $v_{fn}=1$ if country $f$ voted yes to the call $n$, $v_{fn}=0$ if it voted no, and a missing value if the country did not participate in that call or was not a member of the UN at that time. Next, abstention votes will be treated as either negative votes (no) or missing data, as specified in each experiment. The resulting matrix has $F=200$ rows and $N=5429$ columns.

6.2 Methods and setting

State-of-the-art methods. We compare our proposed methods with the following state-of-the-art methods for binary data.

logPCA-K. Probabilistic PCA with Bernoulli likelihood. We use the algorithm presented in Collins et al. (2002). The notation logPCA-K will embed the chosen number of components $K$ (e.g., logPCA-8 signifies $K=8$). We used the R package logisticPCA (Landgraf and Lee, 2015) with default parameters.

bICA-K. The binary ICA method introduced in Kabán and Bingham (2008), which uses uncollapsed mean-field variational inference over the partially augmented model (Eqs. (15)-(18)). This is also a parametric method that requires setting $K$.

Proposed methods. Our proposed methods are as follows.

Beta-Dir GS. Estimation in the Beta-Dir model with collapsed Gibbs sampling. Beta parameters are set to $\alpha_{k}=\beta_{k}=1$. To emulate a nonparametric setting, the Dirichlet parameters are set to $\gamma_{k}=1/K$ and the number of components is set to $K=100$.

Beta-Dir VB. Estimation in the Beta-Dir model with collapsed variational Bayes (CVB0). Beta parameters are set to $\alpha_{k}=\beta_{k}=1$. To emulate a nonparametric setting, we set $\gamma_{k}=1/K$ and $K=100$.

Dir-Dir GS. Estimation in the Dir-Dir model with collapsed Gibbs sampling. To emulate a nonparametric setting, we set $\gamma_{k}=1/K$, $\eta_{k}=1$, and $K=100$.

c-bICA-K. Collapsed bICA. The algorithm corresponds to Beta-Dir VB with $\alpha_{k}=\beta_{k}=\gamma_{k}=1$. It is the collapsed version of bICA-K using CVB0 and without the nonparametric approximation.

Implementation details. For each dataset, we ran some preliminary experiments to assess the number of iterations needed by the algorithms to converge. For the Gibbs samplers, we set a conservative burn-in phase of $4{,}000$ iterations and kept the last $1{,}000$ samples of $\boldsymbol{\mathbf{Z}}$ after burn-in. A total number of $500$ iterations where used for the variational algorithms. In every experiment, we initialized the Gibbs samplers with a random tensor $\boldsymbol{\mathbf{Z}}$ such that $\boldsymbol{\mathbf{z}}_{fn}={\mathbf{e}}_{k}$ with random $k$. We did not find a special sensitivity to the initial state of the Gibbs Sampler, but we chose a random initialization as a good practice. We also tried initializing with some variational steps, which is another common practice, but did not see significant improvements. Similarly, we initialized the variational algorithms with a random $\mathbb{E}[\boldsymbol{\mathbf{Z}}]$ such that $\mathbb{E}[\boldsymbol{\mathbf{z}}_{fn}]={\mathbf{e}}_{k}$ with random $k$. Our variational algorithms are sensitive to the initial state when some components are empty. Because an empty component has a lower probability of being chosen, in practice the variational algorithms are not capable of refilling it again. Initializing from a random state with lots of used components is, therefore, a safer way to avoid these local maxima.

Estimators. The algorithms Beta-Dir GS and Dir-Dir GS return samples from the posterior of $p(\boldsymbol{\mathbf{W}},\boldsymbol{\mathbf{H}}|{\mathbf{V}})$. Point estimates of the dictionary $\boldsymbol{\mathbf{W}}$ and data expectation $\hat{{\mathbf{V}}}=\boldsymbol{\mathbf{W}}\boldsymbol{\mathbf{H}}$ are computed by averaging (posterior mean) and by Eq. (33), respectively. Beta-Dir VB, bICA and c-bICA return variational approximations of the posterior of $\boldsymbol{\mathbf{W}}$ and $\boldsymbol{\mathbf{H}}$. Point estimates of $\boldsymbol{\mathbf{W}}$ and $\hat{{\mathbf{V}}}$ are computed from the variational distribution mean and by Eq. (39). logPCA returns maximum likelihood (ML) estimates $\hat{\boldsymbol{\mathbf{W}}}$ and $\hat{\boldsymbol{\mathbf{H}}}$. The data expectation is computed as $\hat{{\mathbf{V}}}=\sigma(\hat{{\mathbf{W}}}\hat{{\mathbf{H}}})$.

Sensitivity to hyperparameters. When using Dirichlet priors, inferences may be quite sensitive especially for small values of its concentration parameter (Steck and Jaakkola, 2003). For the setting described in Section 6.4, we have tested our algorithms under toy data generated from the model, setting the concentration parameter to 1. We have repeated different inferences with the concentration parameter of the estimator ranging from 0.1 to 10 and observed very small variations on the perplexity (around 0.1, and decreasing as the size of the observed data increases). Thus, for the sake of simplicity, we have therefore decided to set the concentrations parameters to 1 which corresponds to a uniform Dirichlet prior.

Computational cost. The time complexity of the collapsed algorithms is $\mathcal{O}(FKN)$ (assuming, for the sake of simplicity, that the Multinomial random number generator is $\mathcal{O}(1)$). Note that, unlike some non-mean-parameterized Bernoulli models (Zhou, 2015), or all Poisson models, zeros cannot be ignored because they represent another category, not a lack of observation or a zero-count.

6.3 Dictionary learning and data approximation

6.3.1 Experiments with the parliament dataset

First, we want to form an idea of how well the different models can fit original data. We focus on the parliament dataset, which has reasonable size and a clear structure. We applied the three proposed nonparametric methods Beta-Dir GS, Beta-Dir VB, Dir-Dir GS and the state-of-the-art methods bICA and logPCA. For each method, we compute the negative log-likelihood of the data approximation $\hat{{\mathbf{V}}}$, which serves as a measure of fit:

[TABLE]

bICA was run with increasing values of $K$ and the fit ceased increasing for $K=8$ which is the value used in the results (note that in this case $\hat{{\mathbf{V}}}$ is the posterior mean estimate and not the ML estimate, so the likelihood is not meant to increase monotonically). logPCA was run with the same value $K=8$. The data approximations $\hat{{\mathbf{V}}}$ and dictionaries obtained with the different methods are displayed in Figs. 6 and 7, respectively.

In terms of data approximation, Beta-Dir VB achieves the best fit among the mean-parameterized models in terms of negative log-likelihood ($4{,}729$) followed by Beta-Dir GS ($4{,}863$). The dictionaries returned by these two algorithms are very similar, with only nine active components. bICA-8 comes next in terms of fit ($4{,}957$). We also applied bICA with $K>8$ components but this did not substantially improve the likelihood. Dir-Dir GS returns the worst fit ($8{,}930$), with only two active components. Overall logPCA-8 returns the smallest negative log-likelihood ($1{,}783$). This due to its larger flexibility as compared to the mean-parameterized models (real-valued factors $\boldsymbol{\mathbf{W}}$ and $\boldsymbol{\mathbf{H}}$, with product $\boldsymbol{\mathbf{W}}\boldsymbol{\mathbf{H}}$ mapped to $[0,1]$). However, this is at the cost of meaningfulness of the decomposition, as shown in Fig. 7 and explained next.

In dictionary learning, we want to learn a meaningful decomposition of the data. The columns of the dictionary $\boldsymbol{\mathbf{W}}$ are expected to contain patterns or prototypes characteristic of the data. In particular, NMF is known to produce so-called part-based representations (each sample, a column of $\boldsymbol{\mathbf{V}}$, is approximated as a constructive linear combination of building units) (Lee and Seung, 1999). When the rows of the dictionary are given Dirichlet priors, ${\mathbf{w}}_{f}$ can also be interpreted as the probability distribution of feature $f$ over the $K$ components. In Fig. 7, the rows of the dictionaries displayed correspond to members of the parliament (MP). For each MP we show its Twitter username and its political party. The dictionaries returned by the mean-parameterized factorization methods are easily interpretable. In particular, the dictionaries returned by Beta-Dir GS, Beta-Dir VB and to some extent bICA-8 closely reflect the party memberships of the MPs. Dir-Dir GS, which is based on a less flexible model, only captures two sets of MPs, one with the members of Cs (the main opposition party) and the other with members of the remaining parties, regardless of political alignment (left-wing, right-wing, independentist and anti-independentist). In contrast, the dictionary returned by logPCA-8 is much more difficult to interpret.

6.3.2 Experiments with the unvotes dataset

In this section, we consider a subset of unvotes, reduced to the 1946-1990 range which corresponds to the Cold War period. Furthermore, the abstentions are here treated as missing values. Fig. 8 shows the dictionaries learned by the five considered methods. As before, bICA was applied with various values of $K$ and we selected the value that leads to the smallest negative log-likelihood ($K=7$). Accordingly, logPCA was also applied with $K=7$.

Fig. 8 shows that Beta-Dir GS returns the finest dictionary, detecting political blocks that tended to vote similarly in the UN assembly and capturing some nuances that the other algorithms do not find. European countries (and members or allies of NATO such as the USA, Japan, or Australia) are concentrated in one component, denoting similar voting strategies. The former members of the Soviet Union and the Warsaw Pact also form a block of their own, with some allies such as Cuba or the former Yugoslavia. Members of the Non-Aligned Movement (even countries that became members after 1991, such as Guatemala, Thailand, or Haiti), from Egypt to Cuba and from Honduras to Haiti, are split into two blocks: the Latin American group and the Asian-African group. Another detected alliance is between the United States and Israel, which are distributed between the European component and a component of their own. Beta-Dir VB detects the split between the Warsaw and NATO blocks, and the alliance between the USA and Israel, but it fails to detect the two subgroups of the Non-Aligned Movement, which is considered a single block. bICA-7 returns similar results to Beta-Dir GS but fails to detect the alliance between the USA and Israel. Note that the results of bICA are obtained with a well-chosen value of $K$ while Beta-Dir GS automatically detects a suitable value. The underlying assumption of Dir-Dir (one topic per country and one topic per vote) seems too simplistic for this dataset, and the algorithm puts every country in the same component. Again and as somewhat expected, the dictionary learned by logPCA is more difficult to interpret.

6.3.3 Experiments with the paleo dataset

We finally look into the dictionaries returned by the five considered method on the paleo dataset, see Fig. 9. The same strategy was applied to find a suitable value of $K$ for bICA and logPCA, leading to $K=7$. The results can be read as the probability of a genus to be found in a set of prototypical locations. Interestingly, Dir-Dir GS is the method that returns the most detailed dictionary for this dataset. The other methods tend to produce larger clusters of genera. This highlights the importance of choosing the right model for each dataset since they imply different underlying assumptions. Dir-Dir GS assumes one topic per genus and one topic per location. We will see in the following section that Dir-Dir GS also gives the best predictions for this dataset. Again, the dictionary obtained with logPCA is harder to interpret.

6.4 Prediction

6.4.1 Experimental setting

We now evaluate the capability of the five previously considered methods together with c-bICA to predict missing data. For each of the five considered datasets, we applied the algorithms to a $75\%$ random subset of the original data. We here use the full unvotes in which abstentions are treated as negative votes ($v_{fn}=0$). bICA and c-bICA where applied with $K=2,\ldots,8$. logPCA was applied with $K=2,\ldots,4$. Then we computed the perplexity of the test set (the 25% held-out entries) given the estimate $\hat{\boldsymbol{\mathbf{V}}}=\mathbf{E}[\boldsymbol{\mathbf{WH}}|\boldsymbol{\mathbf{V}}_{\text{train}}]$ (for all methods except logPCA-K) or $\hat{\boldsymbol{\mathbf{V}}}=\sigma(\hat{\boldsymbol{\mathbf{W}}}\hat{\boldsymbol{\mathbf{H}}})$ (for logPCA). The perplexity is here simply taken as the negative log-likelihood of the test set (Hofmann, 1999):

[TABLE]

where $T$ is the number of elements in the test set (in our case, $T=0.25FN$).

6.4.2 Prediction performance

Fig. 10 displays the perplexities obtained by all methods from 10 repetitions of the experiment with randomly selected training and test sets, and random initializations (the same starting point is used for bICA and c-bICA). The proposed Beta-Dir VB performs similarly or better (lastfm, parliament) than bICA, while automatically adjusting the number of relevant components. As hinted from the dictionary learning experiments, Dir-Dir GS performs considerably better than the other mean-parameterized methods on the paleo dataset. c-bICA does not specifically improve over bICA (remember they are based on the same model, only inference changes) and performs worse in some cases (lastfm, parliament). However, its performance is more stable, with less variation between different runs, a likely consequence of the collapsed inference.

Despite its flexibility (unconstrained $\boldsymbol{\mathbf{W}}$ and $\boldsymbol{\mathbf{H}}$), logPCA provides marginally better perplexity (except on the animal dataset where it performs worse than almost all other methods), and only given a suitable value of $K$. Its predictive performance can drastically decay with ill-chosen values of $K$. In contrast, our proposed methods do not require tuning $K$ to a proper value. Furthermore, they provide competitive prediction performance together with the interpretability of the decomposition. We have also compared against a standard LDA using the same range of $K$ than logPCA, but the perplexity was much worse (around six) for all datasets; we did not plot it due to the high difference in the scale.

6.4.3 Convergence of the variational inference algorithms

Fig. 11 displays the average perplexity values returned by the variational algorithms Beta-Dir VB, c-bICA-5 and bICA-5 along iterations. As expected, c-bICA tends to converge faster than bICA though not consistently so. Being initialized with a full tensor of dimension $K=100$ (as described in Section 6.2), Beta-Dir VB starts with a relatively higher perplexity but catches up with the two other methods in a reasonable number of iterations.

7 Conclusions

We have presented a unified view for Bayesian mean-parameterized NBMF. The interest of mean-parameterized models in NMF is that they keep factors interpretable since they belong to the same space than the observed data. We have addressed three models that correspond to three possible sets of constraints that each respect mean-parameterization. One model, Dir-Beta, is a Bayesian extension of the Aspect Bernoulli model of Bingham et al. (2009). Another model, Beta-Dir, corresponds to the binary ICA model of Kabán and Bingham (2008). We have proposed a new collapsed Gibbs sampler and a new collapsed variational inference method for estimation in these models. We have proposed a novel, third model, Dir-Dir, and we have designed a collapsed Gibbs sampler for inference with this model. Lastly, we have proposed a nonparametric extension for these three models. Experiments have shown that our nonparametric methods can achieve similar performance than the state-of-the-art methods applied with a suitable value of $K$. As expected, the more flexible logPCA can achieve better data approximation and in some cases prediction, but at the cost of interpretation which of utter importance in some applications.

Acknowledgements

This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 681839 (project FACTORY).

Appendix A Probability distributions functions

A.1 Bernoulli distribution

Distribution over a binary variable $x\in\{0,1\}$, with mean parameter $\mu\in[0,1]$:

[TABLE]

A.2 Beta distribution

Distribution over a continuous variable $x\in[0,1]$, with shape parameters $a>0$, $b>0$:

[TABLE]

A.3 Gamma distribution

Distribution for a continuous variable $x>0$, with shape parameter $a>0$ and rate parameter $b>0$:

[TABLE]

A.4 Dirichlet distribution

Distribution for $K$ continuous variables $x_{k}\in[0,1]$ such that $\sum_{k}x_{k}=1$. Governed by $K$ shape parameters $\alpha_{1},...\alpha_{K}$ such that $\alpha_{k}>0$:

[TABLE]

A.5 Discrete distribution

Distribution for the discrete variable ${\mathbf{x}}\in\{{\mathbf{e}}_{1},\ldots,{\mathbf{e}}_{K}\}$, where ${\mathbf{e}}_{i}$ is the $i^{th}$ canonical vector. Governed by the discrete probabilities $\mu_{1},...,\mu_{K}$ such that $\mu_{k}\in[0,1]$ and $\sum_{k}\mu_{k}=1$:

[TABLE]

The probability mass function can be written as:

[TABLE]

We may write $\text{Discrete}(\boldsymbol{\mathbf{x}}|\boldsymbol{\mathbf{\mu}})=\text{Multinomial}(\boldsymbol{\mathbf{x}}|1,\boldsymbol{\mathbf{\mu}})$.

A.6 Multinomial distribution

Distribution for an integer-valued vector $\boldsymbol{\mathbf{x}}=[x_{1},...,x_{K}]^{T}\in\mathbb{N}^{K}$. Governed by the total number $L=\sum_{k}x_{k}$ of events assigned to $K$ bins and the probabilities $\mu_{k}$ of being assigned to bin $k$:

[TABLE]

Appendix B Derivations for the Beta-Dir model

B.1 Marginalizing out $\boldsymbol{\mathbf{W}}$ and $\boldsymbol{\mathbf{H}}$ from the joint likelihood

We seek to compute marginal joint probability introduced in Eq. (22) and given by:

[TABLE]

Using the expression of the normalization constant of the Dirichlet distribution, the first integral can be computed as follows:

[TABLE]

The second integral in Eq. (22) is computed as follows. In Eq. (80) we use that $p(v_{fn}|\boldsymbol{\mathbf{h}}_{n},\boldsymbol{\mathbf{z}}_{fn})=\text{Bernoulli}(v_{fn}|\prod_{k}h_{kn}^{z_{fkn}})=\prod_{k}\text{Bernoulli}(v_{fn}|h_{kn})^{z_{fkn}}$ (recall that $\mathbf{z}_{fn}$ is an indicator vector). In Eq. (83), we use the expression of the normalization constant of the Beta distribution.

[TABLE]

B.2 Conditional prior and posterior distributions of $\boldsymbol{\mathbf{z}}_{fn}$

Applying the Bayes rule, the conditional posterior of $\boldsymbol{\mathbf{z}}_{fn}$ is given by:

[TABLE]

The likelihood itself decomposes as $p(\boldsymbol{\mathbf{V}}|\boldsymbol{\mathbf{Z}})=\prod_{n}p(\boldsymbol{\mathbf{v}}_{n}|\boldsymbol{\mathbf{Z}}_{n})$ and we may ignore the terms that do not depend on $\boldsymbol{\mathbf{z}}_{fn}$. Using Eq. (24) and the identity $\Gamma(n+b)=\Gamma(n)n^{b}$ where $b$ is a binary variable, we may write:

[TABLE]

The conditional prior term is given by

[TABLE]

Using $p(\boldsymbol{\mathbf{Z}})=\prod_{f}\uline{\boldsymbol{\mathbf{Z}}}_{f}$ and Eq. (23) we have

[TABLE]

Using $\sum_{k}p(\boldsymbol{\mathbf{z}}_{fn}={\mathbf{e}}_{k}|\boldsymbol{\mathbf{Z}}_{\neg fn})=1$, a simple closed-form expression of $p(\boldsymbol{\mathbf{z}}_{fn}|\boldsymbol{\mathbf{Z}}_{\neg fn})$ is obtained as follows:

[TABLE]

Combining Eqs. (84), (89) and (94), we obtain

[TABLE]

Appendix C Alternative Gibbs sampler for the Dir-Dir model

In this appendix, we show how to derive an alternative Gibbs sampler based on a single augmentation, like in the Beta-Dir model. This is a conceptually interesting result, though it does not lead to an efficient implementation.

Likewise the Beta-Dir model, the Dir-Dir model can be augmented using the single indicator variables $\boldsymbol{\mathbf{z}}_{fn}$, as follows:

[TABLE]

Note that compared to Eqs. (15)-(18) only the prior on $\boldsymbol{\mathbf{h}}_{n}$ is changed.

Like in Beta-Dir, we seek in this appendix to derive a Gibbs sampler from the conditional probabilities $p(\boldsymbol{\mathbf{z}}_{fn}|\boldsymbol{\mathbf{Z}}_{\neg fn},\boldsymbol{\mathbf{V}})$ given by

[TABLE]

The conditional prior term is identical to that of Beta-Dir and given by

[TABLE]

Like in Beta-Dir, the likelihood term factorizes as $p(\boldsymbol{\mathbf{V}}|\boldsymbol{\mathbf{Z}})=\prod_{n}p(\boldsymbol{\mathbf{v}}_{n}|\boldsymbol{\mathbf{Z}}_{n})$, and we now derive the expression of $p(\boldsymbol{\mathbf{v}}_{n}|\boldsymbol{\mathbf{Z}}_{n})$. As compared to Beta-Dir, a major source of difficulty lies in the fact that $p({\mathbf{h}}_{n})$ does not fully factorize anymore because of the Dirichlet assumption (and in particular $\sum_{k}h_{kn}=1$). In the following, we use the multinomial theorem to obtain Eq. (107)444Many thanks to Xi’an (Christian Robert) for giving us the trick via StackExchange. and we use the expression of the normalization constant of the Dirichlet distribution to obtain Eq. (110):

[TABLE]

We conclude that, though available in closed form, the expression of $p(\boldsymbol{\mathbf{v}}_{n}|\boldsymbol{\mathbf{Z}}_{n})$ (and thus $p(\boldsymbol{\mathbf{z}}_{fn}|\boldsymbol{\mathbf{Z}}_{\neg fn})$) involves the computation of $K\prod_{k=1}^{K}B_{kn}$ terms involving binomial coefficients, which is impractical in typical problem dimensions.