Bayesian Beta-Bernoulli Process Sparse Coding with Deep Neural Networks
Several approximate inference methods have been proposed for deep discrete
latent variable models. However, non-parametric methods which have previously
been successfully employed for classical sparse coding models have largely been
unexplored in the context of deep models. We propose a non-parametric iterative
algorithm for learning discrete latent representations in such deep models.
Additionally, to learn scale invariant discrete features, we propose local data
scaling variables. Lastly, to encourage sparsity in our representations, we
propose a Beta-Bernoulli process prior on the latent factors. We evaluate our
spare coding model coupled with different likelihood models. We evaluate our
method across datasets with varying characteristics and compare our results to
current amortized approximate inference methods.
Arunesh Mittal, Kai Yang, Paul Sajda, John Paisley
# Bayesian Beta-Bernoulli Process Sparse Coding with Deep Neural Networks

###### Abstract

Several approximate inference methods have been proposed for deep discrete latent variable models. However, non-parametric methods which have previously been successfully employed for classical sparse coding models have largely been unexplored in the context of deep models. We propose a non-parametric iterative algorithm for learning discrete latent representations in such deep models. Additionally, to learn scale invariant discrete features, we propose local data scaling variables. Lastly, to encourage sparsity in our representations, we propose a Beta-Bernoulli process prior on the latent factors. We evaluate our spare coding model coupled with different likelihood models. We evaluate our method across datasets with varying characteristics and compare our results to current amortized approximate inference methods.

## 1 Introduction

Sparse coding (Olshausen and Field, 1996) is an unsupervised latent factor model that has been widely used to uncover sparse discrete latent structure from data. Unlike auto-encoders, where the encoder is a parametric model, the encoder in sparse coding is an optimization algorithm that searches for an optimal encoding \(\mathbf{z}^{*}=\operatorname*{arg\,max}_{\mathbf{z}}p(\mathbf{x},\mathbf{z};\theta)\), which maximizes the joint likelihood of the data \(\mathbf{x}\) and latent encodings \(\mathbf{z}\). An advantage of the non-parametric approach is that it decouples the encoder and decoder such that the generalization error rises entirely from the reconstruction error of the decoder. In such models, sparsity is encouraged in the latent encodings \(\mathbf{z}\) via a prior such as Laplace, Cauchy or factorized Student-t prior (Goodfellow et al., 2016). Sparse coding by optimizing the MAP objective with a Laplace prior allows one to use gradient optimization methods for inferring \(\mathbf{z}\). However, one major drawback using such priors is that the latent factors in the encoding are encouraged to remain close to zero, even when those factors are active, whereas, for inactive elements, under the prior distribution, a factor being exactly zero has zero probability (Goodfellow et al., 2012).

Variational Auto Encoders (Kingma and Welling, 2013, 2019) have been popular deep generative models employed to uncover lower dimensional latent structure in data. Despite the flexibility of the deep likelihood model \(p(\mathbf{x}\mid\mathbf{z})\), VAEs use a parametric encoder network for inferring the latent encoding \(\mathbf{z}\), and hence do not benefit from the same advantages as that of a non-parametric encoding model. In VAEs, the generalization error is linked to both the decoder and the encoder and is difficult to disentangle. In addition to the limitations of using a parametric network for inference, amortized variational inference using parametric neural networks has additional learning constraints due to the amortization and approximation gaps in the variational objective used to train VAEs (Cremer et al., 2018). In principle, a non-parametric encoding model with a MAP-EM optimization can perform better than neural net parameterized ammortized inference, as it does not suffer from the amortization gap or the variational approximation gap. This comes at the cost of losing posterior uncertainty estimates, however, this might be an acceptable trade-off given that the posterior uncertainty in deep generative models via ammortized approximate inference is poorly calibrated and is still an area of active research (Nalisnick et al., 2018). Additionally, utilizing the MAP estimates, we can still potentially approximate posterior uncertainty using a Laplace approximation (Ritter et al., 2018).

VAE models with discrete latent factors (Maddison et al., 2016, Jang et al., 2016), do not work well with continuous data likelihood models, as the discrete sparse latent factors have limited representational capacity, and are unable to adequately represent local scale variations across an entire dataset. In fact, often one desires that the latent encodings only encode underlying latent structure of the data that is invariant to local data point scale variations.

To address the aforementioned issues, we propose a generative model with local scaling variables that decouples the data scaling from the discrete latent representation. We utilize a Beta-Bernoulli process prior on the latent codes that allows us to learn sparse discrete latent factors. For inference
in this model, we propose a MAP-EM greedy pursuit algorithm. We expect the inferred latent codes with true zeroes to have a stronger regularizing effect than the above mentioned sparsity promoting priors, which is especially advantageous in deep generative models with flexible neural network parameterized likelihood models. The primary disadvantage of the non-parametric encoder is that it requires greater time to compute \(\mathbf{z}\) due to the iterative algorithm, however, since the Beta-Bernoulli prior encourages the encodings to be sparse, as training progresses, the time taken to encode each data point significantly decreases over training iterations.

