{'fraction_non_alphanumeric': 0.036961456446472384, 'fraction_numerical': 0.017380138949395822, 'mean_word_length': 4.719868852459016, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Memory-Efficient Topic Modeling

As one of the simplest probabilistic topic modeling techniques, latent Dirichlet allocation (LDA) has found many important applications in text mining, computer vision and computational biology. Recent training algorithms for LDA can be interpreted within a unified message passing framework. However, message passing requires storing previous messages with a large amount of memory space, increasing linearly with the number of documents or the number of topics. Therefore, the high memory usage is often a major problem for topic modeling of massive corpora containing a large number of topics. To reduce the space complexity, we propose a novel algorithm without storing previous messages for training LDA: tiny belief propagation (TBP). The basic idea of TBP relates the message passing algorithms with the non-negative matrix factorization (NMF) algorithms, which absorb the message updating into the message passing process, and thus avoid storing previous messages. Experimental results on four large data sets confirm that TBP performs comparably well or even better than current state-of-the-art training algorithms for LDA but with a much less memory consumption. TBP can do topic modeling when massive corpora cannot fit in the computer memory, for example, extracting thematic topics from 7 GB PUBMED corpora on a common desktop computer with 2GB memory.


Introduction
Latent Dirichlet allocation (LDA) (Blei et al., 2003) is a three-layer hierarchical Bayesian model for probabilistic topic modeling, computer vision and computational biology (Blei, 2012). The collections of documents can be represented as a document-word co-occurrence matrix, where each element is the number of word count in the specific document. Modeling each document as a mixture topics and each topic as a mixture of vocabulary words, LDA assigns thematic labels to explain non-zero elements in the document-word matrix, segmenting observed words into several thematic groups called topics. From the joint probability of latent labels and observed words, existing training algorithms of LDA approximately infers the posterior probability of topic labels given observed words, and estimate multinomial parameters for document-specific topic proportions and topic distributions of vocabulary words. The time and space complexity of these training algorithms depends on the number of non-zero (N N Z) elements in the matrix.
Probabilistic topic modeling for massive corpora has attracted intense interests recently. This research line is motivated by increasingly common massive data sets, such as online distributed texts, images and videos. Extracting and analyzing the large number of topics from these massive data sets brings new challenges to current topic modeling algorithms, particularly in computation time and memory requirement. In this paper, we focus on reducing the memory usage of topic modeling for massive corpora, because the memory limitation prohibits running existing topic modeling algorithms. For example, when the document-word matrix has N N Z = 5 × 10 8 , existing training algorithms of LDA often requires allocating more than 12GBytes memory including space for data and parameters. Such a topic modeling task cannot be done on a common desktop computer with 2GB memory even if we can tolerate the slow speed of topic modeling.
Because computing the exact posterior of LDA is intractable, we must adopt approximate inference methods for training LDA. Modern approximate posterior inference algorithms for LDA fall broadly into three categories: variational Bayes (VB) (Blei et al., 2003), collapsed Gibbs sampling (GS) (Griffiths and Steyvers, 2004), and loopy belief propagation (BP) (Zeng et al., 2011). We may interpret these methods within a unified message passing framework (Bishop, 2006), which infers the approximate marginal posterior distribution of the topic label for each word called message. According to the expectation-maximization (EM) algorithm (Dempster et al., 1977), the local inferred messages are used to estimate the best multinomial parameters in LDA based on the maximum-likelihood (ML) criterion.
VB is a variational message passing algorithm (Winn and Bishop, 2005), which infers the message from a factorizable variational distribution to be close in Kullback-Leibler (KL) divergence to the joint distribution. The gap between variational and true joint distributions cause VB to use computationally expensive digamma functions, introducing biases and slowness in the message updating and passing process Zeng et al., 2011). GS is based on Markov chain Monte Carlo (MCMC) sampling process, whose stationary distribution is the desired joint distribution. GS usually updates its message using the sampled topic labels from previous messages, which does not keep all uncertainties of previous messages. In contrast, BP directly updates and passes the entire messages without sampling, and thus achieves a much higher topic modeling accuracy. Till now, BP is very competitive in both speed and accuracy for topic modeling (Zeng et al., 2011). Similar BP ideas have also been discussed as the zero-order approximation of the collapsed VB (CVB0) algorithm within the mean-field framework Asuncion, 2010).
However, the message passing techniques often require storing previous messages for updating and passing, which leads to the high memory usage increasing linearly with the number of documents or the number of topics. So, to save the memory usage, we propose a novel algorithm for training LDA: tiny belief propagation (TBP). The basic idea of TBP is inspired by the multiplicative update rules of non-negative matrix factorization (NMF) (Lee and Seung, 2001), which absorbs the message updating into passing process without storing previous messages. Extensive experiments demonstrate that TBP enjoys a significantly less memory usage for topic modeling of massive data sets, but achieves a comparable or even better topic modeling accuracy than VB, GS and BP. Moreover, the speed of TBP is very close to BP, which is currently the fastest batch learning algorithm for topic modeling (Zeng et al., 2011). We also extend the proposed TBP using the block optimization framework (Yu et al., 2010) to handle the case when data cannot fit in computer memory. For example, we extend TBP to extract 10 topics from 7GB PUBMED biomedical corpus using a desktop computer with 2GB memory.
There have been two straight-forward machine learning strategies to process large-scale data sets: online and parallel learning schemes. On the one hand, online topic modeling algorithms such as online VB (OVB) (Hoffman et al., 2010) read massive corpora as a data stream composed of multiple smaller mini-batches. Loading each smaller mini-batch into memory, OVB optimizes LDA within the online stochastic optimization framework, theoretically converging to the batch VB's objective function. But OVB still needs to store messages for each mini-batch. When the size of mini-batch is large, the space complexity of OVB is still higher than the batch training algorithm TBP. In addition, the best online topic modeling performance depends highly on several heuristic parameters including the mini-batch size. On the other hand, parallel topic modeling algorithms such as parallel GS (PGS) (Newman et al., 2009) use expensive parallel architectures with more physical memory. Indeed, PGS does not reduce the space complexity for training LDA, but it distributes massive corpora into P distributed computing units, and thus requires only 1/P memory usage as GS. By contrast, the proposed TBP can reduce the space complexity for batch training LDA on a common desktop computer. Notice that we may also develop much more efficient online and parallel topic modeling algorithms based on TBP in order for a significant speedup.
The rest paper is organized as follows. Section 2 compares VB, GS and BP for message passing, and analyzes their space complexity for training LDA. Section 3 proposes the TBP algorithm to reduce the space complexity of BP, and discusses TBP's relation with the multiplicative update rules of NMF. Section 4 shows extensive experiments on four realworld corpora. Finally, Section 5 draws conclusions and envisions future work.

