Revisiting the Probabilistic Latent Semantic Analysis: The Method, Its Extensions and Its Algorithms

1. Introduction

Informally, Information Retrieval (IR) can be defined as the methods to process information to construct collections of documents. Probabilistic latent semantic analysis (PLSA) was firstly formulated as an unsupervised IR technique. This method, also known as Probabilistic latent semantic indexing (PLSI), was introduced in conference proceedings [1,2]. The classical reference is unsupervised learning by probabilistic latent semantic analysis by Hofmann [3]. PLSA is based on the ideas of Latent semantic analysis (LSA) [4] and in fact is a probabilistic remake. LSA uses cross terms and documents of a corpus to obtain a count or a table of co-occurrences. Then, arranging frequencies in a matrix, the Singular value decomposition (SVD) space span is considered a set of latent variables and interpreted as the aspect model [5]. The PLSA uses the frequencies to decompose them as mixtures or aggregate Markov models [3], and adjusted with the Expectation maximization (EM) algorithm.

Although PLSA was formulated as an IR technique, it has been used for diverse purposes. PLSA’s versatility, clarity of results, and solid statistical properties have enabled a wide range of applications, in which the concepts of words and documents are assimilated into other discrete entities, thus enabling justification of the hypotheses on which PLSA relies. However, the PLSA has several problems: $\left(i\right)$ the nature of the data and the underlying hypotheses leads to a rigid model; $\left(ii\right)$ the iterative nature, based on the EM algorithm, has very slow convergence; and $\left(iii\right)$ probabilistic interpretation is lacking for latent variables. Those problems translate to uneven growth, as partly determined by algorithmic and computational advances. These limitations have prompted several reformulations and a myriad of algorithms, the development of related techniques, such as Latent Dirichlet allocation (LDA), and other studies focused on the relationship between PLSA and Non-negative matrix factorization (NMF).

Exists many surveys and review articles including PLSA as a technique for IR, as [6], with a classical perspective. Other works recompile this technique as an alternative to classifying opinions from Twitter [7], or a method to detect fake news [7]. However, few reviews have focused exclusively on PLSA. One such review, by Tian [8], focuses on semantic image analysis.

In this manuscript, we pay special attention to what has been written on PLSA, the extension of this method to less restrictive data structures than co-occurrences or contingency tables, the obtained results by modifying the underlying hypotheses, and the relationship with other techniques. These results makes the PLSA a fundamental character, providing a probabilistic interpretation of the SVD.

This article is structured to reflect this point of view. The Section 3 is dedicated to the solutions of the two PLSA formulations, and the extent of their use as a supervised and semi-supervised technique. Criticism introduced by LDA is the starting point of several contributions, examined in Section 4. Mixtures in which it decomposes the PLSA relate it in a natural way to the NMF. From this viewpoint, extensions are introduced in Section 5, and are considered to represent a qualitative leap. Extensions which are conceptually relevant according to the exposed criteria, are compiled in the Section 6 and are summarized in Table 1. The works dedicated to efficiency improvement are also discussed in the Section 7. Although no studies have provided definitive results, they illustrate the efforts to construct fast and reliable computational solutions. In addition, this article has attempted to use a consistent notation.

2. Literature Review

According to our bibliographic searches (Web of Science, Scopus, Arxiv, and Google Scholar), there are a relatively large number of articles based on the PLSA. It has successfully spaned many research areas, as shown in Table 2.

The widest field of applications is Engineering, which we consider separate from Information Engineering or Computer Sciences. The applications of PLSA to this field rely on the ability of PLSA to handle discrete entities as words. Some examples are [19] seeks communication with machines and proposes to use the method to obtain a distribution. It claims that it is advantageous concerning supervised methods.

Applications of Information Engineering include syntactic structure study, quickly examined in [20]. Collaborative filtering, in which user ratings construct suitable matrices to perform PLSA algorithms [21]; or speech recognition, introducing a score concatenation matrix [22,23], cybersecurity [24] and analysis of keywords from webs related to certain topics and sentiment analysis (involving a system of definitions on which the users’ opinions and other instances, analyzed as co-occurrences) [25,26].

We adopt the Tians’ criteria for PLSA image applications. Tian classifies contributions into three types: image annotation, image retrieval, and image classification. Image annotation involves generating textual words to describe the content [27], with diverse applications including internet image filtering [28]; image retrieval, consisting of a procedure of ranking images in a database according to their posterior probabilities of being relevant to infer which visual patterns describe each object [8]; pioneering work [29]; and use in clinical image diagnosis [30,31] or facial expression semantics [32]. Image classification [33] has also enabled pain recognition [34] or autonomous driving [35]. Several variants exist, such as co-regularized PLSA, for the recognition of observed images with different perspectives [36], among many others. For applications and developments in semantic image analysis, we refer readers to the Tian (2018) review and its selected references [8].

Examples in Life Sciences can be found in bioinformatics [37,38] identificating genomic sequences with documents and some classes of genotype characteristics, such as words; neurodegenerative diseases, identifying common and non-common symptoms [39] and biology, the work [40] is devoted to the nuclear prediction and localization of proteins.

Applications in foundamental sicences include geophysics [41], instrumentation [42], spectroscopy [43]. A comparative study of PLSA, latent Dirichlet allocation (LDA), and other techniques in the framework of spectroscopy, recently published is [44].

The PLSA is a fundamental method that provides probabilistic sense to the SVD. Limitations on the original formulation compromise this assertion. Also, it causes problems of various types (described in the 4 section). We pay special attention to the works involving reformulations and those that relate to other techniques. In this manuscript, we have selected works of a methodological nature that imply a methodological contribution. We do not consider the number of citations to each work, the impact of the publication, or the interpretative orthodoxy. In this way, we equate many publications with interpretability or computability criteria, especially those that are extensions of the method.

3. The Method: PLSA Formulas

The original formulation of the PLSA, according to [3], provides a probabilistic solution to the problem of extracting a set of $z_k(k=1,\cdots ,K)$ latent variables of a data frame $N(d_i,w_j)$, obtained from a corpus of $d_i(i=1,\dots ,m)$ documents when crossed with a thesaurus of $w_j(j=1,\dots ,n)$ words. The relative frequencies

$$\begin{array}{c}\hfill n(d_i,w_j)=\frac{N(d_i,w_j)}{{\sum }_{ij}N(d_i,w_j)}\end{array}$$

(1)

are estimated by the joint probability $P(d_i,w_j)$. A key idea in this method is decomposing this probabilistic approximation as the product of conditional distributions over a set of latent variables. After some manipulations, and using the Bayes rule,

$$\begin{array}{cccc}\hfill P(d_i,w_j)& =P\left(d_i\right)\sum _{k}P\left(w_j\right|z_k)P\left(z_k\right|d_i)\hfill & \hfill & \left(\mathrm{asymmetric}\mathrm{formulation}\right)\hfill \end{array}$$

(2)

$$\begin{array}{cccc}\hfill & =\sum _{k}P\left(z_k\right)P\left(w_j\right|z_k)P\left(d_i\right|z_k)\hfill & \hfill & \left(\mathrm{symmetric}\mathrm{formulation}\right)\hfill \end{array}$$

(3)

where $P\left(d_i\right)$ and $P\left(z_k\right)$ are probabilities of the document $d_i$ and the latent variable $z_k$, respectively. Formulas (2) and (3) are was called by Hofmann the asymmetric and symmetric formulations [13], or formulations I and II [45].

The discrete nature of the documents identifies each one with the probabilities of ${(d_1,\dots ,d_n)}^t$ over the latent variables, and justifies the postulation that the mixtures $P\left(d_i\right|z_k)$ are k independent identically distributed (iid) multinomials. Because the same occurs for the words, the objective is to determine the parameters $\theta$ and $\varphi$ such that the conditional probabilities $P\left(w_j\right|z_k)\sim \mathrm{Multinomial}\left({\theta }_{jk}\right)$ and $P\left(z_k\right|d_i)\sim \mathrm{Multinomial}\left({\varphi }_{ki}\right)$ for the asymmetric formulation (alternatively $P\left(w_j\right|z_k)\sim \mathrm{Multinomial}\left({\theta }_{jk}\right)$ and $P\left(d_i\right|z_k)\sim \mathrm{Multinomial}\left({\varphi }_{ik}\right)$ for the symmetric case), with no hypothesis regarding the number or distribution of $z_k$, which is a set of dummy variables with no probabilistic sense.

