Bayesian mean-parameterized nonnegative binary matrix factorization

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-parameterization 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.

Appendices

A Probability distributions functions

1.1 A.1 Bernoulli distribution

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

$$\begin{aligned} \text {Bernoulli}(x | \mu )&= \mu ^x(1-\mu )^{1-x}. \end{aligned}$$

(68)

1.2 A.2 Beta distribution

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

$$\begin{aligned} \text {Beta}(x | a,b)&= \frac{\Gamma (a+b)}{\Gamma (a)\Gamma (b)}x^{a-1}(1-x)^{b-1}. \end{aligned}$$

(69)

1.3 A.3 Gamma distribution

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

$$\begin{aligned} \text {Gamma}(x | a,b)&= \frac{b^a}{\Gamma (a)}x^{a-1} e^{-bx}. \end{aligned}$$

(70)

1.4 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$:

$$\begin{aligned} \text {Dirichlet}(\varvec{\mathbf {x}} | \varvec{\mathbf {\alpha }})&= \frac{\Gamma (\sum _k \alpha _k)}{\prod _k \Gamma (\alpha _k)}\prod _k x^ {\alpha _k-1}. \end{aligned}$$

(71)

1.5 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$:

$$\begin{aligned} p( {\mathbf {x}} = {\mathbf {e}} _{k}) = \mu _{k} \end{aligned}$$

(72)

The probability mass function can be written as:

$$\begin{aligned} \text {Discrete}(\varvec{\mathbf {x}} | \varvec{\mathbf {\mu }})&= \prod _k \mu _k^{x_k}. \end{aligned}$$

(73)

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

1.6 A.6 Multinomial distribution