The Message Passing Algorithms for Training LDA
LDA allocates a set of semantic topic labels, z = {z k w,d }, to explain non-zero elements in the document-word co-occurrence matrix x W ×D = {x w,d }, where 1 ≤ w ≤ W denotes the word index in the vocabulary, 1 ≤ d ≤ D denotes the document index in the corpus, and 1 ≤ k ≤ K denotes the topic index. Usually, the number of topics K is provided by users. The topic label satisfies z k w,d = {0, 1}, K k=1 z k w,d = 1. After inferring the topic labeling configuration over the document-word matrix, LDA estimates two matrices of multinomial parameters: topic distributions over the fixed vocabulary φ W ×K = {φ ·,k }, where θ ·,d is a Ktuple vector and φ ·,k is a W -tuple vector, satisfying k θ k,d = 1 and w φ w,k = 1. From a document-specific proportion θ ·,d , LDA independently generates a topic label z k ·,d = 1, which further combines φ ·,k to generate a word index w, forming the total number of observed word counts x w,d . Both multinomial vectors θ ·,d and φ ·,k are generated by two Dirichlet distributions with hyperparameters α and β. For simplicity, we consider the smoothed LDA with fixed symmetric hyperparameters provided by users (Griffiths and Steyvers, 2004). To illustrate the generative process, we refer the readers to the original three-layer graphical µ w,d,n (1) representation for LDA (Blei et al., 2003) and the two-layer factor graph for the collapsed LDA Zeng et al. (2011).
Recently, there have been three types of message passing algorithms for training LDA: GS, BP and VB. These message passing algorithms have space complexity as follows, Total memory usage = data memory + message memory + parameter memory, (1) where the data memory is used to store the input document-word matrix x W ×D , the message memory is allocated to store previous messages during passing, and the parameter memory is used to store two output parameter matrices φ W ×K and θ K×D . Because the input and output matrices of these algorithms are the same, we focus on comparing the message memory consumption among these message passing algorithms.