We demonstrate the efficacy of our model by proposing three different instantiations of our general model. We evaluate our models on discrete and continuous data by examining the representational capacity of our model by measuring the data reconstruction error, as well as the sparsity of our learned representations. We compare our models to widely used VAE (Kingma and Welling, 2013) and its discrete variant the Gumbel Softmax VAE (Jang et al., 2016). Not only does our model perform better in terms of reconstruction errors, it also learns substantially sparser latent encodings.

## 2 Related Work

We briefly review the VAE model that has been widely used to learn latent representations. In the typical VAE generative model, \(\mathbf{z}_{n}\) is drawn from a Gaussian prior, then given \(\mathbf{z}_{n}\), \(\mathbf{x}_{n}\) is then drawn from a distribution parametrized by a deep neural network \(f_{\theta}(\cdot)\), which maps \(\mathbf{z}_{n}\) to the sufficient statistics of the likelihood function \(p(\mathbf{x}_{n}\mid\mathbf{z}_{n})\):

\[\mathbf{z}_{n} \sim p(\mathbf{z}_{n})\] \[\mathbf{x}_{n} \sim p(\mathbf{x}_{n}\mid f_{\theta}(\mathbf{z}_{n});\mathbf{\theta})\]

Inference in this model is then performed using variational inference, however, unlike free form optimization used with mean field variational inference (Jordan et al., 1999), the VAE models, parametrize the variational \(q(\mathbf{z}_{n};\mathbf{\phi})\) distribution also with a neural network, that maps the data \(\mathbf{x}_{n}\) to the sufficient statistics of the \(q(\cdot)\) distribution. Then the posterior inference is performed by optimizing the Evidence Lower Bound ELBO using gradient methods:

\[\ln p(\mathbf{x})\geq\mathrm{ELBO}=\sum_{n}\mathbb{E}_{q(\mathbf{z}_{n}|\mathbf{x}_{n}; \mathbf{\phi})}\left[\ln\frac{p(\mathbf{x}_{n},\mathbf{z}_{n};\mathbf{\theta})}{q(\mathbf{z}_{n}| \mathbf{x}_{n};\mathbf{\phi})}\right]\]

## 3 Beta-Bernoulli Generative Process

We propose the following generative model with Beta-Bernoulli process prior. Given observed data \(\mathbf{x}_{n}\), the corresponding latent encoding \(\mathbf{z}_{n}\) is drawn from a Bernoulli process (BeP) parameterized by a beta process (BP), where, the Bernoulli process prior over each of the \(k\) factors \(\mathbf{z}_{nk}\in\mathbf{z}_{n}\), is parameterized by \(\mathbf{\pi}_{k}\) drawn from a finite limit approximation to the beta process (Griffiths and Ghahramani, 2011; Paisley and Carin, 2009). Since \(\mathbf{z}_{nk}\) is drawn from \(\mathrm{Bern}(\mathbf{\pi}_{k})\), where \(\mathbf{\pi}_{k}\sim\mathrm{Beta}\left(\alpha\gamma/K,\alpha\left(1-\gamma\,K \right)\right)\), \(k\in\{1,\ldots,K\}\), the random measure \(G_{n}^{K}=\sum_{k=1}^{K}\mathbf{z}_{nk}\delta_{f_{\theta}(\mathbf{z}_{nk})}\), \(\lim_{K\rightarrow\infty}G_{n}^{K}\) converges to a Bernoulli process (Paisley and Jordan, 2016).

Then given a latent binary vector \(\mathbf{z}_{n}\in\{0,1\}^{K}\), the observed data point \(\mathbf{x}_{n}\) is drawn from an exponential family distribution with a local scaling factor \(\mathbf{\lambda}_{n}\), also drawn from an appropriate exponential family distribution (sections 3.2 to 3.1.1). The natural parameters of this data distribution are parametrized by a \(L\) layered neural network \(f_{\theta}(\cdot)\). The neural network \(f_{\theta}(\cdot)\), maps the binary latent code \(\mathbf{z}_{n}\in\{0,1\}^{K}\) to \(\mathbb{R}^{D}\). This corresponds to the following generative process:

\[\mathbf{\pi}_{k} \sim\mathrm{Beta}\left(\alpha(\gamma/K),\alpha\left(1-(\gamma/K) \right)\right)\] \[\mathbf{z}_{nk} \sim\mathrm{Bern}(\mathbf{\pi}_{k})\] \[\mathbf{\lambda}_{n} \sim\mathrm{ExpFam}(\mathbf{\phi})\] \[\mathbf{x}_{n} \sim\mathrm{ExpFam}(f_{\theta}(\mathbf{z}_{n});\mathbf{\lambda}_{n})\]