The adjustment of mixtures, given by Formulas (2) and (3), is the other key idea for obtaining a reliable probabilistic interpretation by maximizing the likelihood of the parameters. A method widely used for this purpose is the EM algorithm, which always converges [46]. The use of the EM algorithm is roughly equivalent to the problem of fitting $P(d_i,w_j)$ to $n(d_i,w_j)$, but ensuring a maximum likelihood estimation of the sufficient (not necessarily minimal) parameters $\theta$ and $\varphi$.

In fact, the EM algorithm is a consequence of Jensen inequality [47]. For a function Q such that

$$\begin{array}{c}\hfill Q\left(M\right(\theta \left)\right|\theta )\ge M(Q\left(\theta \right|\theta \left)\right)\end{array}$$

(4)

where M is a map, and in statistics usages is the expectation, usually written as E. Then, for the log-likelihood $\mathcal{L}$, $\mathcal{L}\left(M\right(\theta \left)\right)\ge M\left(\mathcal{L}\right(\theta \left)\right)$ occurs, defining a monotonically increasing sequence reaching the limit if $M\left(\theta \right)=\theta$. In the PLSA case, the parameters (which are not provided by the model in a closed manner) are the mixtures of relations (2) or (3).

The EM algorithm supposes two steps: expectation and maximization. Expectation (E-step) is computed on the log-likelihood

$$\begin{array}{c}\hfill \mathcal{L}=\sum _{ij}n(d_i,w_j)logP(d_i,w_j)\end{array}$$

(5)

and for parametrization (2) or (3) takes the forms

$$\begin{array}{cc}\hfill \mathcal{L}& =\sum _{ij}n(d_i,w_j)log\left\{P\left(d_i\right)\sum _{k}P\left(w_j\right|z_k)P\left(z_k\right|d_i)\right\}\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & =\sum _{ij}n(d_i,w_j)log\left\{\sum _{k}P\left(z_k\right)P\left(w_j\right|z_k)P\left(d_i\right|z_k)\right\}\hfill \end{array}$$

(7)

for the asymmetric and the symmetric cases, respectively. In both cases, the expectation of $\mathcal{L}$ is the posterior

$$\begin{array}{c}\hfill E\left(\mathcal{L}\right)=P\left(z_k\right|d_i,w_j)\end{array}$$

(8)

and after several manipulations

$$\begin{array}{c}\hfill P\left(z_k\right|d_i,w_j)=\frac{P(z_k,d_i,w_j)}{P(d_i,w_j)}\end{array}$$

(9)

The expressions for $E\left(\mathcal{L}\right)$ for both formulations are shown in the Table 3.

The calculation of expectation $E\left(\mathcal{L}\right)$ presents several complications related to the meaning of the primed index appearing in the formulas of Table 3. Interpretation requires consideration of the expression $P\left(z_k\right|d_i,w_j)$ of Formula (8). For computational purposes, the object supporting this data structure is an array containing the matrices with the estimates of $P(d_i,w_j)$, fixing for each one the values of $z_k$ . Then each of the elements of the array is a matrix taking the form

$$\begin{array}{c}\hfill {\left[P(d_i,w_j)\right]}_{ij{k}^{\prime }}=\mathrm{vec}\left[P(·|z_{k^{\prime }})\right]\mathrm{vec}{\left[P(·|z_{k^{\prime }})\right]}^t(k^{\prime }=1,2,\cdots )\end{array}$$

(10)

indicating the primed index that is fixed (it should be noticed that a vector multiplied by its transpose is a matrix. In this case are $k^{\prime }=k$ matrices). The vec notation has been used to better identify the scalar products of the vectors of probabilities $P(·|z_{k^{\prime }})$ obtained by varying $z_k$. The entire array is

$$\begin{array}{c}\hfill {\left[P(d_i,w_j)\right]}_{ijk}=\left[{[P(d_i,w_j)]}_{ij1}]\left|{[P(d_i,w_j)]}_{ij2}\right|\dots |{[P(d_i,w_j)]}_{ijK}\right]\end{array}$$

(11)

Maximization (M-step) uses Lagrange multipliers, the correspondent derivatives, to obtain the solutions maximizing probabilities after eliminating them. These solutions for each formulation yield the generative models for the figures shown in Figure 1.

The execution of adjustment of probabilities, in both formulations, involves selecting a value for k, initializing the distributions appearing in (2) or (3), and computing the E-step and M-step in an iterative process in which $P(d_i,w_j)$ is recalculated until a certain condition is achieved. Hofmann has noted that the iterative process can end when there are no changes in the qualitative inputs, a condition called early stop [3]. A detailed, accessible derivation of the PLSA formulas and an introductory discussion of the EM algorithm convergence can be found in [45].

Another point to consider is what PLSA solutions are. Table 3 provides the formulas leading to the solutions. In many cases, providing words or documents best identify each aspect or latent variable would be more appropriate. Then, the numerical values of the columns of the involved matrices are ordered, and the corresponding labels are substituted, thus revealing the most relevant items in the respective latent class. Because the type of the desired solution has no ambiguity in the respective context, it is rarely explicit but does not confuse which result to provide. As an example, we provide two cases related to image study. For classification purposes, qualitative solutions are more suitable, and the numerical solutions are more suitable for spatial co-occurrence analysis on image regions.

Example 1. 

An example provided by Hofmann is reproduced below to illustrate the sense of word rank for interpretation of "the 4 aspects most likely generate the word “segment,” derived from a K=128 aspect model of the CLUSTER document collection. The displayed word stems are the most probable words in the class-conditional distributions $P\left(w_j\right|z_k)$, from top to bottom in descending order" [3].

Aspect 1	Aspect 2	Aspect 3	Aspect 4
imag	video	region	speaker
SEGMENT	sequenc	contour	speech
color	motion	boundari	recogni
tissu	frame	descript	signal
Aspect1	scene	imag	train
brain	SEGMENT	SEGMENT	hmm
slice	shot	precis	sourc
cluster	imag	estim	speakerindepend
mri	cluster	pixel	SEGMENT
algorithm	visual	paramet	sound

In addition, we provide an artificial example to illustrate the effects of selection of K, consisting of a corpus of 5 ($d1$ to $d5$) documents containing letters $\{a,b,c,d,e,f\}$, which we assimilate into words in a thesaurus. The co-occurrences’ data frame N isand the frequency matrix n

$$\begin{array}{c}\hfill n(d_i,w_j)=\left(\begin{array}{cccccc}.086& .114& 0& 0& 0& 0\\ .086& .086& 0& 0& 0& 0\\ .029& .086& .114& .029& 0& 0\\ 0& 0& .057& .114& 0& 0\\ 0& 0& 0& 0& .086& .114\end{array}\right)\end{array}$$

If in this example, the objective is to classify documents by subject (or specialized words with the correspondent matters). Simple visual inspection indicates that they are 3. For the symmetric case formulas, running $p=1000$ iterations in each case, the results are

$$\begin{array}{cccccc}\hfill & \mathrm{for} k=2P\left(d_i\right|z_k)\hfill & \hfill & =\left(\begin{array}{cc}0& .411\\ 0& .588\\ .333& 0\\ .278& 0\\ .167& 0\\ .222& 0\end{array}\right)\hfill & \hfill & \left(\begin{array}{cc}c& b\\ d& a\\ f& -\\ e& -\\ -& -\\ -& -\end{array}\right)\hfill \\ \hfill & \mathrm{for} k=3P\left(d_i\right|z_k)\hfill & \hfill & =\left(\begin{array}{ccc}0& 0& .462\\ 0& 0& .538\\ .545& 0& 0\\ .454& 0& 0\\ 0& .429& 0\\ 0& .571& 0\end{array}\right)\hfill & \hfill & \left(\begin{array}{ccc}c& f& b\\ d& e& a\\ -& -& -\\ -& -& -\\ -& -& -\\ -& -& -\end{array}\right)\hfill \\ \hfill & \mathrm{for} k=5P\left(d_i\right|z_k)\hfill & \hfill & =\left(\begin{array}{ccccc}.007& 0& 0& .462& 0\\ .347& 0& 0& .538& 0\\ .646& .333& 0& 0& 0\\ 0& .667& 0& 0& 0.\\ 0& 0& .538& 0& .419\\ 0& 0& .462& 0& .581\end{array}\right)\hfill & \hfill & \left(\begin{array}{ccccc}c& d& e& b& f\\ b& c& f& a& e\\ a& -& -& -& -\\ -& -& -& -& -\\ -& -& -& -& -\\ -& -& -& -& -\end{array}\right)\hfill \end{array}$$