Collapsed Gibbs Sampling (GS)
After integrating out the multinomial parameters {φ, θ}, LDA becomes the collapsed LDA in the collapsed hidden variable space {z, α, β}. GS (Griffiths and Steyvers, 2004) is a Markov Chain Monte Carlo (MCMC) sampling technique to infer the marginal distribution or message, µ w,d,n (k) = p(z k w,d,n = 1), where 1 ≤ n ≤ x w,d is the word token index. The message update equation is where z k ·,d,−n = w z k w,d,−n , z k w,·,−n = d z k w,d,−n , and the notation −n denotes excluding the current topic label z k w,d,n . After normalizing the message k µ w,d,n (k) = 1, GS draws a random number u ∼ Uniform[0, 1] and checks which topic segment will be hit as shown in Fig. 1A, where K = 4 for example. If the topic index k = 3 is hit, then we assign z 3 w,d,n = 1. The sampled topic label will be used immediately to estimate the message for the next word token. If we view the sampled topic labels as particles, GS can be interpreted as a special case of non-parametric belief propagation (Sudderth et al., 2003), in which only particles rather than complete messages are updated and passed at each iteration. Eq.
(2) sweeps all word tokens for 1 ≤ t ≤ T training iterations until the convergence criterion is satisfied. To exclude the current topic label z k w,d,n in Eq.
(2), we need to store all topic labels, z k w,d,n = 1, ∀w, d, n, in memory for message passing. In a common 32-bit desktop computer, GS generally uses the integer type (4 bytes) for each topic label, so the approximate message memory in bytes can be estimated by where w,d x w,d is the total number of word tokens in the document-word matrix. For example, 7GB PUBMED corpus has 737, 869, 083 word tokens, occupying around 2.75GB message memory according to Eq. (3). Based on inferred topic configuration z k w,d,n over word tokens, the multinomial parameters can be estimated as follows, These equations look similar to Eq. (2) except including the current topic label z k w,d,n in both numerator and denominator.

Loopy Belief Propagation (BP)
Similar to GS, BP (Zeng et al., 2011) performs in the collapsed hidden variable space of LDA called collapsed LDA. The basic idea is to integrate out the multinomial parameters {θ, φ}, and infer the marginal posterior probability in the collapsed space {z, α, β}. The collapsed LDA can be represented by a factor graph, which facilitates the BP algorithm for approximate inference and parameter estimation. Unlike GS, BP infers messages, µ w,d (k) = p(z k w,d = 1), without sampling in order to keep all uncertainties of messages. The message update equation is where The notation −w and −d denote all word indices except w and all document indices except d. After normalizing k µ w,d (k) = 1, BP updates other messages iteratively. Fig. 1B illustrates the message passing in BP when K = 4, slightly different from GS in Fig. 1A. Eq. (6) differs from Eq. (2) in two aspects. First, BP infers messages based on word indices rather than word tokens. Second, BP updates and passes complete messages without sampling. In this sense, BP can be viewed as a soft version of GS. Obviously, such differences give Eq. (6) two advantages over Eq.
(2). First, it keeps all uncertainties of messages for high topic modeling accuracy. Second, it scans a total of N N Z word indices for message passing, which is significantly less than the total number of word tokens w,d x w,d in x. So, BP is often faster than GS by scanning a significantly less number of elements (N N Z ≪ w,d x w,d ) at each training iteration (Zeng et al., 2011). Eq. (6) scans N N Z in the document-word matrix for 1 ≤ t ≤ T training iterations until the convergence criterion is satisfied. However, BP has a higher space complexity than GS. Because BP excludes the current message µ w,d (k) in message update (6), it requires storing all K-tuple messages. In the widely-used 32-bit desktop computer, we generally use the double type (8 bytes) to store all messages with the memory occupancy in bytes, which increases linearly with the number of topics K. For example, 7GB PUBMED corpus has N N Z = 483, 450, 157. When K = 10, BP needs around 36GB for message passing.
Notice that when K is large, Eq. (7) is significantly higher than Eq.
(3). Based on the normalized messages, the multinomial parameters can be estimated by These equations look similar to Eq. (6) except including the current message µ w,d (k) in both numerator and denominator.

Variational Bayes (VB)
Unlike BP in the collapsed space, VB (Blei et al., 2003;Winn and Bishop, 2005) passes variational messages,μ w,d (k) =p(z k w,d = 1), derived from the approximate variational distributionp to the true joint distribution p by minimizing the KL divergence, KL(p||p). The variational message update equation is , and the notation exp and Ψ are exponential and digamma functions, respectively. After normalizing the variational message kμ w,d (k) = 1, VB passes this message to update other messages. There are two major differences between Eq. (10) and Eq. (6). First, Eq. (10) involves computationally expensive digamma functions. Second, it include the current variational messageμ w,d in the update equation. The digamma function significantly slows down VB, and also introduces bias in message passing Zeng et al., 2011). Fig. 1C shows the variational message passing in VB, where the dashed line illustrates that the variational message is derived from the variational distribution. Because VB also stores the variational messages for updating and passing, its space complexity is the same as BP in Eq. (7). Based on the normalized variational messages, VB estimates the multinomial parameters as These equations are almost the same as Eqs. (8) and (9) but using variational messages.