where, \(\mathbf{\pi}_{k}\) is the global prior on \(\mathbf{z}_{nk}\), which corresponds to the \(k^{th}\) dimension of the latent encoded vector \(\mathbf{z}_{n}\), and \(\mathbf{z}_{n}\) is the local latent encoding for the \(n^{th}\) data point \(\mathbf{x}_{n}\). The likelihood model for \(\mathbf{x}\) is parametrized by local parameters \(\{\mathbf{\lambda}_{n}\}_{n=1}^{N}\) and the global parameters \(\mathbf{\theta}\).

During inference, the Beta-Bernoulli process prior on \(\mathbf{z}\), encourages the model to learn sparse latent encodings. As we would like the binary encodings to be scale invariant, modeling local data point specific scale distribution \(p(\mathbf{\lambda}_{n})\), allows us to marginalize out the scale variations in data when when inferring the latent code \(\mathbf{z}_{n}\). We demonstrate the utility of this non-parametric encoding model by coupling the Beta-Bernoulli process sparse encoding prior with three distinct exponential family likelihood models in the following sections.

### Scale Invariant Models

Given two data points \(\mathbf{x}_{m}\) and \(\mathbf{x}_{n}\), where \(\mathbf{x}_{m}\) is just a scaled version of \(\mathbf{x}_{n}\), we would want these data points to have the same latent embedding \(\mathbf{z}\). To disentangle the scale of the data points from the latent discrete representation, we introduce a local scale distribution for the Gaussian and Poisson likelihood models.

#### 3.1.1 GaussBPE

For real valued data, we use a Gaussian likelihood model, where \(f_{\theta}(\cdot)\) parametrizes the mean of the Gaussian distribution. We model the local data point scale \(\mathbf{\lambda}_{n}\) with a univariate Gaussian:

\[\mathbf{\lambda}_{n} \sim\mathcal{N}(0,c)\] \[\mathbf{x}_{n} \sim\mathcal{N}(\mathbf{\lambda}_{n}f_{\theta}(\mathbf{z}_{n}),\sigma^{2}I)\]

The likelihood given the local encoding \(\mathbf{z}_{n}\), depends on both local parameters \(\mathbf{\lambda}_{n}\) and the global neural net parameters
\(\mathbf{\theta}\). For \(\mathbf{\lambda}=1\), this is equivalent to the isotropic Gaussian likelihood with Gaussian prior generative model, employed by Kingma and Welling (2013).

#### 3.1.2 PoissBPE

For count data we use a Poisson likelihood model, where \(f_{\theta}(\cdot)\) parametrizes the rate of the Poisson distribution. We model the local data point rate \(\mathbf{\lambda}_{n}\) with a Gamma distribution. Additionally, we introduce a global parameter \(\mathbf{\beta}\):

\[\mathbf{\lambda}_{n} \sim\mathrm{Gamma}(a,b)\] \[\mathbf{x}_{n} \sim\mathrm{Poiss}(\mathbf{\lambda}_{n}\mathbf{\beta}f_{\theta}(\mathbf{z}_{ n}))\]

The likelihood then depends on local parameters \(\{\mathbf{\lambda}_{n}\}_{n=1}^{N}\) and the global neural net parameters \(\mathbf{\varphi}\), and \(\mathbf{\beta}\), a \(W\times T\) matrix, where each column \(\mathbf{\beta}_{,t}\in\Delta_{W-1}\). The global parameters include both the neural net parameters \(\mathbf{\varphi}\) and \(\mathbf{\beta}\), \(\mathbf{\theta}=\{\mathbf{\varphi},\mathbf{\beta}\}\)

In the context of topic modeling, \(W\) corresponds to number of words in the vocabulary, \(T\) corresponds to the number of topics, and \(\mathbf{\beta}_{:,t}\), corresponds to the \(t^{th}\) topic distrubution over words. Then \(\mathbf{x}_{n}^{(d)}\) is the number of occurances of word \(d\) in the \(n^{th}\) document.

### BernBPE

To evaluate a likelihood model, where we do not need to explicitly model the local scale, such a binary data, we use a Bernoulli likelihood model, where \(f_{\theta}(\cdot)\) parametrizes the mean of the Bernoulli distribution, without any local scaling variables:

\[\mathbf{x}_{n}\sim\mathrm{Bern}(f_{\theta}(\mathbf{z}_{n}))\]

Given the local encoding \(\mathbf{z}_{n}\), the likelihood model only depends on the global neural net parameters \(\mathbf{\theta}\). This model is equivalent to the Bernoulli likelihood model with Bernoulli prior employed by Jang et al. (2016).

## 4 Inference

We propose a MAP-EM algorithm to perform inference in this model. We compute point estimates for local latent encodings \(\{\mathbf{z}_{n}\}_{n=1}^{N}\) and the global parameters \(\mathbf{\theta}\) and compute posterior distributions over \(\mathbf{\pi}\) and \(\mathbf{\lambda}\). Since, \(\mathbf{\pi}\perp\mathbf{x}_{n}\mid\mathbf{z}_{n}\), utilizing the conjugacy in the model, we can analytically compute the conditional posterior \(q(\mathbf{\pi})\). Similarly, for the local scaling variables in the Gaussian likelihood and Poisson likelihood models, we can analytically compute the conditional posterior \(q(\mathbf{\lambda}_{n})\).