The characters’ matrices are the ordination of the most likely words identifying each latent variable (informally, the subjects in our toy example). Lines represent probabilities close to zero and are not useful for classification. The effect of selecting K is clear in the comparison of columns 3 and 5, which are equivalent (for $k=5$).

3.1. Training and Prediction

The PLSA algorithm can be executed for the entire data set, providing results in the same manner as probabilistic clustering methods [48, Chap. 3]. However, to exploit the predictive power of PLSA, the model must be fitted on the available data (or training phase). Predictions for new observations are made by simply comparing them with the trained data set.

In the prediction phase, cannot assign probabilities for documents that are not in the training phase, because non-zero probabilities are needed. This problem has been solved in [49] by splitting the data set into a training group with the $d_i$ observed documents and the new unobserved documents $q\in \mathcal{Q}$. By using probabilities $P\left(z_k\right|d_i)$ instead of $P\left(d_i\right|z_k)$ in (2) and expanding the logarithm, equation (6) can be rewritten as

$$\begin{array}{c}\hfill \mathcal{L}=\sum _{ij}n(d_i,w_j)logP\left(d_i\right)+\sum _{ij}n(d_i,w_j)logP\left(w_j\right|d_i)\end{array}$$

(12)

To avoid a zero probability of the unseen documents in the training phase, Brants has introduced $P\left(d_i\right)>0$, stating that the log-likelihood can be maximized, taking into account only the second term of (12), and

$$\begin{array}{c}\hfill P^{\left(new\right)}\left(Q\right)=\prod _{ij}P\left(w_j\right|q_i)\end{array}$$

(13)

Brants has highlighted that equation (13) does not represent the true likelihood, but if the goal is likelihood maximization, the same parameter setting is found as that when the true likelihood had been maximized [49]. The same article has proposed other methods for estimating likelihood on the basis of marginalization and splitting. Brants has also proposed PLSA folding-in, a more refined derivation of this technique [50]. A further improvement, which is more computationally efficient and is protected by a patent is [51], involves estimating the log-likelihood by spiting the data set in the training set, denoted $n^{\prime }(d_i,w_j)$, and introducing the unknown documents one by one as the second term of

$$\begin{array}{c}\hfill \mathcal{L}\propto \sum _{ij}{n}^{\prime }(d_i,w_j)logP\left(d\right)+\sum _{ij}logP\left(w_j\right|d_i)\end{array}$$

(14)

In the symmetric formulation, after training on the documents by using the formulas given in Table 3, new documents can be classified by simply alternating the expressions given by [38]

$$\begin{array}{cc}\hfill P\left(z_k\right|d_i,w_j)& =\frac{P\left(z_k\right|d_i)P\left(w_j\right|z_k)}{{\sum }_{k^{\prime }}P(z_k\prime |d_i)P\left(w_j\right|z_k\prime )}\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill P(z_k,d_i)& =\frac{{\sum }_{i}n(d_i,w_j)P\left(z_k\right|d_i,w_j)}{{\sum }_{ij}n(d_i,w_j)P\left(z_k\right|d_i,w_j)}\hfill \end{array}$$

(16)

In this case, binary data can be handled by entering a matrix $\mathbf{A}$ such that