Synchronous and Asynchronous Message Passing
Message passing algorithms for LDA first randomly initialize messages, and then pass messages according to two schedules: the synchronous and the asynchronous update schedules (Tappen and Freeman, 2003). The synchronous message passing schedule uses all messages at t − 1 training iteration to update current messages at t training iteration, while the asynchronous schedule immediately uses the updated messages to update other remaining messages within the same t training iteration. Empirical results demonstrate that the asynchronous schedule is slightly more efficient than the synchronous schedule (Zeng et al., 2011) for topic modeling. However, the synchronous schedule is much easier to extend for parallel computation. GS is naturally an asynchronous message passing algorithm. The sampled topic label will immediately influence the topic sampling process at the next word token. Both synchronous and asynchronous schedules of BP work equally well in terms of topic modeling accuracy, but the asynchronous schedule converges slightly faster than the synchronous one (Elidan et al., 2006). VB is a synchronous variational message passing algorithm, updating messages at iteration t using messages at iteration t − 1.

Tiny Belief Propagation
In this section, we propose TBP to save the message memory and data memory usage of BP in section 2.2. Generally, the parameter memory of BP takes a relatively smaller space when the number of topics K is small. For example, as far as 7GB PUBMED data set is concerned (D = 8, 200, 000 and W = 141043), when K = 10, the parameter θ K×D occupies around 0.6GB memory, while the parameter φ W ×K occupies around 0.01GB memory. For simplicity, we assume that the parameter memory is enough for topic modeling.

Message Memory
The algorithmic contribution of TBP is to reduce the message memory of BP to almost zero during message passing process. Combining Eqs. (6), (8) and (9) yields the approximate message update equation, where the current message µ w,d (k) is added in both numerator and denominator in Eq. (6). Notice that such an approximation does not distort the message update very much because the message µ w,d (k) is significantly smaller than the aggregate of other messages in both numerator and denominator. Eq. (13) has the following intuitive explanation. If the wth word has a higher likelihood in the topic k and the topic k has a larger proportion in the dth document, then the topic k has a higher probability to be assigned to the element x w,d , i.e., z k w,d = 1. The normalized message can be written as the matrix operation, where (φθ) w,d is the element at {w, d} after matrix multiplication φθ. Within the probabilistic framework, LDA generates the word token at index {w, d} using the likelihood (φθ) w,d , which satisfies w (φθ) w,d = 1, so that w,d (φθ) w,d = D is a constant. Replacing the normalized messages by Eq. (14), we can re-write Eqs. (8) and (9) as where the denominators play normalization roles to constrain k θ k,d = 1, θ k,d ≥ 0 and w φ w,k = 1, φ w,k ≥ 0. We absorb the message update equation into the parameter estimation in Eqs. (15) and (16), so that we do not need to store the previous messages during message passing process. We refer to these matrix update algorithm as TBP.
If we discard the hyperparameters α and β in Eqs. (16) and (15), we find that these matrix update equations look similar to the following multiplicative update rules in nonnegative matrix factorization (NMF) (Lee and Seung, 2001), where the objective of NMF is to minimize the following divergence, under the constraints φ w,d ≥ 0 and θ k,d ≥ 0. First, Eqs. (16) and (15) are different from Eqs. (17) and (18) In conclusion, if we discard hyperparameters in Eqs. (15) and Eq. (16), the proposed TBP algorithm becomes a special NMF algorithm: x ≈ φθ, where TBP focuses only on approximating non-zero elements x w,d = 0 by φθ in terms of the KL divergence. Notice that the hyperparameters play smoothing roles in avoiding zeros in the factorized matrices in Eqs. (15) and Eq. (16), where zeros are major reasons for worse performance in predicting unseen words in the test set (Blei et al., 2003). Conventionally, different training algorithms for LDA can be fairly compared by the perplexity metric (Blei et al., 2003;Hoffman et al., 2010), which has been previously interpreted as the geometric mean of the likelihood in the probabilistic framework. Comparing (24) with (20), we find that the perplexity metric can be also interpreted as a KL divergence between the document-word matrix x and the multiplication of two factorized matrices φθ. Because the TBP algorithm directly minimizes the KL divergence (22), it often has a much lower predictive perplexity on unseen test data than both GS and VB algorithms in Section 2 for better topic modeling accuracy. This theoretical analysis has also been supported by extensive experiments in Section 4. Indeed, LDA is a full Bayesian counterpart of the probabilistic latent semantic analysis (PLSA) (Hofmann, 2001), which is equivalent to the NMF algorithm with the KL divergence (Gaussier and Goutte, 2005). Moreover, the inference objective functions between LDA and PLSA are very similar, and PLSA can be viewed a maximum-a-posteriori (MAP) estimated LDA model (Girolami and Kabán, 2003). For example, two recent studies Zeng et al., 2011) find that the CVB0 and the simplified BP algorithms for training LDA resemble the EM algorithm for training PLSA. Based on these previous works, it is a natural step to connect the NMF algorithms with those message passing algorithms for training LDA. More generally, we speculate that such intrinsic relations also exist between finite mixture models such as LDA and latent factor models such as NMF (Gershman and Blei, 2012). As an example, the NMF algorithm in theory has been recently justified to learn topic models such as LDA with a polynomial time (Arora et al., 2012). Notice that TBP and other NMF algorithms do not need to store previous messages within the message passing framework in Section 2, and thus save a lot of memory usage.
Based on (15) and (16), we implement two types of TBP algorithms: synchronous TBP (sTBP) and asynchronous TBP (aTBP), similar to the synchronous and asynchronous message passing algorithms in Section 2.4. Because the denominator of (16) is a constant, it does not influence the normalized message (14). So, we consider only the unnormalized θ during the matrix factorization. However, the denominator of (15) depends on k, so we use a K-tuple vector λ K to store the denominator, and use the unnormalized φ during the matrix factorization. The normalization can be easily performed by a simple division (φ w,k + β)/(λ k + W β). Fig. 2 shows the synchronous TBP (sTBP) algorithm. We use three temporary matriceŝ φ,θ,λ in Line 2 to store numerators of (15), (16) and denominator of (15) for synchronization. Form Line 4 to 9, we randomly initializeφ,θ,λ by rand(K), which generates a random integer k, 1 ≤ k ≤ K. At each training iteration t, 1 ≤ t ≤ T , we copy the temporary matrices to φ, θ, λ and clear the temporary matrices to zeros from Line 12 to 13. Then, for each non-zero element in the document-word matrix, we accumulate the numerators of (15), (16) and the denominator of (15) by the K-tuple message η k in the temporary matricesφ,θ,λ from Line 15 to 19. In the synchronous schedule, the update of elements in the factorized matrices does not influence other elements within each iteration t. Fig. 3 shows the asynchronous TBP (aTBP) algorithm. Unlike the sTBP algorithm, aTBP does not require temporary matricesφ,θ,λ. After the random initialization from Line 4 to 9, aTBP reduces the matrices φ, θ, λ in a certain proportion from Line 13 to 15, which can be compensated by the updated K-tuple message η k from Line 17 to 20. In the asynchronous schedule, the change of elements in the factorized matrices φ, θ will immediately influence the update of other elements. In anticipation, the asynchronous schedule is more efficient to pass the influence of the updated elements in matrices than the synchronous schedule. The sTBP and aTBP algorithms will iterate until the convergence condition is satisfied or the maximum iteration T is reached. The time complexity of TBP is O(N N Z × KT ), where N N Z is the number of non-zero elements in the document-word matrix, K is the number of topics and T is the number of training iterations. sTBP has the space complexity O(3 × N N Z + 2 × KW + 2 × KD), but aTBP has the space complexity O(3 × N N Z + KW + KD). Generally, we use 3 × N N Z to store data in the memory including indices of non-zero elements in the document-word matrix, and also use KW + KD memory to store matrices φ and θ. Because sTBP uses additional matricesφ,θ for synchronization, it uses 2 × KW + 2 × KD for all matrices.