Distribution for an integer-valued vector $\varvec{\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:

$$\begin{aligned} \text {Multinomial}(\varvec{\mathbf {x}} | L, \varvec{\mathbf {\mu }})&= \frac{L!}{x1!...x_K!}\prod _k \mu _k^{x_k}. \end{aligned}$$

(74)

B Derivations for the Beta-Dir model

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

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

$$\begin{aligned} p(\mathbf {V}, \mathbf {Z}) = \prod _f \overbrace{ \int p(\varvec{\mathbf {w}}_{f}) \prod _{n} p(\varvec{\mathbf {z}}_{fn} | \varvec{\mathbf {w}}_f) \,d\mathbf {w}_f }^{p(\underline{\varvec{\mathbf {Z}}}_{f}) } \prod _n \overbrace{ \int \prod _{k} p(h_{kn}) \prod _{f} p(v_{fn} | \varvec{\mathbf {h}}_n, \varvec{\mathbf {z}}_{fn}) \, d {\mathbf {h}} _{n} }^{p(\varvec{\mathbf {v}}_n | \varvec{\mathbf {Z}}_n) }. \end{aligned}$$

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

$$\begin{aligned} p(\underline{\varvec{\mathbf {Z}}}_{f})&= \int p(\varvec{\mathbf {w}}_{f}) \prod _{n} p(\varvec{\mathbf {z}}_{fn} | \varvec{\mathbf {w}}_f) \,d\mathbf {w}_{f}\end{aligned}$$

(75)

$$\begin{aligned}&= \int \frac{\Gamma (\sum _k \gamma _k)}{\prod _k \Gamma (\gamma _k)}\prod _k w_{fk}^{\gamma _k-1}\prod _n w_{fk}^{z_{fkn}} \,d\mathbf {w}_{f}\end{aligned}$$

(76)

$$\begin{aligned}&= \frac{\Gamma (\sum _k \gamma _k)}{\prod _k \Gamma (\gamma _k)} \int \prod _k w_{fk}^{\gamma _k + L_{fk}-1} \,d\mathbf {w}_{f} \end{aligned}$$

(77)

$$\begin{aligned}&= \frac{\Gamma (\sum _k \gamma _k)}{\prod _k \Gamma (\gamma _k)} \frac{\prod _k \Gamma (\gamma _k + L_{fk})}{\Gamma (\sum _k \gamma _k + L_{fk})}. \end{aligned}$$

(78)

The second integral in Eq. (22) is computed as follows. In Eq. (80) we use that $p(v_{fn} | \varvec{\mathbf {h}}_n, \varvec{\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.

$$\begin{aligned} p(\varvec{\mathbf {v}}_n | \varvec{\mathbf {Z}}_n)&= \int \prod _k p(h_{kn}) \prod _{f} p(v_{fn} | \varvec{\mathbf {h}}_n, \varvec{\mathbf {z}}_{fn}) \, {d \mathbf {h}_n} \end{aligned}$$

(79)

$$\begin{aligned}&= \int \prod _k \left[ \frac{\Gamma (\alpha _k + \beta _k)}{\Gamma (\alpha _k)\Gamma (\beta _k)} h_{kn}^{\alpha _k-1}(1-h_{kn})^{\beta _k-1} \right] \prod _{fk} \left[ h_{kn}^{v_{fn}}(1-h_{kn})^{1-v_{fn}}\right] ^{z_{fkn}} {d \mathbf {h}_n} \end{aligned}$$

(80)

$$\begin{aligned}&= \prod _k \int \frac{\Gamma (\alpha _k + \beta _k)}{\Gamma (\alpha _k)\Gamma (\beta _k)} h_{kn}^{\alpha _k-1}(1-h_{kn})^{\beta _k-1} \prod _{f} \left[ h_{kn}^{v_{fn}}(1-h_{kn})^{1-v_{fn}}\right] ^{z_{fkn}} dh_{kn} \end{aligned}$$

(81)

$$\begin{aligned}&= \prod _k \frac{\Gamma (\alpha _k + \beta _k)}{\Gamma (\alpha _k)\Gamma (\beta _k)} \int h_{kn}^{\alpha _k + A_{kn} -1}(1-h_{kn})^{\beta _k + B_{kn}-1} dh_{kn} \end{aligned}$$

(82)

$$\begin{aligned}&= \prod _k \frac{\Gamma (\alpha _k + \beta _k)}{\Gamma (\alpha _k)\Gamma (\beta _k)} \frac{\Gamma (\alpha _k + A_{kn})\Gamma (\beta _k + B_{kn})}{\Gamma (\alpha _k + \beta _k + M_{kn})}. \end{aligned}$$

(83)

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

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

$$\begin{aligned} p(\varvec{\mathbf {z}}_{fn} | \varvec{\mathbf {Z}}_{\lnot fn}, \varvec{\mathbf {V}}) \propto p(\varvec{\mathbf {V}} | \varvec{\mathbf {Z}})p(\varvec{\mathbf {z}}_{fn} | \varvec{\mathbf {Z}}_{\lnot fn}). \end{aligned}$$

(84)

The likelihood itself decomposes as $p(\varvec{\mathbf {V}} | \varvec{\mathbf {Z}}) = \prod _{n} p(\varvec{\mathbf {v}}_n | \varvec{\mathbf {Z}}_n)$ and we may ignore the terms that do not depend on $\varvec{\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:

$$\begin{aligned} p(\varvec{\mathbf {v}}_n | \varvec{\mathbf {Z}}_n)&= \prod _{k} \frac{\Gamma (\alpha _k + \beta _k)}{\Gamma (\alpha _k)\Gamma (\beta _k)} \frac{\Gamma (\alpha _k + A_{kn}) \Gamma (\beta _k + B_{kn}) }{\Gamma (\alpha _k + \beta _k + M_{kn})} \end{aligned}$$

(85)

$$\begin{aligned}&\propto \prod _{k} \frac{\Gamma (\alpha _k + A_{kn}) \Gamma (\beta _k + B_{kn}) }{\Gamma (\alpha _k + \beta _k + M_{kn})} \end{aligned}$$

(86)

$$\begin{aligned}&= \prod _{k} \frac{\Gamma (\alpha _k + A_{kn}^{\lnot fn} + z_{fkn}v_{fn}) \Gamma (\beta _k + B_{kn}^{\lnot fn} + z_{fkn}\bar{v}_{fn}) }{\Gamma (\alpha _k + \beta _k + M_{kn}^{\lnot fn} + z_{fkn})} \end{aligned}$$

(87)

$$\begin{aligned}&\propto \prod _k \frac{ \Gamma (\alpha _k + A_{kn}^{\lnot fn}) (\alpha _k + A_{kn}^{\lnot fn})^{z_{fkn}v_{fn}} \Gamma (\beta _k + B_{kn}^{\lnot fn}) (\beta _k + B_{kn}^{\lnot fn})^{z_{fkn}\bar{v}_{fn}} }{ \Gamma (\alpha _k + \beta _k + M_{kn}^{\lnot fn}) (\alpha _k + \beta _k + M_{kn}^{\lnot fn})^{z_{fkn}} } \end{aligned}$$

(88)

$$\begin{aligned}&\propto \prod _k \left[ \frac{ (\alpha _k + A_{kn}^{\lnot fn})^{v_{fn}} (\beta _k + B_{kn}^{\lnot fn})^{\bar{v}_{fn}} }{ (\alpha _k + \beta _k + M_{kn}^{\lnot fn}) } \right] ^{z_{fkn}}. \end{aligned}$$

(89)

The conditional prior term is given by

$$\begin{aligned} p(\varvec{\mathbf {z}}_{fn} | \varvec{\mathbf {Z}}_{\lnot fn}) = p(\varvec{\mathbf {Z}})/p(\varvec{\mathbf {Z}}_{\lnot fn}). \end{aligned}$$

(90)

Using $p(\varvec{\mathbf {Z}}) = \prod _{f} \underline{\varvec{\mathbf {Z}}}_{f}$ and Eq. (23) we have

$$\begin{aligned} p(\varvec{\mathbf {z}}_{fn} | \varvec{\mathbf {Z}}_{\lnot fn})&\propto p(\underline{\varvec{\mathbf {Z}}}_f) \end{aligned}$$

(91)

$$\begin{aligned}&\propto \prod _k\Gamma (\gamma _k + L_{kn}^{\lnot fn} + z_{fkn}) \end{aligned}$$

(92)

$$\begin{aligned}&= \prod _k \Gamma (\gamma _k + L_{kn}^{\lnot fn}) (\gamma _k + L_{kn}^{\lnot fn})^{z_{fkn}} \end{aligned}$$

(93)

$$\begin{aligned}&\propto \prod _k (\gamma _k + L_{kn}^{\lnot fn})^{z_{fkn}}. \end{aligned}$$

(94)

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

$$\begin{aligned} p(\varvec{\mathbf {z}}_{fn} = {\mathbf {e}} _{k} | \varvec{\mathbf {Z}}_{\lnot fn})&= \frac{\gamma _k + L_{kn}^{\lnot fn}}{\sum _{k} (\gamma _k + L_{kn}^{\lnot fn})} \end{aligned}$$

(95)

$$\begin{aligned}&= \frac{\gamma _k + L_{kn}^{\lnot fn}}{\sum _{k} \gamma _k + N-1}. \end{aligned}$$

(96)

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

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

(97)

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 $\varvec{\mathbf {z}}_{fn}$, as follows:

$$\begin{aligned} \varvec{\mathbf {h}}_{n}&\sim \text {Dirichlet}(\varvec{\mathbf {\eta }}) \end{aligned}$$

(98)

$$\begin{aligned} \varvec{\mathbf {w}}_{f}&\sim \text {Dirichlet}(\varvec{\mathbf {\gamma }}) \end{aligned}$$

(99)

$$\begin{aligned} \varvec{\mathbf {z}}_{fn} | \varvec{\mathbf {w}}_f&\sim \text {Discrete}(\varvec{\mathbf {w}}_f) \end{aligned}$$

(100)

$$\begin{aligned} v_{fn} | \varvec{\mathbf {h}}_{n}, \varvec{\mathbf {z}}_{fn}&\sim \text {Bernoulli}\left( \prod _k h_{kn}^{z_{fkn}}\right) \end{aligned}$$

(101)

Note that compared to Eqs. (15)–(18) only the prior on $\varvec{\mathbf {h}}_{n}$ is changed. Like in Beta-Dir, we seek in this appendix to derive a Gibbs sampler from the conditional probabilities $p(\varvec{\mathbf {z}}_{fn} | \varvec{\mathbf {Z}}_{\lnot fn}, \varvec{\mathbf {V}})$ given by

$$\begin{aligned} p(\varvec{\mathbf {z}}_{fn} | \varvec{\mathbf {Z}}_{\lnot fn}, \varvec{\mathbf {V}}) \propto p(\varvec{\mathbf {V}} | \varvec{\mathbf {Z}})p(\varvec{\mathbf {z}}_{fn} | \varvec{\mathbf {Z}}_{\lnot fn}). \end{aligned}$$

(102)

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

$$\begin{aligned} p(\varvec{\mathbf {z}}_{fn} | \varvec{\mathbf {Z}}_{\lnot fn}) \propto \prod _k (\gamma _k + L_{kn}^{\lnot fn})^{z_{fkn}}. \end{aligned}$$

(103)

Like in Beta-Dir, the likelihood term factorizes as $p(\varvec{\mathbf {V}} | \varvec{\mathbf {Z}}) = \prod _{n} p(\varvec{\mathbf {v}}_n | \varvec{\mathbf {Z}}_n)$, and we now derive the expression of $p(\varvec{\mathbf {v}}_n | \varvec{\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)⁴ and we use the expression of the normalization constant of the Dirichlet distribution to obtain Eq. (110):

$$\begin{aligned} p(\varvec{\mathbf {v}}_n | \varvec{\mathbf {Z}}_n)&= \int p(\varvec{\mathbf {h}}_n) \prod _{f} p(v_{fn} | \varvec{\mathbf {h}}_n, \varvec{\mathbf {z}}_{fn}) \,d\mathbf {h}_n\end{aligned}$$

(104)

$$\begin{aligned}&=\int \frac{\Gamma (\sum _k \eta _k)}{\prod _k \Gamma (\eta _k)} \prod _k h_{kn}^{\eta _k-1} \prod _{f}\prod _k \left[ h_{kn}^{v_{fn}} (1-h_{kn})^{1-v_{fn}}\right] ^{z_{fkn}} \,d\mathbf {h}_n\end{aligned}$$

(105)

$$\begin{aligned}&= \frac{\Gamma (\sum _k \eta _k)}{\prod _k \Gamma (\eta _k)} \int \prod _k h_{kn}^{\eta _n + A_{kn}-1} (1-h_{kn})^{B_{kn}} \,d\mathbf {h}_n \end{aligned}$$

(106)

$$\begin{aligned}&= \frac{\Gamma (\sum _k \eta _k)}{\prod _k \Gamma (\eta _k)} \int \prod _k h_{kn}^{\eta _n + A_{kn}-1} \sum _{j_k=0}^{B_{kn}} \left( {\begin{array}{c}B_{kn}\\ j_k\end{array}}\right) (-h_{kn})^{j_k} \,d\mathbf {h}_n \end{aligned}$$

(107)

$$\begin{aligned}&=\frac{\Gamma (\sum _k \eta _k)}{\prod _k \Gamma (\eta _k)} \int \sum _{j_1=0}^{B_{1n}} ... \sum _{j_K=0}^{B_{Kn}} \prod _k h_{kn}^{\eta _k + A_{kn}-1} \left( {\begin{array}{c}B_{kn}\\ j_k\end{array}}\right) (-h_{kn})^{j_k} \,d\mathbf {h}_n \end{aligned}$$

(108)

$$\begin{aligned}&= \frac{\Gamma (\sum _k \eta _k)}{\prod _k \Gamma (\eta _k)} \sum _{j_1=0}^{B_{1n}} ... \sum _{j_K=0}^{B_{Kn}} \prod _k (-1)^{j_k}\left( {\begin{array}{c}B_{kn}\\ j_k\end{array}}\right) \int \prod _k h_{kn}^{\eta _k + A_{kn} + j_k -1} \,d\mathbf {h}_n\end{aligned}$$

(109)

$$\begin{aligned}&= \frac{\Gamma (\sum _k \eta _k)}{\prod _k \Gamma (\eta _k)} \sum _{j_1=0}^{B_{1n}} ... \sum _{j_K=0}^{B_{Kn}} \prod _k (-1)^{j_k} \left( {\begin{array}{c}B_{kn}\\ j_k\end{array}}\right) \frac{\Gamma (\eta _k + A_{kn} + j_k)}{\Gamma (\sum _k \eta _k + A_{kn} + j_k)}. \end{aligned}$$

(110)

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