$$\begin{array}{c}\hfill {\left[\mathbf{A}\right]}_{ij}=\left\{\begin{array}{c}1\mathrm{if}i\mathrm{is}\mathrm{annotated}\mathrm{to}j\hfill \\ 0\mathrm{otherwise}\hfill \end{array}\right.\end{array}$$

(17)

substituting $n(d_i,w_j)$ in equations of Table 3.

PLSA can also be used as a semi-supervised learning tool, in a process known as semi-supervised PLSA [52]. Using this mode requires entering labeled and non-labeled data in the EM iterative process, and being able to split the data set into a portion in which the labels are assigned and a portion in which the labels are not assigned. A measure of similarity performs the rest of the task. Another related strategy involves introducing the link functions must-link and cannot-link in the training phase [53].

4. Criticism: The LDA and Reformulations

The work of Hofmann is not a closed contribution. The rigidity of the distributional hypotheses that limit applicability has been exploited for different purposes. In addition, several studies have related Hofmann’s work to other techniques.

One of the first criticisms was noted by Blei, who has argued that Hofmann’s work is incomplete in that it provides no probabilistic model at the level of documents. This incompleteness leads to several problems: $\left(i\right)$ the number of parameters grows linearly with the size of the corpus, thus resulting in severe problems with overfitting, and $\left(ii\right)$ how to assign probabilities to a document outside the training set is unclear; LDA has been proposed to solve this problem [9].

4.1. Latent Dirichlet Allocation

LDA introduces a generative Bayesian model that maps documents on topics such that the words of each document are captured by these topics. Each document is described by a topic distribution, and each topic is described by a word distribution. Introducing $\theta$, a k dimensional Dirichlet with parameter ${\alpha }_k$, and $\beta$ as an array of initialization with values $P\left(w\right|z)$, and maintaining the notation of Formulas (2) and (3),

$$\begin{array}{c}\hfill P(\theta ,\mathbf{z},\mathbf{w}|\alpha ,\beta )=P\left(\theta \right|\alpha )\prod _{k}P\left(z_h\right|\theta )P\left(w_j\right|z_k,\beta )\end{array}$$

(18)

The LDA is also a generative model. The probabilities of a document and a corpus are obtained by marginalizing (integrating) over the complete collection. Further improvements to the model, also provided by Blei, include hierarchical LDA [9] and dynamic LDA [54].

The LDA is a closely related technique that is different from PLSA. Its criticisms provided a starting point for several developments. Formal equivalence with PLSA has been shown by [55] and has led to several proposed solutions to those problems in the case of the PLSA. Although LDA is not the objective of our review, we indicate exists further developments of this technique. We underline the works of Teh in which he proposes a non-parametric approach of mixture components with a hierarchical Bayesian distribution [56]. A hierarchical nested topic model is described in [57], and more recently in [58].

4.2. Other Formulations

There are reformulations of the PLSA, mainly arising from these criticisms. These developments have the objective of relaxing distributional hypotheses and overfitting problems.

4.2.1. Extension to Continuous Data

A generalization of the PLSA for continuously evaluated responses has also been provided by Hofmann, in the context of collaborative filtering, as an alternative to the neighbor regression method [10]. The method construction assumes a set of person items $y_j$ rated v for a set of persons $u_i$. Then,

$$\begin{array}{c}\hfill P\left(v\right|u,y)=\sum _{z}P\left(z\right|u)P(v;{\mu }_{yz},{\sigma }_{yz})\end{array}$$

(19)

where $\mu$ and $\sigma$ are the expectation and variance, respectively, and assuming $P\left(z\right|u\left)P\right(v;\mu ,\sigma )$

$$\begin{array}{c}\hfill P(v;u,\sigma )=\frac{1}{\sqrt{2\pi }\sigma }exp\left\{-\frac{{(v-\mu )}^2}{2{\sigma }^2}\right\}\end{array}$$

(20)

which is fitted with the EM algorithm.

Within the semantic image analysis field, the visual entities from a database are assimilated with the words from a thesaurus [59], but as discrete entities. This variant constitutes the Gaussian mixture model PLSA [60], and assumes a normal distribution of the descriptors $f_j$ (the most relevant visual words) such that $f_h\sim N\left(f_h\right|{\mu }_k,{\Sigma }_k)(h\le j)$. Horster has noted that this expression is difficult to train, and has proposed the alternative models shared Gaussian words PLSA and fixed shared Gaussian words PLSA. A more general treatment, in which normality is postulated for the mixtures $P\left(w_j\right|d_i)$, has been reported [61].

4.2.2. Randomized Probabilistic Latent Semantic Analysis

Randomized PLSA arose to address the problem of overfitting [62]. Taking a random fraction of the trained data sets, the method proceeds by folding the training data set $\mathcal{T}=\{T_1,\dots ,T_{\mathsf{\Omega }}\}$ and the fraction $T^l(T^l<\mathsf{\Omega })$ to run the PLSA algorithm with the l samples. The average of the results is the provided output.

The basis for this statement is the work by Ho [63] on the subspace method. This method takes random subsets of the support vector machine to avoid computational complexity. In addition, the derived algorithm has been reported to be slower than the conventional PLSA implementation.

4.2.3. Tensorial Approach

Non-negative tensor factorization was introduced by [64] for n-way data structures. Peng has established the relationship with the PLSA in [65], noting that the case $n=2$ corresponds to the NMF and illustrating that $n=3$ allows for handling more complex data structures. The objective is to better estimate the number of latent variables, or clueters.

Peng has introduced a structure of the type

$$\begin{array}{c}\hfill {\left[\underline{\mathbf{F}}\right]}_{ijl}\approx P(d_i,w_j,x_l)\end{array}$$

(21)

called a tensor, and being now $x_l(l=1,\dots L)$ other probabilistic observations. The extension of these ideas to the PLSA is obtained by considering the factorizations

$$\begin{array}{cc}\hfill P(d_i,w_j,x_l)& =\sum _{p}P\left(d_i\right|x_p)P\left(w_j\right|z_r)P\left(z_k\right|z_r)P(x_p,y_q,z_k)(k<r)\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & =\sum _{r}P\left(d_i\right|x_r)P\left(w_j\right|x_r)P\left(z_k\right|x_r)P\left(x_r\right)\hfill \end{array}$$

(23)

Those decompositions are the tensorial cases of the asymmetric and symmetric formulations given by Formulas (2) and (3).

Two methods exist for adjusting Formulas (22) and (23): parafrac [66] (parallel factor analysis), assuming a linear approximation of the fibers (the one-dimensional structures that can be extracted from $P(d_i,w_j,z_k)$) and Tucker [67], a multiway Principal component analysis (PCA). Peng has noted that both methods provide different results even when the objective function is the same [65], and has indicated that the method is useful to determine the number of latent factors. An alternative formulation has been proposed by [68].

This decomposition has several implications, mainly related to neural network applications. For details and examples of applications, we refer readers to [69, Chap. 7].

5. NMF Point of View

The algebraic object that supports the probabilities of (2) or (3) are matrices with restricted entries to the set $[0,1]$, and they are non-negative. To construct such matrices, the transformation (1) involves identifying $N(d_i,w_j)$ to a multivariate matrix $\mathbf{X}$. The matrix $\mathbf{Y}$ containing the probabilities $P(d_i,w_j)$ is obtained with the transformation $\mathbf{Y}=\mathbf{X}/{\sum }_{ij}\mathbf{X}$. This is a special case of probabilistic transformation in which probabilities are in Laplace’s sense or relative frequencies.

For matrices obtained in this way, the standard formulation of the NMF is [69, p. 131]

$$\begin{array}{cc}\hfill {\left[\mathbf{Y}\right]}_{ij}& ={\left[\mathbf{W}\right]}_{ik}{\left[\mathbf{H}\right]}_{kj}+{\left[\mathbf{E}\right]}_{ij}(k=1,2,\dots )\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill & \approx {\left[\mathbf{W}\right]}_{ik}{\left[\mathbf{H}\right]}_{kj}\hfill \end{array}$$

(25)

where $\mathbf{E}$ is the error matrix, and makes the NMF suitable for its use in alternative formulations of the PLSA [70]1.

Many authors attribute the introduction of this technique to the work of Paatero [74], while others do so with the publication of Lee and Seung [75]. Both approaches are not equivalent. While Paatero uses the Euclidean norm as the objective function, Lee and Seung use the I-divergence (distances d are maps that satisfy, for vectors $\mathbf{a}$, $\mathbf{b}$, and $\mathbf{c}$, the following axioms: $\left(i\right)$ symmetry, $d(\mathbf{a},\mathbf{b})=d(\mathbf{b},\mathbf{a})$; $\left(ii\right)$ identity $d(\mathbf{a},\mathbf{b})=0$ if $\mathbf{a}=\mathbf{b}$; and $\left(iii\right)$ (triangular inequality) $d(\mathbf{a},\mathbf{b})\le d(\mathbf{a},\mathbf{c})+d(\mathbf{c},\mathbf{b})$. A divergence D does not satisfy one of these axioms, usually symmetry, which is more suitable for measuring how densities are similar.). Furthermore, the work of Lee and Seung is focused on the clustering problem. This attribution creates conceptual errors in many works, identifying NMF techniques to classification. A previous and algebraically rigorous and sound formulation of the NMF is a debt of Chen [76]. A brief introduction to NMF as an optimization problem can be found in [77, Chap. 6]; a more standard introduction is provided in [48, Chap. 7].

On the other hand, the SVD arose from the efforts of several generations of mathematicians, from the nineteenth-century works of Beltrami [78], and independently Jordan [79], to more recent contributions regarding inequalities between eigenvalues and matrix norms by Ky-Fan [80,81]. Currently, the SVD plays a central role in algebra, constituting a field known as eigenanalysis, and serves as a basis for matrix function theory [82], and is also the basis of many multivariate methods. This research field remains active. Currently it is formulated as [77, p. 275]

Theorem 1. 

Let $\mathbf{X}\in {\Re }^{m\times n}$ (or ${\mathbb{C}}^{m\times n}$); then orthogonal (or unitary) matrices $\mathbf{U}\in {\Re }^{m\times m}$ (or $\mathbf{U}\in {\mathbb{C}}^{m\times m}$) and $\mathbf{V}\in {\Re }^{n\times n}$ (or $\mathbf{V}\in {\mathbb{C}}^{m\times m}$) exist such that

$$\begin{array}{c}\hfill \mathbf{X}=\mathbf{U}\Sigma {\mathbf{V}}^t(or\mathbf{X}=\mathbf{U}\Sigma {\mathbf{V}}^H)\Sigma =\left[\begin{array}{cc}{\Sigma }_1& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right]\end{array}$$

(26)

where ${\Sigma }_1=diag({\sigma }_1,\dots ,{\sigma }_r)$ with diagonal entries

$$\begin{array}{c}\hfill {\sigma }_1\ge ,\dots ,\ge {\sigma }_r>0r=rank\left(\mathbf{A}\right)\end{array}$$

One of the first proofs can be found in [83]. The theorem as given is known as full rank SVD. The approximation for $r^{\prime }<r$ is known as low-rank approximation, assuming an approximation for (26) [84]. In the PLSA context, connected with probabilities, only real matrices are used.

Hofmann has related PLSA (in the symmetric formulation case) to SVD in conference proceedings [1] and [3], writing the formal equivalence

$${\left[\mathbf{U}\right]}_{ik}\sim P\left(d_i\right|z_k)$$

(27a)

$$\mathrm{diag}{\left(\Sigma \right)}_k\sim P\left(z_k\right)$$

(27b)

$${\left[\mathbf{V}\right]}_{kj}\sim P\left(w_j\right|z_k)$$

(27c)

where $\mathbf{U}$, $\mathrm{diag}\left(\Sigma \right)$, and $\mathbf{V}$ are related to the SVD of the matrix $\mathbf{Y}$.

The relationships between the PLSA and the SVD and these relations have severe restrictions because the data are frequencies obtained from counts obeying multinomial laws, whereas SVD exists for every matrix of real entries. In addition, the conditions for the degree of adjustment of $n(d_i,w_j)$ to $\mathbf{Y}$ are unclear since is not defined the approximation bound. Also, the popssibility of the use of smooth tecniques for cases that data are not frequencies, is omitted. The relations (27a)-(27c), first written by Hofmann, was considered a mere formal equivalence [3,38].

Several attempts focusing on the equivalence between PLSA and SVD, in light of NMF, have aimed to build a more rigorous relation. The explicit relationship between PLSA and NMF, stated by Gaussier, minimizes I-divergence 2 with non-negative constraints

$$\begin{array}{cc}\hfill & {[·]}_{ij}\ge 0\hfill \\ \hfill & \nabla D(·)\ge 0\hfill \\ \hfill & {[·]}_{ij}\odot \nabla D_I=0\hfill \end{array}$$

known as Karush-Kuhn-Tucker (KKT) conditions, where ⊙ is the Hadamard or element-wise product. KKT conditions are a widespread optimization method when divergences are used.

Solutions are [75]

$$\begin{array}{c}\hfill {\left[\mathbf{W}\right]}_{ik}\leftarrow {\left[\mathbf{W}\right]}_{ik}\odot \frac{{\left[\mathbf{Y}{\mathbf{H}}^t\right]}_{ik}}{{\displaystyle {\left[\mathbf{WH}{\mathbf{H}}^t\right]}_{ik}}}\end{array}$$

(28)

$$\begin{array}{c}\hfill {\left[\mathbf{H}\right]}_{kj}\leftarrow {\left[\mathbf{H}\right]}_{kj}\odot \frac{{\left[{\mathbf{W}}^t\mathbf{Y}\right]}_{kj}}{{\displaystyle {\left[{\mathbf{W}}^t\mathbf{WH}\right]}_{kj}}}\end{array}$$

(29)

and the matrix quotient is the element-wise entry division.

After adjusting equation (25) in an iterative process, consisting of selecting a value of k, switching between (28) and (29) until a satisfactory approximation degree is achieved, Gaussier has introduced diagonal matrices ${\mathbf{D}}_1$ and ${\mathbf{D}}_2$ of suitable dimension

$$\begin{array}{cc}\hfill {\left[\mathbf{W}\mathbf{H}\right]}_{ij}& ={\left[\left({\mathbf{WD}}_1^{-1}{\mathbf{D}}_1\right)\right]}_{ik}{\left[\left({\mathbf{D}}_2{\mathbf{D}}_2^{-1}\mathbf{H}\right)\right]}_{kj}\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill & ={\left[\left(\mathbf{W}{\mathbf{D}}_1^{-1}\right)\right]}_{ik}\mathrm{diag}\left[\left({\mathbf{D}}_1{\mathbf{D}}_2\right)\right]{\left[\left({\mathbf{D}}_2^{-1}\mathbf{H}\right)\right]}_{kj}\hfill \end{array}$$

(31)

stating that any (local) maximum solution of PLSA is a solution of the NMF with KL-divergence (I-divergence according to the nomenclature herein) [11].

Further work by Ding [86], with the same divergence, has introduced normalization for matrices $\mathbf{W}$ and $\mathbf{H}$, such that the column stochastic matrix $\tilde{\mathbf{W}}=[{\tilde{\mathbf{w}}}_1,\cdots ,{\tilde{\mathbf{w}}}_K]$ and the row stochastic matrix $\tilde{\mathbf{H}}=[{\tilde{\mathbf{h}}}_1,\cdots ,{\tilde{\mathbf{h}}}_K]$ are obtained as

$$\begin{array}{cc}\hfill {\tilde{\mathbf{w}}}_k& =\frac{{\mathbf{w}}_k}{{\sum }_i{w}_{ik}}=1\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\tilde{\mathbf{h}}}_k& =\frac{{\mathbf{h}}_k}{{\sum }_j{h}_{kj}}=1\hfill \end{array}$$