Data Memory
When the corpus data is larger than the computer memory, traditional algorithms cannot train LDA due to the memory limitation. We assume that the hard disk is large enough to store the corpus file. Recently, reading data from hard disk into memory as blocks is a promising method (Yu et al., 2010) to handle such problems. We can extend the TBP algorithms in Figs. 2 and 3 to read the corpus file as blocks, and optimize each block sequentially. For example, we can read each document in the corpus file at one time into memory and perform the TBP algorithms to refine the matrices {φ, θ}. After scanning all documents in the corpus data file, TBP finishes one iteration of training in Figs. 2 and 3. Similarly, we can also store the matrices {φ, θ} in the file on the hard disk when they are larger than computer memory. In such cases, TBP consumes almost no memory to do topic modeling. Because loading data into memory requires additional time, TBP running on files is around twice slower than that running on memory. For example, for the 7GB PUBMED corpus and K = 10, aTBP requires 259.64 seconds to scan the whole data file on the hard disk, while it requires only 128.50 seconds to scan the entire data in the memory at each training iteration. Another choice is to extend TBP to the online learning (Hoffman et al., 2010), which partitions the whole corpus file into mini-batches and optimizes each minibatch after one look sequentially. Although some online topic modeling algorithms like OVB can converge to the objective of corresponding batch topic modeling algorithm, we find that the best topic modeling accuracy depends on several heuristic parameters including the mini-batch size (Hoffman et al., 2010). In contrast, TBP is a batch learning algorithm that can handle large data memory with better topic modeling accuracy. Reading block data from hard disk to memory can be also applied to both GS and VB algorithms for LDA.