### Inference for local scale parameters

#### 4.1.1 \(q(\mathbf{\lambda})\) for Gaussian Likelihood

Since the conditional posterior \(q(\mathbf{\lambda})\triangleq p(\mathbf{\lambda}\mid\mathbf{x},\mathbf{\theta},\mathbf{z})\) factorizes as \(\prod_{n}q(\mathbf{\lambda}_{n})=\prod_{n}p(\mathbf{\lambda}_{n}\mid\mathbf{x}_{n},\mathbf{ \theta},\mathbf{z}_{n})\) and the posterior distribution over \(\mathbf{\lambda}_{n}\) is also a Gaussian, we can analytically compute the posterior \(q(\mathbf{\lambda}_{n})\):

\[q(\mathbf{\lambda}_{n}) =\mathcal{N}(\mathbf{\lambda}_{n}\mid\mu_{\mathbf{\lambda}_{n}|\mathbf{x}_{n },\mathbf{z}_{n},\mathbf{\theta}},\sigma^{2}_{\mathbf{\lambda}_{n}|\mathbf{x}_{n},\mathbf{z}_{n}, \mathbf{\theta}}) \tag{1}\] \[\sigma^{2}_{\mathbf{\lambda}_{n}|\mathbf{x}_{n},\mathbf{z}_{n},\mathbf{\theta}} =\big{(}c^{-1}+f_{\theta}(\mathbf{z}_{n})^{\top}f_{\theta}(\mathbf{z}_{n })/\sigma^{2}\big{)}^{-1}\] \[\mu_{\mathbf{\lambda}_{n}|\mathbf{x}_{n},\mathbf{z}_{n},\mathbf{\theta}} =(\sigma^{2}_{\mathbf{\lambda}_{n}|\mathbf{x}_{n},\mathbf{z}_{n},\mathbf{\theta}}) (f_{\theta}(\mathbf{z}_{n})^{\top}\mathbf{x}_{n})/\sigma^{2}\]

### \(q(\mathbf{\lambda})\) for Poisson Likelihood

The conditional posterior \(q(\mathbf{\lambda})\triangleq p(\mathbf{\lambda}\mid\mathbf{x},\mathbf{z},\mathbf{\theta},\mathbf{ \beta})\) factorizes as \(\prod_{n}q(\mathbf{\lambda}_{n})=\prod_{n}p(\mathbf{\lambda}_{n}\mid\mathbf{x}_{n},\mathbf{z }_{n},\mathbf{\theta},\mathbf{\beta})\). Given the Gamma prior on \(\mathbf{\lambda}_{n}\), the posterior distribution over \(\mathbf{\lambda}_{n}\) is also Gamma distributed, hence, we can analytically compute the posterior \(q(\mathbf{\lambda}_{n})\):

\[q(\mathbf{\lambda}_{n}) \propto\mathbf{\lambda}_{n}^{a-1}e^{-b\mathbf{\lambda}_{n}}\prod_{w}(\bm {\lambda}_{n})^{\mathbf{x}_{n}^{(w)}}e^{(-\mathbf{\lambda}_{n}\mathbf{\phi}_{w}^{(w)})} \tag{2}\] \[=\mathrm{Gamma}\left(\sum_{w}\mathbf{x}_{n}^{(w)}+a,b+1\right)\]

Where \(\mathbf{\phi}_{n}^{(d)}\triangleq[\mathbf{\beta}f_{\theta}(\mathbf{z}_{n})]_{d}\). Since \(\sum_{t}[f_{\theta}(\mathbf{z}_{n})]_{t}\triangleq 1\) and \(\sum_{w}\mathbf{\beta}_{w,t}\triangleq 1\), the sum over the random vector,

Figure 1: Greedy pursuit for \(\mathbf{z}\). For each greedy sparse coding step, first all bits in the active set \(\mathbf{\Omega}\) are turned on, then individually all \(j\in\{0,\dots,K\}\setminus\mathbf{\Omega}\) are turned on. The bit \(j^{*}\) that leads to maximal increase in the bound \(\mathcal{L}(\mathbf{z}_{\mathbf{\Omega}\cup\{j^{*}\}})\) is added to the active set \(\mathbf{\Omega}\). If adding a \(j^{*}\) leads to a decrease in the bound, the search is terminated and the sparse vector \(\mathbf{z}_{\mathbf{\Omega}}\) is returned as the sparse code. Here \(\mathcal{L}(\mathbf{z}^{*})\) is the optimal encoding that could be recovered if an exhaustive search were performed over \(2^{K}\) possible codes.
\(1\). Hence, the posterior \(q(\mathbf{\lambda}_{n})\) does not depend on \(f_{\theta}(\mathbf{z}_{n})\) or \(\mathbf{\beta}\). In practice, we only need to compute \(q(\mathbf{\lambda})\) for the entire dataset just once during training.