(33)

calling those conditions probabilistic normalizations, and writing

$$\begin{array}{cc}\hfill \mathbf{Y}& =\tilde{\mathbf{W}}{\mathbf{D}}_W\tilde{\mathbf{H}}{\mathbf{D}}_H\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill & =\tilde{\mathbf{W}}\mathbf{S}\tilde{\mathbf{H}}(\mathrm{s}.\mathrm{t}.\mathbf{S}={\mathbf{D}}_W{\mathbf{D}}_H)\hfill \end{array}$$

(35)

where the ${\mathbf{D}}_W$ and ${\mathbf{D}}_H$ diagonal matrices contain the column sums of the respective sub-index matrices. Ding has arrived at similar conclusions to Gaussier, and assimilated the latent variables into the space span of matrix factorization [87].

Conditions for the reverse result are shown in [18] by the KL divergence

$$\begin{array}{cc}\hfill D_{KL}(\mathbf{Y}\parallel \mathbf{W}\mathbf{H})& =\sum _{ij}{\left[\mathbf{Y}\right]}_{ij}log\frac{{\left[\mathbf{Y}\right]}_{ij}}{{\displaystyle {\left[\mathbf{WH}\right]}_{ij}}}\hfill \end{array}$$

(36)

obtaining the solutions

$$\begin{array}{cc}\hfill {\left[\mathbf{W}\right]}_{ik}& \leftarrow {\left[\mathbf{W}\right]}_{ik}\odot \left(\frac{{\left[\mathbf{Y}\right]}_{ij}}{{\left[\mathbf{W}\mathbf{H}\right]}_{ij}}{\left[\mathbf{H}\right]}_{kj}^t\right)\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill {\left[\mathbf{H}\right]}_{kj}& \leftarrow {\left[\mathbf{H}\right]}_{kj}\odot \left({\left[\mathbf{W}\right]}_{ik}^t\frac{{\left[\mathbf{Y}\right]}_{ij}}{{\displaystyle {\left[\mathbf{W}\mathbf{H}\right]}_{ij}}}\right)\hfill \end{array}$$

(38)

after proof that $\mathbf{WH}\to \mathbf{Y}$ if $k\ge min(m,n)$, choosing the diagonal matrix as

$$\begin{array}{c}\hfill \mathbf{t}=\frac{\mathrm{diag}{\left({\left[\mathbf{WH}\right]}_{ij}^t{\left[\mathbf{WH}\right]}_{ij}\right)}^{1/2}}{\mathrm{trace}{\left({\left[\mathbf{WH}\right]}_{ij}^t{\left[\mathbf{WH}\right]}_{ij}\right)}^{1/2}}\end{array}$$

(39)

and arranging the entries of $\mathbf{t}$ in decreasing order, with the same permutation on the columns of $\mathbf{W}$ and the rows of $\mathbf{H}$, and obtaining the respective column and row stochastic matrices $\tilde{\mathbf{F}}$ and $\tilde{\mathbf{G}}$, indicating that

$$\begin{array}{c}\hfill {\left[\mathbf{WH}\right]}_{ij}={\left[\tilde{\mathbf{F}}\right]}_{ik}\mathrm{diag}\left(\mathbf{t}\right){\left[\tilde{\mathbf{G}}\right]}_{kj}\end{array}$$

(40)

In this case, factorization 37 reaches the SVD of the orthonormalization of $\mathbf{Y}$ (see [88, p. 24] for orthonormalization process).

This procedure keeps matrix norms (also row or columns norms) [63]. Moreover, minimization of KL divergence is equivalent to maximization of the likelihood in certain cases (as can easily be seen by expanding the logarithm of the KL divergence as a difference; while the first term is a constant, the second term is the log-likelihood), but this is not exact. The minimization of the $KL$ divergence is known as the em algorithm. In many cases, the results obtained with both methods are similar. Amari has shown that in the general case, the em solutions are asymptotes of the EM algorithm case [89].

6. Extensions

The possibility to formulate the PLSA from the NMF was early noticed by [3]. For the symmetric case, modifying the hypotheses on the nature of the data can constitute the basis for other techniques, furnishing probabilistic sense. In this section, we present them in the chronological order of appearance.

6.1. Kernelization

The dot product is used to measure similarity among instances. The transformation of the scalar products of the observations to a different space (not necessarily of the same dimension) is called kernelization and in fact is a generalization of the dot product, transforming $\langle x_{i_1},x_{i_2}\rangle$ to $K(x_{i_1},x_{i_2})$. The PLSA symmetric formulation allows for building a Fisher kernel. This approach, proposed by [3], despite computational difficulties in supporting messy data, has found practical applications in analysis of document similarity [90].

The Fisher kernel is defined as [91]

$$\begin{array}{c}\hfill K(\mathbf{y},{\mathbf{y}}^t)=U_{\theta }\left(\mathbf{y}\right){\mathcal{I}}_F^{-1}{U}_{\theta }\left(\mathbf{y}\right)\end{array}$$

(41)

where

$$\begin{array}{c}\hfill U_{\theta }\left(\mathbf{y}\right)=-\frac{\partial }{\partial \theta }logP\left(\mathbf{y}\right|\theta )\end{array}$$

(42)

the Fisher scores, and

$$\begin{array}{c}\hfill {\mathcal{I}}_F=E_{\mathbf{Y}}\left[U_{\theta }{U}_{\theta }^t\right]\end{array}$$

(43)

the Fisher information matrix. This kernel provides a consistent estimator of the posterior [92]. Hofmann’s proposal in [13] is