Relationship to Previous Algorithms
The proposed TBP connects the training algorithm of LDA to the NMF algorithm with KL divergence. The intrinsic relation between probabilistic topic models (Hofmann, 2001;Blei et al., 2003) and NMF (Lee and Seung, 2001) have been extensively discussed in several previous works (Buntine, 2002;Gaussier and Goutte, 2005;Girolami and Kabán, 2003;Zeng et al., 2011). A more recent work shows that learning topic models by NMF has a polynomial time (Arora et al., 2012). Generally speaking, learning topic models can be formulated within the message passing framework in Section 2 based on the generalized expectation-maximization (EM) (Dempster et al., 1977) algorithm. The objective is to maximize the joint distribution of PLSA or LDA in two iterative steps. At the E-step, we approximately infer the marginal distribution of a topic label assigned to a word called message. At the M-step, based on the normalized messages, we estimate two multinomial parameters according to the maximum-likelihood criterion. The EM algorithm iterates until converges to the local optimum. On the other hand, the NMF algorithm with KL divergence has a probabilistic interpretation (Lee and Seung, 2001), which views the multiplication of two factorized matrices as the normalized probability distribution. Notice that the widely-used performance measure perplexity (Blei et al., 2003;Hoffman et al., 2010) for topic models follows exactly the same KL divergence in NMF, which implies that the NMF algorithm may achieve a lower perplexity in learning topic models. Therefore, connecting NMF with LDA may inspire more efficient algorithms to learn topic models. For example, in this paper, we show that the proposed TBP can avoid storing messages to reduce the memory usage. More generally, we speculate that finite mixture models and latent factor models (Gershman and Blei, 2012) may share similar learning techniques, which may inspire more efficient training algorithms to each other in the near future.

Experimental Results
Our experiments aim to confirm the less memory usage of TBP compared with the state-ofthe-art batch learning algorithms such as VB (Blei et al., 2003), GS (Griffiths and Steyvers, 2004) and BP (Zeng et al., 2011) algorithms, as well as online learning algorithm such as online VB (OVB) (Hoffman et al., 2010). We use four publicly available document data sets (Porteous et al., 2008;Hoffman et al., 2010): ENRON, NYTIMES, PUBMED and WIKI. Previous studies (Porteous et al., 2008) revealed that the topic modeling result is relatively insensitive to the total number of documents in the corpus. Because of the mem- We randomly partition each data set into halves with one for training set and the other for test set. The training perplexity (24) is calculated on the training set in 500 iterations. Usually, the training perplexity will decrease with the increase of number of training iterations. The algorithm often converges if the change of training perplexity at successive iterations is less than a predefined threshold. In our experiments, we set the threshold as one because the decrease of training perplexity is very small after satisfying this threshold. The predictive perplexity for the unseen test set is computed as follows ). On the training set, we estimate φ from the same random initialization after 500 iterations. For the test set, we randomly partition each document into 80% and 20% subsets. Fixing φ, we estimate θ on the 80% subset by training algorithms from the same random initialization after 500 iterations, and then calculate the predictive perplexity on the rest 20% subset, where x 20% w,d denotes word counts in the the 20% subset. The lower predictive perplexity represents a better generalization ability.