### Inference for latent variables

#### 4.3.1 Stochastic update for \(q(\mathbf{\pi})\)

For scalable inference, given a batch of data \(\{\mathbf{x}_{n}\}_{n\in\mathcal{B}}\), we first compute the latent codes \(\{\mathbf{z}_{n}\}_{n\in\mathcal{B}}\), then we can efficiently compute the posterior \(q(\mathbf{\pi}\mid\{\mathbf{z}_{n}\}_{n\in\mathcal{B}})\) using natural gradient updates (Hoffman et al., 2013). This posterior parameter update is a stochastic gradient step along the natural gradient, which is equivalent to following updates to the posterior sufficient statistics \(\{\mathbf{a}_{k}\}_{k=1}^{K}\) and \(\{\mathbf{b}_{k}\}_{k=1}^{K}\) with step size \(\eta\):

\[q(\mathbf{\pi}) =\prod_{k}\operatorname{Beta}(\mathbf{\pi}_{k}\mid\mathbf{a}_{k},\mathbf{b}_ {k}) \tag{3}\] \[\mathbf{a}_{k}^{\prime} =\alpha\frac{\gamma}{K}+\frac{N}{|S|}\sum_{n\in S}\mathbf{z}_{k}^{(n)}\] \[\mathbf{b}_{k}^{\prime} =\alpha\left(1-(\gamma/K)\right)+\frac{N}{|S|}\sum_{n\in S} \left(1-\mathbf{z}_{k}^{(n)}\right)\] \[\mathbf{a}_{k} \leftarrow(1-\eta)\mathbf{a}_{k}+\eta\;\mathbf{a}_{k}^{\prime}\] \[\mathbf{b}_{k} \leftarrow(1-\eta)\mathbf{b}_{k}+\eta\;\mathbf{b}_{k}^{\prime}\]

### Greedy pursuit for \(\mathbf{z}\)

For each model we marginalize the local scale variables \(\mathbf{\lambda}\) and global \(\mathbf{\pi}\) to compute the complete data joint likelihood lower bounds \(\mathcal{L}_{\mathcal{G}},\mathcal{L}_{\mathcal{P}},\mathcal{L}_{\mathcal{B}}\) for Gaussian, Poisson and Bernoulli likelihood models respectively, which include terms that only depend on \(\mathbf{z}\):

\[\mathcal{L}_{\mathcal{G}}(\mathbf{z}) =\ln\int p(\mathbf{x},\mathbf{\lambda}\mid\mathbf{\theta},\mathbf{z})d\lambda+ \mathbb{E}_{q(\mathbf{\pi})}\left[\ln p(\mathbf{z}\mid\mathbf{\pi})\right]\] \[\mathcal{L}_{\mathcal{P}}(\mathbf{z}) =\mathbb{E}_{q(\mathbf{\lambda})}\left[\ln p(\mathbf{x},\mathbf{\lambda} \mid\mathbf{\theta},\mathbf{z})\right]+\mathbb{E}_{q(\mathbf{\pi})}\left[\ln p(\mathbf{z} \mid\mathbf{\pi})\right]\] \[\mathcal{L}_{\mathcal{B}}(\mathbf{z}) =\ln p(\mathbf{x}\mid\mathbf{\theta},\mathbf{z})+\mathbb{E}_{q(\mathbf{\pi})} \left[\ln p(\mathbf{z}\mid\mathbf{\pi})\right]\]

The expected log prior is given by:

\[\mathbb{E}_{q(\mathbf{\pi})}\left[\log p(\mathbf{z}_{n}\mid\mathbf{\pi})\right] =\sum_{k}\mathbf{z}_{nk}[\psi(\mathbf{a}_{k}-\psi(\mathbf{a}_{k}+\mathbf{b}_{k})]\] \[+(1-\mathbf{z}_{nk})[\psi(\mathbf{b}_{k})-\psi(\mathbf{a}_{k}+\mathbf{b}_{k})]\]

where \(\psi(\cdot)\) is the digamma function.

For the Gaussian likelihood model we can marginalize \(\mathbf{\lambda}_{n}\) when maximizing \(p(\mathbf{x}_{n},\mathbf{\lambda}_{n},\mathbf{z}_{n})\):

\[\mathcal{L}_{\mathcal{G}}(\mathbf{z}_{n})=\ln\int p(\mathbf{x}_{n},\mathbf{ \lambda}_{n}\mid\mathbf{\theta},\mathbf{z}_{n})d\mathbf{\lambda}_{n}+\] \[\mathbb{E}_{q(\mathbf{\pi})}\left[\ln p(\mathbf{z}_{n}|\mathbf{\pi})\right]\]

where the marginal log likelihood can be calculated:

\[\ln\int p(\mathbf{x}_{n},\mathbf{\lambda}_{n}\mid\mathbf{\theta},\mathbf{z}_{n })d\mathbf{\lambda}_{n}=\] \[-\frac{1}{2} \bigg{[}\ln\left(1+\frac{c}{\sigma^{2}}f_{\theta}(\mathbf{z}_{n})^{ \top}f_{\theta}(\mathbf{z}_{n})\right)\] \[+\mathbf{x}_{n}^{\top}\left(\sigma^{-2}I-\frac{\frac{1}{\sigma^{2}}f_{ \theta}(\mathbf{z}_{n})f_{\theta}(\mathbf{z}_{n})^{\top}}{e^{-1}\sigma^{2}+f_{\theta}( \mathbf{z}_{n})^{\top}f_{\theta}(\mathbf{z}_{n})}\right)\mathbf{x}_{n}\bigg{]}\]

For the Poisson likelihood model we compute the expectation:

\[\mathbb{E}_{q(\mathbf{\lambda})}[\ln p(\mathbf{x},\mathbf{\lambda}\mid\mathbf{ \theta},\mathbf{z})]=\] \[\mathbf{x}_{n}^{(w)}\ln\mathbf{\phi}_{n}^{(w)}-\frac{a+\sum_{w}\mathbf{x}_{n}^ {(w)}}{b+1}\mathbf{\phi}_{n}^{(w)}\]

To optimize \(\mathcal{L}_{(\cdot)}(\mathbf{z})\), we employ a greedy pursuit algorithm, which is similar to the matching pursuit used by K-SVD (Aharon et al., 2006). We use \(\mathbf{z}_{\Omega_{n}}\) to denote a \(k\)-vector, corresponding to the latent vector for the \(n^{\text{th}}\) data point, where, \(\forall j\in\mathbf{\Omega}_{n},\mathbf{z}_{nj}=1\) and \(\forall j\notin\mathbf{\Omega}_{n},\mathbf{z}_{nj}=0\). To compute the sparse code given a data point \(\mathbf{x}_{n}\), we start with an empty active set \(\mathbf{\Omega}_{n}\), then \(\forall j\in\{1,\dots,K\}\), we individually set each \(\mathbf{z}_{nj}=1\) to find \(j^{*}\in\{1,\dots,K\}\setminus\mathbf{\Omega}_{n}\) that maximizes \(\mathcal{L}_{(\cdot)}(\mathbf{z}_{\Omega_{n}\cup\{j^{*}\}})\). We compute the scores \(\zeta^{+}\triangleq\mathcal{L}_{(\cdot)}(z_{\Omega_{n}\cup\{j^{*}\}})\) and \(\zeta^{-}\triangleq\mathcal{L}_{(\cdot)}(z_{\Omega_{n}},\mathbf{\theta})\). We add \(j^{*}\) to \(\mathbf{\Omega}_{n}\) only if \(\zeta^{+}>\zeta^{-}\), this step is necessary because unlike matching pursuit, the neural net \(f_{\theta}(\cdot)\) is a non-linear mapping from \(\mathbf{z}_{\Omega_{n}}\), hence, adding \(j^{*}\) to \(\mathbf{\Omega}_{n}\) can decrease \(\mathcal{L}_{(\cdot)}(\mathbf{z}_{\Omega_{n}})\). For each \(\mathbf{x}_{n}\), we repeat the preceding greedy steps to sequentially add factors to \(\mathbf{\Omega}_{n}\) till \(\mathcal{L}_{(\cdot)}(\mathbf{z}_{\Omega_{n}})\) ceases to monotonically increase.

The expected log prior on \(\mathbf{z}\) imposes an approximate beta process penalty. Low probability factors learned through \(q(\mathbf{\pi})\) lead to negative scores, and hence eliminate latent factors, encouraging sparse encodings \(\mathbf{z}_{n}\). During optimization as \(q(\mathbf{\pi}_{k})\) for a given dimension \(k\) decreases, the likelihood
that the \(k^{th}\) dimension will be utilized to encode the data point also decreases. Consequently, as training progresses, this allows for speed up of the sparse coding routine over iterations.

### Update for \(\boldsymbol{\theta}\)

To update the global parameters \(\boldsymbol{\theta}\), for each model we marginalize the local scale variables \(\boldsymbol{\lambda}\) and global \(\boldsymbol{\pi}\) to compute the complete data joint likelihood lower bounds \(\widehat{\mathcal{L}}_{\mathcal{G}},\widehat{\mathcal{L}}_{\mathcal{P}}, \widehat{\mathcal{L}}_{\mathcal{B}}\) for Gaussian, Poisson and Bernoulli likelihood models respectively, which include terms that only depend on \(\boldsymbol{\theta}\):

\[\widehat{\mathcal{L}}_{\mathcal{G}}(\boldsymbol{\theta}) =\mathbb{E}_{q(\boldsymbol{\lambda})}\left[\ln p(\boldsymbol{x}, \boldsymbol{\lambda}\mid\boldsymbol{\theta},\boldsymbol{z})\right]\] \[\widehat{\mathcal{L}}_{\mathcal{P}}(\boldsymbol{\theta}) =\mathbb{E}_{q(\boldsymbol{\lambda})}\left[\ln p(\boldsymbol{x}, \boldsymbol{\lambda}\mid\boldsymbol{\theta},\boldsymbol{z})\right]\] \[\widehat{\mathcal{L}}_{\mathcal{B}}(\boldsymbol{\theta}) =\ln p(\boldsymbol{x}\mid\boldsymbol{\theta},\boldsymbol{z})\]

For Gaussian:

\[\mathbb{E}_{q(\boldsymbol{\lambda})}\left[\ln p(\boldsymbol{x}, \boldsymbol{\lambda}\mid\boldsymbol{\theta},\boldsymbol{z})\right]=\| \boldsymbol{x}_{n}-\mu_{\boldsymbol{\lambda}_{n}|\boldsymbol{x}_{n}, \boldsymbol{z}_{n},\boldsymbol{\theta}}\cdot f_{\theta}(\boldsymbol{z}_{n}) \|^{2}\] \[+\sigma_{\boldsymbol{\lambda}_{n}|\boldsymbol{x}_{n},\boldsymbol{ z}_{n},\boldsymbol{\theta}}^{2}\cdot f_{\theta}(\boldsymbol{z}_{n})^{\top}f_{ \theta}(\boldsymbol{z}_{n})+c\]

For Poisson and Bernoulli the likelihood is same as that in sparse coding step. We use stochastic optimization to update \(\boldsymbol{\theta}\) using ADAM Kingma and Ba (2014). First order gradient methods with moment estimates such as ADAM, can implicitly take into account the rate of change of natural parameters \((\boldsymbol{a}_{k},\boldsymbol{b}_{k})\) for \(q(\boldsymbol{\pi})\) when optimizing the neural net parameters. The full sparse coding algorithm is outlined in Algorithm 1.

## 5 Empirical study

We demonstrate the potential of our beta process sparse encoding models in a variety of settings. We evaluate the Gaussian likelihood Beta-Bernoulli Process Encoder (Gauss-BPE) on scaled MNIST (LeCun et al., 2010) and CIFAR-10 (Krizhevsky, 2009) datasets. The scaled MNIST dataset consists of MNIST images that are randomly scaled using a scaling factor sampled from \(\mathcal{U}(-\text{scale max},\text{scale max})\). We evaluate the BernBPE on MNIST data. To compare Gauss-BPE to Gaussian VAE and BernBPE to Gumbel-Softmax VAE, we compare the sparsity of the learned encodings, as well as the reconstruction error on held-out data. We utilize the following metrics:

Lastly, we present qualitative results for the PoissonBPE on 20-Newsgroup dataset (Joachims, 1996) to uncover latent distributions over topics.

Figure 2: (A) Left: samples from scaled MNIST data. Middle: the corresponding sparse codes (reshaped as matrix for visualization), and the reconstructions using GaussBPE. Right: reconstructions using the VAE. (B) The probability of activation of the most class discriminative latent dimensions for randomly chosen classes from the CIFAR-10 dataset. As expected, the CIFAR-10 dataset utilizes more latent factors relative to the simpler MNIST dataset. (C) Sorted mean activation probabilities for all latent dimensions, for MNIST and CIFAR-10 datasets. (D) Time duration per epoch during training for all models.
### Sparsity

We quantify the sparsity of the inferred latent encodings using the Hoyer extrinsic metric [10], which is \(0\) for a fully dense vector and \(1\) for a fully sparse vector. For a set of latent encodings \(\{\mathbf{z}_{n}\}_{n=1}^{N}\), the sparsity is defined as:

\[\mathrm{Sparsity}(\{\mathbf{z}_{n}\}_{n=1}^{N})=\frac{1}{N}\sum_{n} \mathrm{Hoyer}(\mathbf{z}_{n})\] \[\mathrm{Hoyer}(\mathbf{z}_{n})=\frac{\sqrt{K}-\|\mathbf{z}_{n}\|_{1}/\| \mathbf{z}_{n}\|_{2}}{\sqrt{K}-1}\in[0,1]\]

For the VAE models, we use the encoding means \(\{\mathbb{E}_{q(\mathbf{z}_{n}|\mathbf{z}_{n})}|\mathbf{z}_{n}|\}_{n=1}^{N}\) in lieu of \(\{\mathbf{z}_{n}\}_{n=1}^{N}\).

### Reconstruction Error

For GaussBPE and Gaussian likelihood VAE, we report the reconstruction mean squared error (MSE). For the Gaussian likelihood VAE, \(\mathbb{E}[\mathbf{\lambda}_{n}]=1\), and we use \(\mathbb{E}_{q(\mathbf{z}_{n}|\mathbf{x}_{n})}[f_{\theta}(\mathbf{z}_{n})]\) instead of \(f_{\theta}(\mathbf{z}_{n})\):

\[\mathrm{MSE}(\{\mathbf{x}_{n},\mathbf{z}_{n}\}_{n=1}^{N})=\frac{1}{N}\|\mathbf{x}_{n}- \mathbb{E}[\mathbf{\lambda}_{n}]f_{\theta}(\mathbf{z}_{n})\|^{2}\]

For the BernBPE and Bernoulli likelihood Gumbel Softmax VAE, we report the negative log likelihood (NLL):

\[\mathrm{NLL}(\{\mathbf{x}_{n},\mathbf{z}_{n}\}_{n=1}^{N})=-\frac{1}{N}\ln p(\mathbf{x}_{n} |\mathbf{z}_{n})\]

For the VAE models, we use the same recognition network architecture as the original papers. For the VAE likelihood models and the GaussBPE and BernBPE likelihood models, we use the same architecture as that used by the Gumbel Softmax VAE paper. Notably, the last layer is linear for Gaussian VAE, however, sigmoid for GaussBPE, as in our model, \(\mathbf{\lambda}_{n}\) decouples the scaling of individual data points. A summary of all the hyperparameters used for all models can be found in the supplementary material.

We evaluate the PoissBPE model on 20-Newsgroup data. We pre-process the data by removing headers, footers and quotes, as well as English stop words to get a \(512\) dimensional vocabulary. We then vectorize each document to a \(1142\) dimensional vector, where each dimension represents the number of occurrences of a particular word in the vocabulary, within the document. For the PoissBPE, we choose a \(W\times T\)\(\mathbf{\beta}\) matrix, with \(W=1142\) and \(T=15\), this corresponds to a topic model with \(15\) topics, where each topic vector \(\mathbf{\beta}_{,t}\), is a distribution over the \(1142\) words. The last layer non-linearity is a softmax, hence, the \(f_{\theta}(\mathbf{z}_{n})\), maps \(\mathbf{z}_{n}\) to a probability distribution over the \(T\) topics.

## 6 Results

On binary MNIST data, where the scale of the data points does not affect the latent encodings, we found the BernBPE model to be comparable to the Gumbel Softmax VAE in terms of reconstruction error, however, it does so by utilizing substantially fewer latent dimensions. For real valued MNIST data, the GaussBPE significantly outperformed the Gaussian likelihood VAE in terms of both the reconstruction error as well as the sparsity of the latent codes. For randomly scaled MNIST data, the relative improvement in sparsity was similar to the improvement observed over VAE

\begin{table}
\begin{tabular}{l c c c c c c c c} \hline \hline Dataset & \multicolumn{2}{c}{MNIST} & \multicolumn{2}{c}{Scaled MNIST} & \multicolumn{2}{c}{MNIST} & \multicolumn{2}{c}{CIFAR 10} \\ Model & GS-VAE & BernBPE & VAE & GaussBPE & VAE & GaussBPE & VAE & GaussBPE \\ \cline{2-10} NLL & **81.55** & 82.16 & MSE & 32.92 & **9.18** & 16.94 & **8.51** & 79.15 & **75.88** \\ Sparsity & 0.40 & **0.93** & Sparsity & 0.72 & **0.86** & 0.83 & **0.86** & 0.81 & **0.96** \\ \hline \hline \end{tabular}
\end{table}
Table 1: Comparison of reconstruction errors and latent code sparsity on held-out data.

Figure 3: Learned topic associations through sparse codes. Each node represents a topic. The words within each node are the most representative words chosen from the top 15 most probable words from that topic. Groups of nodes connected by edges denote topics activated by the same sparse code, where the size of the node is proportional to the probability of the topic.
on real valued MNIST data, however, the reconstruction error was markedly better. Lastly, on the CIFAR-10 dataset, the GaussBPE performed better than VAE in terms of reconstruction error and sparsity. We summarize our experimental results in Table 1.

## 7 Discussion

We evaluated our models across four datasets to explore the effects of the different variables we introduce in our generative model. On the binary MNIST data, where scale is not a factor, as expected, we observed similar performance in terms of reconstruction error, however, the Beta-Bernoulli process prior encouraged sparsity in the latent representation, which lead to our model to learn much sparser representaions. For real valued MNIST data, with variations in intensity across images, the local scaling variable allowed the model to learn sparser encodings, while also improving the reconstruction error. We further explored this effect by exaggerating the local scale variations by randomly perturbing the intensity of the MNIST images. As we expected, this lead to significant deterioration in image reconstructions by the VAE. Our explicit modeling of local variations decoupled the local data scaling from the encoding process, which allowed the model to learn scale invariant encodings, resulting in substantially improved performance over the VAE. On natural image datasets such as CIFAR-10, we expect more variation in image intensity relative to the standard MNIST dataset. Since our model performs well even under random perturbation local data scaling, we expected the GaussBPE to perform well on the CIFAR-10 dataset. As we had hoped, the our model learned sparser encodings while also improving reconstruction error on the CIFAR-10 dataset.