$$\begin{array}{cc}\hfill K(d_i,d_i^{\prime })& =\langle u(d_i;\widehat{\theta }){\mathcal{I}}_F{\left(\widehat{\theta }\right)}^{-1};u\left(d_j\widehat{\theta }\right)\rangle \hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill & =\sum _j\widehat{P}(d_i,w_j)\widehat{P}(d_i^{\prime },w_j)\sum _k\frac{P\left(z_k\right|d_i,w_j)P\left(z_k\right|d_i,w_j)}{P(w_j,z_k)}\hfill \end{array}$$

(45)

by direct computation, and $\widehat{P}$ denote the documents in which the distance is measured. A later version is [90], assuming only iid mixtures.

A related technique is graph-regularized PLSA [93]. The objective of this method is to classify entities according topics, according to probabilistic criteria, to measure similarity.

In addition, NMF enables use of the generalization of the dot product to measure similarity and preserve consistency [94].

6.2. Principal Component Analysis

PCA is one of the most extended multivariate and data analysis tools. It can be considered as a particular case of the SVD, being the terms SVD and PCA sometimes interchanged. The objective is to find an orthogonal axis system in Euclidean space maximizing the variance. From a statistical point of view, this representation is a descriptive method. In addition, several attempts have been made to provide a probabilistic sense for PCA [95,96], and establishing a relationship with PLSA seems natural.

From relation (26), and restricted to the case of real matrices

$$\begin{array}{c}\hfill {\Sigma }^2={\left({\left[\mathbf{U}\right]}_{ik}\Sigma {\left[\mathbf{V}\right]}_{ik}\right)}^t\left({\left[\mathbf{U}\right]}_{ik}\Sigma {\left[\mathbf{V}\right]}_{ik}\right)\end{array}$$

(46)

relates ${\Sigma }^2$ of the SVD theorem to the variance matrix $\mathbf{S}$ when ${\mathbf{X}}_c$ is the centered matrix obtained from $\mathbf{X}$ as

$$\begin{array}{c}\hfill {\left[{\mathbf{X}}_c\right]}_{ij}={\left[\mathbf{X}\right]}_{ij}-{\left[\mathbf{J}{\overline{x}}_1|\dots |\mathbf{J}{\overline{x}}_n\right]}_{ij}\left(\mathrm{with}{\overline{x}}_j=\frac{1}{m}\sum _i{x}_{ij}\mathrm{for}\mathrm{all}j\right)\end{array}$$

(47)

where $\mathbf{J}$ is a $m\times 1$ dimension matrix of ones, and the second term relation (47) is the expectation of $\mathbf{X}$. It is immediate

$$\begin{array}{cc}\hfill {\Sigma }^2& ={\left[{\mathbf{X}}_c\right]}_{ij}^t{\left[{\mathbf{X}}_c\right]}_{ij}\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill & =E{\left({\left[\mathbf{X}\right]}_{ij}-E\left({\left[\mathbf{X}\right]}_{ij}\right)\right)}^{t}E\left({\left[\mathbf{X}\right]}_{ij}-E\left({\left[\mathbf{X}\right]}_{ij}\right)\right)\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill & =\mathbf{S}\hfill \end{array}$$

(50)

PCA provides graphical representations for considering the orthogonal projections of observations on the planes formed by the consecutive pairs of columns of the matrix $\mathbf{V}$, which are orthogonal as a consequence of the SVD (i.e., ${\Sigma }^{1/2}\mathbf{V}$).

The relationship between the planes in which the PCA and the PLSA project entities was not immediately clear but has been determined in light of the NMF. However, the column vectors of the matrices of Formula (40) are not necessarily orthogonal but are non-negative. Interpreting probabilities as coordinates, Klingenberg has introduced simplicial cones $\mathsf{\Gamma }$ [15]

$$\begin{array}{c}\hfill \mathsf{\Gamma }=\left\{{\mathbf{y}}_j\mathrm{s}.\mathrm{t}.\mathbf{y}=\sum _j{\alpha }_j{\mathbf{h}}_j\mathrm{with}{\alpha }_j\ge 0\mathrm{and}{\mathbf{h}}_j\in {\left[\mathbf{H}\right]}_{kj}\right\}\end{array}$$

(51)

which is a convex region in the positive orthant. Figure 2 illustrates this transformation.

A formulation known as logistic PCA, described in [97], formulates the likelihood optimization problem of $\mathcal{L}\left(\mathbf{WH}\right)$