Comparison with Batch Learning Algorithms
We compare TBP with other batch learning algorithms such as GS, BP and VB. For all data sets, we fix the same hyperparameters as α = 2/K, β = 0.01 (Porteous et al., 2008). The CPU time per iteration is measured after sweeping the entire data set. We report the average CPU time per iteration after T = 500 iterations, which practically ensures that GS, BP and VB converge in terms of training perplexity. For a fair comparison, we use the same random initialization to examine all algorithms with 500 iterations. To repeat our experiments, we have made all source codes and data sets publicly available (Zeng, 2012). These experiments are run on the Sun Fire X4270 M2 server with two 6-core 3.46 GHz CPUs and 128 GB RAMs. Table 2 compares the the message memory usage during training. VB and BP consumes more than 1GB memory to message passing when K = 100. VB and BP even require more  Figure 4: Predictive perplexity as K = {100, 300, 500, 700, 900} on ENRON, NYTIMES, PUBMED and WIKI data sets. The notation 0.8x and 0.4x denote the predictive perplexity is multiplied by 0.8 and 0.4, respectively.
than 9GB for message passing when K = 900, because their message memory increases linearly with the number of topics K in Eq. (7). In contrast, GS needs only 10 ∼ 20MB memory for message passing. The advantage of GS is that its memory occupancy does not depend on the number of topics K in Eq. (3). Therefore, PGS (Newman et al., 2009) can handle the relatively large-scale data set containing thousands of topics without memory problems using the parallel architecture. However, PGS still requires message memory for message passing at the distributed computing unit. Clearly, both aTBP and sTBP do not need memory space to store previous messages, and thus save a lot of memory usage. This is a significant improvement especially compared with VB and BP algorithms. In conclusion, TBP is our first choice for batch topic modeling when memory is limited for topic modeling of massive corpora containing a large number of topics. Fig. 4 shows the topic modeling accuracy measured by the predictive perplexity on the unseen test set. The lower predictive perplexity implies a better topic modeling performance. Obviously, VB performs the worst among all batch learning algorithms with the highest predictive perplexity. For a better illustration, we multiply VB's perplexity by 0.8 on ENRON and NYTIMES, and by 0.4 on PUBMED data sets, respectively. Also, we find that VB shows an overfitting phenomenon, where the predictive perplexity increases with the increase of the number of topics K on all data sets. The basic reason is that VB optimizes an approximate variational distribution with the gap to the true distribution. When the number of topics is large, this gap cannot be ignored, leading to serious biases. We see that GS performs much better than VB on all data sets, because it theoretically approximates the true distribution by sampling techniques. BP always achieves a much lower predictive perplexity than GS, because it retains all uncertainty of messages without sampling. Both sTBP and aTBP perform equally well on ENRON and PUBMED data sets, which also achieve the lowest predictive perplexity among all batch training algorithms. However, BP outperforms both sTBP and aTBP on NYTIMES and WIKI data sets. Also, aTBP outperforms both sTBP and GS, while sTBP performs slightly worse than GS. Because aTBP has consistently better topic modeling accuracy than GS on all data sets, we advocate aTBP for topic modeling in limited memory. As we discussed in Section 3.1, BP/TBP has the lowest predictive perplexity mainly because it directly minimizes the KL divergence between x and φθ from the NMF perspective. Fig. 5 shows the CPU time per iteration of all algorithms. All these algorithms has a linear time complexity of K. VB is the most time-consuming because it involves complicated digamma function computation Zeng et al., 2011). For a better illustration, we multiply the VB's training time by 0.3. Although BP runs faster than GS when K is small (K ≤ 100) (Zeng et al., 2011), it is sometimes slower than GS when K is large (K ≥ 100), especially on ENRON and PUBMED data sets. The reason lies in that GS often randomly samples a topic label without visiting all K topics, while BP requires searching all K topics for the message update. When K is very large, this slight difference will be enlarged. sTBP runs as fast as BP in most cases, but aTBP runs slightly slower than both sTBP and BP. Comparing two algorithms in Figs. 2 and 3, we find that aTBP uses more division operations than sTBP at each training iteration, which accounts for aTBP's slowness. As a summary, TBP has a comparable topic modeling speed as GS and BP but with reduced memory usage. Fig. 6 shows the training perplexity as a function of training iterations. All algorithms converge to a fixed point given enough training iterations. On all data sets, VB usually uses 110 ∼ 170 iterations, GS uses around 400 ∼ 470 iterations, and BP/TBP uses 180 ∼ 230 iterations for convergence. Although the digamma function calculation is slow, it reduces the number of training iterations of VB to reach convergence. GS is a stochastic sampling method, and thus requires more iterations to approximate the true distribution. Because BP/TBP is a deterministic message passing method, it needs less iterations to achieve convergence than GS. Overall, BP/TBP consumes the least training time until convergence according to Figs. 5 and 6.  7 shows the top ten words of K = 10 topics extracted by VB (red), GS (blue), BP (black), aTBP (green) and sTBP (magenta). We see that most topics contain similar top ten words but with a different order. More formally, we can adopt subjective measures such as the word intrusion in topics and the topic intrusion in documents (Chang et al., 2009) to evaluate extracted topics. PUBMED is a biomedical corpus. According to our prior knowledge in biomedical domain, we find these topics are all meaningful. Under this condition, we advocate TBP for topic modeling with reduced memory requirements.