$$\begin{array}{cc}\hfill \mathcal{L}& =P\left(\mathbf{Y}\right|\mathbf{WH})\hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill & {=\sigma (\left[\mathbf{WH}\right]}_{ij}^{{\sum }_{ij}{\left[\mathbf{Y}\right]}_{ij}}(1-\sigma {\left({\left[\mathbf{WH}\right]}_{ij}\right)}^{1-{\sum }_{ij}{\left[\mathbf{Y}\right]}_{ij}}\hfill \end{array}$$

(53)

being $\sigma \left(\mathbf{WH}\right)={(1+exp\left(\mathbf{WH}\right))}^{-1}$. Optimizing the likelihood as a Bernoulli pdf with parameters in $[0,1]$, leads to a model for dichotomous variables.

A comparison between NMF and the PCA has been provided by [98], who have noted that the PCA is capable of finding the global minimum, whereas NMF (interpreting the PCA as a dimension reduction problem and not in the full rank case) does not. In addition, the ranking of factors in the NMF is not ordered, and all are equally important. Moreover, non-negative constraints are violated by PCA.

6.3. Clustering

The relationships between the PLSA and clustering techniques, which have been satisfactorily studied, represent the classification ability of the PLSA. PLSA in fact is a probabilistic clustering method, when latent variables are identified to clusters, as has been done in several studies [12,99].

Probabilistic clustering implies that all entities belong to each cluster with different probabilities (including zero) [100], an idea shared with fuzzy clustering methods [101]. However, in the current state of the art, gaps remain concerning overlapping.

In addition, PLSA can be used for partitional clustering, relating PLSA and k-means. This process involves introducing a Bayesian classifier in the matrix $\mathbf{W}$ of Formula (25) [14], after proof, in the conference paper [12], the connection between NMF and PLSA, and relaxing the assumptions of non-negativity assumptions on the basis matrix [14]. Using this technique, Ding has obtained graphical representations close to the centroids of the k-means [86]. In addition, these ideas have been used to build a simplex model based on topics containing normalized data points [102].

6.4. Information Theory Interpretation

The link between PLSA and information theory is apparent when divergences are used to evaluate the similarity between distributions. It is convenient to recall the introduction of distances or divergences induces metrics when the Cauchy-Schwartz inequality ($\parallel \mathbf{AB}\parallel \le \parallel \mathbf{A}\parallel \parallel \mathbf{B}\parallel$) is satisfied [84]. Hofmann has noted that Euclidean distance implies Gaussian distributions [13]. Although this topic is complicated and beyond the scope of this article, it notably has several implications in the symmetric PLSA interpretation.

The frequentist framework does not provides reliable estimations in some cases, as noted by Rao in population diversity studies [103,104]. Divergences satisfying the identity axiom are related to entropy after the introduction of a monotonically decreasing function J of the differences (or quotient) with parameters $\theta$ and $\varphi$ of the same class of densities (Rao has used Jensen’s difference in [105], defined as $J(\theta ,\varphi )=H(\theta ,\varphi )-\lambda H\left(\theta \right)-\mu H\left(\varphi \right)$, where H is an entropy, and $\lambda$ y $\mu$ scalars such that $\lambda +\mu =1$). On the basis of the assumption that the parameter space is a sufficiently differentiable manifold, the (dis)similarity between populations can be locally approximated by $\varphi \approx \theta +d\theta$, known as the geodesic distance

$$\begin{array}{cc}\hfill g_{ij}& =\sum _{ij}\frac{{\partial }^2}{\partial {\theta }_i\partial {\theta }_j}J(\theta ,\theta +d\theta )\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill & =\sum _{ij}\mathcal{H}\left(\mathrm{s}.\mathrm{t}.\mathcal{H}\in {\mathbb{Re}}^{p\times p}\right)\hfill \end{array}$$

(55)

where $\mathcal{H}$ is the Hessian.

Considering the development for $\varphi =\theta$

$$\begin{array}{c}\hfill J(\theta ,\theta +d\theta )=J(\theta ,\varphi =\theta )+\frac{\partial }{\partial \theta }J(\theta ,\varphi =\theta )+\frac{1}{2!}\frac{{\partial }^2}{\partial {\theta }_i\partial {\theta }_j}J(\theta ,\varphi =\theta )+\cdots \end{array}$$

(56)

the first two terms vanish, and the expectation of $\mathcal{J}$ is the Fisher information matrix (or the inverse variance matrix). The connection with $\Sigma$ of (26) and/or (46) is

$$\begin{array}{cc}\hfill E\left(\mathcal{H}\right)& ={\mathcal{I}}_F\hfill \\ \hfill & =\mathrm{diag}\left(1/{\sigma }_1^2,\dots ,1/{\sigma }_K^2\right)\hfill \end{array}$$

(57)

and as consequence of the Jensen inequality, the bound ${\mathcal{I}}_F\ge 1/\mathbf{S}$ appears.

The general treatment for connecting the divergences and underlying distributions is provided in [105]. This article reproduces the relationships between metrics and distributions obtained by [106]. A more recent treatment based on the concept of kernelization is [107].

The conference paper [16] explicitly relates the PLSA, when the KL divergence is used, to Shannon’s information, as a result of expanding the logarithm of KL divergence

$$\begin{array}{c}\hfill D_{KL}(\mathbf{Y}\parallel \mathbf{W}\mathbf{H})=\sum _{ij}{\left[\mathbf{Y}\right]}_{ij}log{\left[\mathbf{Y}\right]}_{ij}-{\left[\mathbf{Y}\right]}_{ij}log{\left[\mathbf{WH}\right]}_{ij}\end{array}$$

(58)

and identifying terms

$$\begin{array}{c}\hfill I\left(\mathbf{Y}\right|\mathbf{WH})=H(\mathbf{Y})-H(\mathbf{W}\left|\mathbf{H}\right)\end{array}$$

(59)

where I is the mutual information. In this context, there are $r!$ representations (if the entries are labeled) that correspond to the indistinguishable entities (different entities that have the same values for all observational variables). A geometric interpretation of the information appears when the equivalence of the maximization of the likelihood is considered with the EM algorithm and the KL divergence. These results provide a stronger foundation for the probability space projection than for the orthogonal projection [108], as Hofmann has noted, where the divergence is the loss of information [109, p. 185].

Although Chaudiri has used this result for k-means error classification purposes, obtaining a bound for the variance expectation, the general consequence of relating divergences is a parametric estimation of the variance.

6.5. Transfer Learning

Transfer learning can be defined as the machine learning problem of trying to transfer knowledge from a source domain to a target domain [110].

PLSA can be used from the point of view of neural networks for transfer learning purposes, by solving the problem in the case in which the source domain shares only a subset of its classes (column vectors of the data matrix) for an unlabeled target data domain [17]. The log-likelihood expression is thus [17]

$$\begin{array}{cc}\hfill \mathcal{L}=& \sum _{ij}n(d^S,w_j)log\sum _{ij}P\left(d^S\right|z_k^S)P\left(w_j\right|z_k^S)P\left(z_k\right|z_k^S)\hfill \\ \hfill & +\sum _{ij}n(d^T,w_j)log\sum _{ij}P\left(d^T\right|z_k^T)P\left(w_j\right|z_k^T)P\left(z_k\right|z_k^T)\hfill \end{array}$$

(60)

where S indicates that a document is in the source, and T indicates the target domain. A detailed survey providing an introduction to neural networks is [111]. A similar work on the problem of transfer learning is [112].

6.6. Open Questions

The work of Gaussier and Ding has been important in relating the PLSA and the SVD; translating the conceptual framework to the context of the probabilistic interpretations of the NMF; and extending the data class domain from the non-negative integers to the non-negative reals, relaxing distributional assumptions, since no-hypothesis is done on the parameter space. This is a non-parametric method based on NMF algebra. In addition, NMF techniques mainly focus on symmetric formulation. It does not seem to be any objection preventing its use for the asymmetric one, although problems in assigning probabilities $P\left(d_i\right)$ in equation (2) could complicate the problem.

In addition, Ding has stated that the difference between the results of the SVD and NMF (and thus PLSA) depends on the convergence to different local optima, and it is true if $k<min(m,n)$. Ding’s work has considered SVD as a dimensional reduction problem or low-rank decomposition. In this case, the matrix $\mathbf{WH}$ will not be $\mathbf{Y}$, but an approximation, and the SVD is not achievable with NMF. In addition, because PCA decomposition is related to the geometric multiplicities of eigenvalues, in the case of the relations obtained with NMF, it is not so clear, and it should be faced with algebraic dimensionality 3.

However, when the discussion is restricted to the symmetric formulation, questions arise, and the results depend on k and determining the equivalence of the solutions. Leaving aside the type of convergence of $\mathbf{WH}\to \mathbf{Y}$, which holds in the case $k\ge min(m,n)$, as indicated by several authors [15,18,114], suboptimality occurs when this condition is not fulfilled and implies that the SVD low-rank approximation is an ill-conditioned problem [15].

On the other hand, the similitude of PLSA and SVD presents a greater issue. Whereas SVD exists for every matrix, NMF requires the matrix entries to be restricted to non-negative values. Thus, the problem can be stated as the triplet $(\mathsf{\Omega },\mathcal{A},\mathcal{P})$ (as usual, $\mathsf{\Omega }$ is the sample space, $\mathcal{A}$ is a $\sigma$-algebra, and $\mathcal{P}$ is the measure or probability) . In the current state of the art, $\mathsf{\Omega }$ is restricted to ${\mathbb{R}}_{+}^{m\times n}$. A treatment to expand the sample space to real values could involve well-known transformations (as suggested by $\phi \left(\mathbf{X}\right)=1/(1+exp\{\mathbf{X}\left\}\right)$) and may be explored in the future.

7. PLSA Processing Steps and State of the Art of Solutions

PLSA is considered an effective technique but has a notable drawback in its high consumption of computing resources, in terms of both execution and internal memory. This drawback has limited its practical applications [115] and additionally makes the relationship between the SVD and PLSA curious. In the SVD case, the typical blackboard exercise of obtaining eigenvalues and eigenvectors is simple but does not occur in the same manner for moderate and large data sets. Methods for its effective computation have arisen from numerous studies and sustained efforts over several decades [116]. Currently, many program languages implement the Linear Algebra Package to facilitate SVD computation [117]. Also, solutions for PLSA are hard to obtain.

Beyond the EM algorithm problems, PLSA is highly dependent on the initialization values [118,119]. This leads to several algorithms for computational efficiency purposes, based on certain initialization conditions, and others on alternative versions of the EM algorithm, apart from those that strictly use computational techniques.

Herein, contributions to increasing the computational efficiency are examined according to the concepts on which they are based, their initialization conditions and the use of EM algorithm variants. Efforts using purely computational techniques are also discussed.

7.1. Algorithm Initialization

The dependence of the PLSA results on the initialization conditions has led to several variations. One possibility applicable only in the symmetric formulation, proposed by [118], initializes the algorithm with LSA solutions, which are the SVD solutions. Because some values can be negative, correction may be necessary (typically setting values to zero). Another strategy applicable in both formulations is execution for several random initialization distributions of the considered algorithm; after running, the higher log-likelihood value offers the best solution [119].

One algorithm is On-line belief propagation (OBP), which is based on a sequence of initializations on subsets of the data frame [119,120]. OBP segments the data frame into several parts. After initialization of the first segmentation, solutions are obtained and used in the next initialization, and so on. This technique enables the use of PLSA on large data sets.

A fundamental of the OBP is stochastic initialization [121], which consists of defining a learning function as a risk function for which the difference in conditional distributions describes a decreasing sequence between iterations [121]. The execution of this algorithm requires at least one iteration for the complete data set and selection of the most significant contributions for the first partition.

7.2. Algorithms Based on Expectation Maximization Improvement

The EM convergence rate is [122]

$$\begin{array}{c}\hfill \parallel {\theta }^{(p+1)}-{\theta }^{★}\parallel \le \lambda \parallel {\theta }^{\left(p\right)}-{\theta }^{★}\parallel \end{array}$$

(61)

where $\lambda$ is the largest eigenvalue of the data matrix. Several methods are used to accelerate convergence, such as the descendant gradient. However, PLSA must preserve maximum log-likelihood solutions. To improve computational efficiency in such conditions, some variants and alternative algorithms have been proposed. The EM algorithm is one of the most studied in statistical environments, and many variants and simplifications exist [123]. A general description of the algorithm used in this section is [124]. An overview of the efforts is provided in Table 2. The EM algorithm is the classic optimization technique for PLSA, and some versions or modifications have been exploited to achieve PLSA solutions.

7.2.1. Tempered EM

Tempered EM uses classical concepts of statistical mechanics for computational purposes [125]. Aside from the significance in physics, the primary idea is how to achieve a posterior (E-step) close to a uniform distribution. An objective function is introduced

$$\begin{array}{cc}\hfill {\mathcal{F}}_{\beta }& =-\beta \sum _{ij}n(d_i,w_j)\sum _k\tilde{P}(z_k;d_i,w_j)log\left[P\left(d_i\right|z_k)P\left(w_j\right|z_k)P\left(z_k\right)\right]\hfill \\ \hfill & +\sum _{ij}n(d_i,w_j)\sum _k\tilde{P}(z_k;d_i,w_j)logP(z_k;d_i,w_j)\hfill \end{array}$$

(62)

where $\tilde{P}(z_k;d_i,w_j)$ is a variational parameter defined as

$$\begin{array}{c}\hfill \tilde{P}(z_k;d_i,w_j)=\frac{{\left[P\left(z_k\right)P\left(d_i\right|z_k)P\left(w_j\right|z_k)\right]}^{\beta }}{{\sum }_k{\left[P\left(z_k\right)P\left(d_i\right|z_k)P\left(w_j\right|z_k)\right]}^{\beta }}\end{array}$$

(63)

and for $\beta <1$, the convergence is faster [13].

7.2.2. Sparse PLSA

A proposal to improve the convergence speed has been based on sparse EM [126]. Assuming that only a subset of values is plausible for latent variables (in terms of probabilities), freezing non-significant avoids many calculations. PLSA is considered as an algebraic optimization problem of the matrix $\mathbf{Y}\in {\mathbb{R}}^{m\times n}$ (which in this case is the data frame containing the relative frequencies $n(d_i,w_j)$) restricted to the constraint ${\sum }_r{\lambda }_r{\mathbf{y}}_r{\mathbf{y}}_r^t(r<m)$, or unknown parameters, minimizing [127]

$$\begin{array}{c}\hfill D_q(\mathbf{Y}\parallel \sum _r{\lambda }_r{\mathbf{y}}_r{\mathbf{y}}_{r^{\prime }})(\mathrm{with}\sum _r{\lambda }_r=\parallel {\mathbf{y}}_r{\parallel }_1=\parallel {\mathbf{y}}_r^{\prime }{\parallel }_1=1,{\mathbf{y}}_r\in \mathbf{Y}\mathrm{and}{r}^{\prime }\ne r<n)\end{array}$$

(64)

named Tsallis divergence [128], and computed for the r non-freezing column vectors of $\mathbf{Y}$ as [69, p.97]

$$\begin{array}{c}\hfill D_q({\mathbf{y}}_j\parallel {\lambda }_r{\mathbf{y}}_r{\mathbf{y}}_r^t)=\frac{1}{\kappa }\sum _i\left({\mathbf{y}}_j\left({\mathbf{y}}_j^{\kappa }-{\left({\lambda }_r{\mathbf{y}}_r{\mathbf{y}}_r^t\right)}^{\kappa }\right)\right)-\sum _i\left({\mathbf{y}}_j^{\kappa }\left({\mathbf{y}}_j-{\lambda }_r{\mathbf{y}}_r{\mathbf{y}}_r^t\right)\right)(\mathrm{s}.\mathrm{t}.\kappa \ne 0)\end{array}$$

(65)

This divergence solves the optimization problem of adjusting $n(d_i,w_j)$ to $P(d_i,w_j)$ [129]. After adjustment, probabilistic factorizations of the considered parametrization must again be obtained.

7.2.3. Incremental PLSA

Instead of global maximization, simpler contributions can be maximized. This update procedure used in the E-step for the PLSA gave rise to the incremental PLSA algorithm [130], with which results can be obtained twice as quickly. Applications in image classification can be found in [131,132].

A recursive algorithm, called recursive probabilistic latent semantic analysis, is based on the computation of the likelihood of a subset of words, as well as other words, recursively [133]. The performance has been reported to be highly similar to that obtained with the incremental PLSA.

7.3. Use of Computational Techniques

Difficulties in obtaining fast and reliable solutions for PLSA have also been approached through purely computational techniques. These advancements are a consequence of developments in computer architecture in recent decades: processing capabilities have been increased, thus resulting in a new branch of algorithms to reduce the computational time of the PLSA. The introduction of multicore processors by Intel and Sun Microsystems, in 2005, for portable machines enabled a major step toward parallel computing [134], which is now the dominant paradigm.

Parallel computing involves simultaneous execution of tasks. It requires dividing a problem into independent pieces and executing each one in a separate processing unit. The use of parallel computing techniques for the PLSA has been proposed in [135]. A current and widespread technique to support parallel capabilities is Map Reduce [136]. This technique essentially consists of dividing tasks into two phases. The first phase is a map that partitions the input data set and assigns labels to each one. The reduce phase supposes the execution of an operation on a set of previously labeled partitions. An algorithm exploiting the possibilities of Map Reduce for PLSA results has been proposed in [137].

Furthermore, graphic processing units have increased the range of capabilities, and they are useful for a broad variety of applications, particularly the simulation of complex models [138]. These capabilities have been transferred to the PLSA algorithm [139] but have not yielded definitive results.

7.4. Open Questions

The described methods do not provide fully satisfactory results, and perhaps it is one of the causes of the division between enthusiasts and indifferent regarding PLSA. The PLSA algorithms inherit the problems of EM, especially the slow convergence. Those problems are independent of the computational efficiency. There are no comparative studies on computational capability or execution time. Surprisingly, despite the EM algorithm is one of the most studied, many versions and works devoted to accelerating convergence (a recent one is [140]), there are no comparative studies.

8. Discussion

Hofmann does not provide indications as to when each formulation is applicable, the value of PLSA is clear and relatively simple: for the asymmetric formulation, it is an unsupervised learning method to train a set of latent variables with unknown number and distribution, when the data include co-occurrence or contingency tables. Also, the model can be extended to continuously evaluated entities. The use of the EM algorithm to adjust probabilities provides a maximum likelihood estimation of the parameters. This technique can be extended to semi-supervised cases. However, the symmetric formulation pales the PLSA within the fundamentals of algebra and multivariate analysis: it is the probabilistic companion of the SVD. However, there are some gaps to make full sense of this statement.

Convergence occurs for $k\ge min(m,n)$, however, for the case $k<min(m,n)$. It has been said that the convergence limit does not necessarily occur at a global optimum [141] and does not necessarily converge to a point but can converge on a compact set [142], thus providing sub-optimal results [143]. In addition, sparse data structures can cause failures in convergence [144]. These statements are due to the interpretation of the SVD, and therefore the PCA as a low-rank approximation. However, these statements do not take into account the low-rank approximation, based on Schmidt’s approximation Theorem [145]. Establishing a bound with the help of this result is an open question.

The LDA techniques and those built on them are hierarchical models whose construction corresponds to particular fields of application. These constructions are also possible in the symmetric formulation assuming additional hypotheses on the data, like distribution, qualitative categories, or hierarchical relations. This type of treatment supposes a pre-processing of the data, preserving the content of the PLSA significance.

Furthermore, although the concept of probabilistic learning is sound, based on Vaillant’s work [146], symmetric PLSA is especially apt in the context of transfer learning. In this way, the certainty depends on the available data, as suggested by the relation (60). Then, reanalyzing a problem with new (or complementary) data or observational variables can provide learning sequences.

Although the applications of the PLSA are numerous, more surprises are likely to be encountered. However, theoretical studies from the perspective of the relationship of PLSA to SVD are scarce. Such studies are necessary for broader interpretation. In particular, data matrices of observed phenomena with real values must be related to probabilities to extend the scope of PLSA applicability.

9. Conclusion

The PLSA is a technique with a quarter of a century of existence. It has been spanned to many research areas with good results. Despite the formal equivalence of the formulations, the asymmetric formulation is an IR technique, while the symmetric formulation also allows for establishing a probabilistic relationship with the SVD. Some consequences of this relationship are that they constitu the foundation for probabilistic construction of other techniques, like kernelization, PCA, clustering, or Transfer Learning, as well as the possibility of building a Fisher kernel. However, there are some open questions, of which we highlight the approximation error when using the low-rank approximation and the poor computational efficiency.