Comparison with Online Algorithms
We compare the topic modeling performance between TBP and the state-of-the-art online topic modeling algorithm OVB (Hoffman et al., 2010) 1 on a desktop computer with 2GB memory. The complete 7GB PUBMED data set (Porteous et al., 2008) contains a total of D = 820, 000, 000 documents with a vocabulary size W = 141, 043. Currently, only TBP and online topic modeling methods can handle 7GB data set using 2GB memory. OVB (Hoffman et al., 2010) uses the following default parameters: κ = 0.5, τ 0 = 1024, and the mini-batch size S = 1024. We randomly reserve 40, 000 documents as the test set, and use the remainder 8, 160, 000 documents as the training set. The number of topics K = 10. The hyperparameters α = 2/K = 0.05 and β = 0.01. Fig. 8 shows the predictive perplexity as a function of training time (seconds in log scale). OVB converges slower than TBP because it reads input data as a data stream, discarding each mini-batch sequentially after one look. Notice that, for each mini-batch, OVB still requires allocating message memory for computation. In contrast, TBP achieves a much lower perplexity using less memory usage and training time. There are two major reasons. First, TBP directly optimizes the perplexity in terms of the KL divergence in Eq (24), so that it can achieve a much lower perplexity than OVB. Second, OVB involves computationally expensive digamma functions which significantly slow down the speed. We see that sTBP is a bit faster than aTBP because it does not perform the division operation at each iteration (see Figs. 2 and 3). Because aTBP influences matrix factorization immediately after the matrix update, it converges at a slightly lower perplexity than sTBP. rat level effect day plasma concentration activity normal cell hormone level rat serum day plasma normal concentration liver effect animal level rat concentration plasma serum liver blood normal effect glucose rat level plasma effect serum concentration day normal liver animal level rat plasma effect day serum concentration normal administration animal Topic 4 effect concentration rat activity level transport glucose sodium acid muscle muscle activity nerve stimulation response neuron effect cat responses potential pressure muscle heart dog cardiac nerve activity response effect stimulation effect activity response muscle brain stimulation responses rat drug nerve effect response muscle brain activity stimulation responses nerve rat drug Topic 5 children patient data problem population test effect age factor program children test patient data problem subject population factor program hospital patient disease children age syndrome infant clinical normal year incidence children patient age infant subject population problem test data factor children test subject problem data population age factor program hospital Topic 6 patient pressure lesion normal clinical disease pulmonary coronary cardiac dog patient pressure heart cardiac pulmonary dog normal left coronary arterial test system data problem patient subject program new analysis clinical patient dog pressure renal heart cardiac normal blood pulmonary lung pressure patient dog cardiac heart pulmonary normal left right flow Topic 7 effect neuron cell activity response stimulation rat cat responses brain effect rat concentration drug level administration action insulin dose glucose effect rat drug level administration activity response treatment action dose activity concentration enzyme effect acid inhibition ph transport cell uptake activity acid concentration effect liver enzyme rat uptake inhibition transport Topic 8 concentration effect measurement system acid drug water solution determination values acid concentration ph solution temperature effect reaction compound degrees*c water concentration effect ph transport solution temperature membrane degrees*c rate water measurement system determination solution model water temperature analysis technique concentration solution temperature concentration measurement determination ph water model reaction degrees*c Topic 9 protein cell acid activity enzyme dna rna fraction ph isolated protein activity enzyme dna fraction acid binding rna synthesis enzymes activity protein acid enzyme fraction enzymes amino*acid binding synthesis isolated protein dna acid rna fraction mutant isolated molecular*weight synthesis chain protein dna binding fraction rna enzyme mutant activity acid molecular*weight Topic 10 patient serum test level antigen antibody antibodies infection normal effect cell tissue number normal surface human membrane chromosome electron*microscopy formation cell dna rna mutant number region growth nuclei normal structure cell tissue number normal surface layer large development membrane section antigen infection strain virus antibody test serum antibodies mice human Figure 7: Top ten words in ten topics extracted from the subset of PUBMED: VB (red), GS (blue), BP (black), aTBP (green) and sTBP (magenta). Most topics contain similar words but with a different order. The subjective measures (Chang et al., 2009) such as word intrusions in topics and topic intrusions in documents are comparable among different algorithms.

Conclusions
This paper has presented a novel tiny belief propagation (TBP) algorithm for training LDA with significantly reduced memory requirements. The TBP algorithm reduces the message memory required by conventional message passing algorithms including GS, BP and VB. We also discuss the intrinsic relation between the proposed TBP and NMF (Lee and Seung, 2001) with KL divergence. We find that TBP can be approximately viewed as a special NMF algorithm for minimizing the perplexity metric, which is a widely-used evaluation method for different training algorithms of LDA. In addition, we confirm the superior topic modeling accuracy of TBP in terms of predictive perplexity on extensive experiments. For example, when compared with the state-of-the-art online topic modeling algorithm OVB, the proposed TBP is faster and more accurate to extract 10 topics from 7GB PUBMED corpus using a desktop computer with 2GB memory. Recently, the NMF algorithm has been advocated to learn topic models such as LDA with a polynomial time (Arora et al., 2012). The proposed TBP algorithm also suggests that the NMF algorithms can be applied to training topic models like LDA with a high accuracy in terms of the perplexity metric. We hope that our results may inspire more and more NMF algorithms (Lee and Seung, 2001) to be extended to learn other complicated LDA-based topic models (Blei, 2012) in the near future.