Survey on Probabilistic Models of Low-Rank Matrix Factorizations

Abstract

Low-rank matrix factorizations such as Principal Component Analysis (PCA), Singular Value Decomposition (SVD) and Non-negative Matrix Factorization (NMF) are a large class of methods for pursuing the low-rank approximation of a given data matrix. The conventional factorization models are based on the assumption that the data matrices are contaminated stochastically by some type of noise. Thus the point estimations of low-rank components can be obtained by Maximum Likelihood (ML) estimation or Maximum a posteriori (MAP). In the past decade, a variety of probabilistic models of low-rank matrix factorizations have emerged. The most significant difference between low-rank matrix factorizations and their corresponding probabilistic models is that the latter treat the low-rank components as random variables. This paper makes a survey of the probabilistic models of low-rank matrix factorizations. Firstly, we review some probability distributions commonly-used in probabilistic models of low-rank matrix factorizations and introduce the conjugate priors of some probability distributions to simplify the Bayesian inference. Then we provide two main inference methods for probabilistic low-rank matrix factorizations, i.e., Gibbs sampling and variational Bayesian inference. Next, we classify roughly the important probabilistic models of low-rank matrix factorizations into several categories and review them respectively. The categories are performed via different matrix factorizations formulations, which mainly include PCA, matrix factorizations, robust PCA, NMF and tensor factorizations. Finally, we discuss the research issues needed to be studied in the future.

1. Introduction

In many practical applications, a commonly-used assumption is that the dataset approximately lies in a low-dimensional linear subspace. Low-rank matrix factorizations (or decompositions, the two terms are used interchangeably) are just a type of method for recovering the low-rank structure, removing noise and completing missing values. Principal Component Analysis (PCA) [1], a traditional matrix factorization method, copes effectively with the situation that the dataset is contaminated by Gaussian noise. If the mean vector is set to be a zero vector, then PCA is transformed into Singular Value Decomposition (SVD) [2]. Classical PCA approximates the low-rank representation according to eigen decomposition of the covariance matrix of a dataset to be investigated. When there are outliers or large sparse errors, the original version of PCA does not work well for obtaining the low-rank representations. In this case, many robust methods of PCA have been extensively studied such as L1-PCA [3], L1-norm maximization PCA [4], L21-norm maximization PCA [5] and so on.

In this paper, low-rank matrix factorizations refer to a general formulation of matrix factorizations and they mainly consist of PCA, matrix factorizations, robust PCA [6], Non-negative Matrix Factorization (NMF) [7] and tensor decompositions [8]. As a matter of fact, matrix factorizations are usually a special case of PCA and they directly factorize a data matrix into the product of two low-rank matrices without consideration of the mean vector. Matrix completion aims to complete all missing values according to the approximate low-rank structure and it is mathematically described as a nuclear norm minimization problem [9]. To this end, we regard matrix completion as a special case of matrix factorizations. In addition, robust PCA decomposes the data matrix into the superposition of a low-rank matrix and a sparse component, and recovers both the low-rank and the sparse matrices via principal component pursuit [6].

Low-rank matrix factorizations suppose that the low-rank components are corrupted by certain random noise and the low-rank matrices are deterministic unknown parameters. Maximum Likelihood (ML) estimation and Maximum a posteriori (MAP) are two most popular methodologies used in estimating the low-rank matrices. A prominent advantage of the aforementioned point estimate methods is that they are simple and easy to implement. However, we can not obtain the probability distributions of the low-rank matrices that are pre-requisite in exploring the generative models. Probabilistic models of low-rank matrix factorizations can approximate the true probability distributions of the low-rank matrices and they have attracted a great deal of attention in the past two decades. These models have been widely applied in the fields of signal and image processing, computer vision, machine learning and so on.

Tipping and Bishop [10] originally presented probabilistic PCA in which the latent variables are assumed to be a unit isotropic Gaussian distribution. Subsequently, several other probabilistic models of PCA were proposed successively, such as Bayesian PCA [11], Robust L1 PCA [12], Bayesian robust PCA [13] and so on. As a special case of probabilistic models of PCA, Bayesian matrix factorization [14] placed zero-mean spherical Gaussian priors on two low-rank matrices and it was applied to collaborative filtering. Probabilistic models of matrix factorizations also include probabilistic matrix factorization [15], Bayesian probabilistic matrix factorization [16], robust Bayesian matrix factorization [17], probabilistic robust matrix factorization [18], sparse Bayesian matrix completion [19], Bayesian robust matrix factorization [20] and Bayesian L1-norm low-rank matrix factorizations [21].

As to robust PCA, we take the small dense noise into account. In other words, the data matrix is decomposed into the sum of a low-rank matrix, a sparse noise matrix and a dense Gaussian noise matrix. Hence, the corresponding probabilistic models are robust to outliers and large sparse noise, and they are mainly composed of Bayesian robust PCA [22], variational Bayesian robust PCA [23] and sparse Bayesian robust PCA [19]. As another type of low-rank matrix factorizations, NMF decomposes a non-negative data matrix into the product of two non-negative low-rank matrices. Compared with PCA, NMF is a technique which learns holistic, not parts-based representations [7]. Different probabilistic models of NMF were proposed in [24,25,26,27,28]. Recently, probabilistic models of low-rank matrix factorizations are also extended to the case of tensor decompositions. Tucker decomposition and CANDECOMP/PARAFAC (CP) decomposition are two most important tensor decompositions. By generalizing the subspace approximation, some new low rank tensor decompositions have emerged, such as the hierarchical Tucker (HT) format [29,30], the tensor tree structure [31] and the tensor train (TT) format [32]. Among them, the TT format is a special case of the HT and the tensor tree structure [33]. The probabilistic models of the Tucker were presented in [34,35,36] and that of the CP were developed in [37,38,39,40,41,42,43,44,45,46,47,48].

For probabilistic models of low-rank matrix factorizations, there are three most frequently used statistical approaches for inferring the corresponding probability distributions or parameters, i.e., Expectation Maximization (EM) [49,50,51,52,53,54], Gibbs sampling(or a Gibbs sampler) [54,55,56,57,58,59,60,61,62] and variational Bayesian (VB) inference [54,62,63,64,65,66,67,68,69,70,71]. EM is an iterative algorithm with guaranteed the local convergence for ML estimation and does not require explicit evaluation of the likelihood function [70]. Although it has many advantages over ML, EM tends to be limited in applications because of its serious requirements for the posterior of the hidden variables and can not be used to solve complex Bayesian models [70]. However, VB inference relaxes some limitations of the EM algorithm and ameliorates its shortcomings. As a means of VB, Gibbs sampling is another method used to infer the probability distributions of all parameters and hidden variables.

This paper provides a survey on probabilistic models of low-rank matrix factorizations. Firstly, we review the significant probability distributions commonly-used in probabilistic low-rank matrix factorizations and introduce the conjugate priors that are essential to Bayesian statistics. Then we present two most important methods for inferring the probability distributions, that is, Gibbs sampling and VB inference. Next, we roughly divide the available low-rank matrix factorization models into five categories: PCA, matrix factorizations, robust PCA, NMF and tensor decompositions. For each category, we provide a detailed overview of the corresponding probabilistic models.

A central task for probabilistic low-rank matrix factorizations is to predict the missing or incomplete data. For the sake of concise descriptions, we do not consider the missing entries in all models except the sparse Bayesian model of matrix completion. The remainder of this paper is listed as below. Section 2 introduces the commonly-used probability distributions and the conjugate priors. Section 3 presents two frequently used inferring methods: Gibbs sampling and VB inference. Probabilistic models of PCA and matrix factorizations are reviewed respectively in Section 4 and Section 5. Section 6 and Section 7 survey probabilistic models of robust PCA and NMF, respectively. Section 8 provides other probabilistic models of low-rank matrix factorizations and probabilistic tensor factorizations. The conclusions and future research directions are drawn in Section 9.

Notation: Let $ℝ$ be the set of real numbers and ${ℝ}_{+}$ the set of non-negative real numbers. We denote scalars with italic letters (e.g., $x$), vectors with bold letters (e.g., $x$), matrices with bold capital letters (e.g., $X$) and sets with italic capital letters (e.g., $X$). Given a matrix $X$, its ith row, jth column and (i, j)th element are expressed as $x_{i\cdot }$, $x_{\cdot j}$ and $X_{ij}$, respectively. If $X$ is square, let $\mathrm{Tr}(X)$ and $|X|$ be the trace and the determinant of $X$, respectively.

2. Probability Distributions and Conjugate Priors

This section will introduce some probability distributions commonly adopted in probabilistic models of low-rank matrix factorizations and discuss the conjugate priors for algebraic and intuitive convenience in Bayesian statistics.

2.1. Probability Distributions

We consider some significant probability distributions that will be needed in the following sections. Given a univariate random variable $x$ or a random vector $x$, we denote its probability density/mass function by $p(x)$ or $p(x)$. Let $\mathbb{E}(x)$ be the expectation of $x$ and $\mathrm{Cov}(x,x)$ the covariance matrix of $x$.

Several probability distributions such as Gamma distribution and Beta distribution deal with the gamma function defined by

$$\Gamma (x)={\displaystyle {\int }_0^{+\infty }{u}^{x-1}}\exp (-u)\mathrm{d}u.$$

(1)

We summarize the probability distributions used frequently in probabilistic models of low-rank matrix factorizations, as shown in

Table 1. Commonly-used probability distributions.

Probability Distribution	Notation	Probability Density/Mass Function	Expectation	Variance/Covariance
Bernoulli distribution	$\mathrm{Bern}(x|\mu )$	${\mu }^x{(1-\mu )}^{1-x},x\in \{0,1\}$	$\mu$	$\mu (1-\mu )$
Poisson distribution	$\mathrm{Poiss}(x|\lambda )$	$\frac{{\lambda }^x}{x!}\exp (-\lambda ),x=0,1,\dots$	$\lambda$	$\lambda$
Uniform distribution	$\mathrm{U}(x|a,b)$	$\frac{1}{b-a},x\in (a,b)$	$\frac{a+b}{2}$	$\frac{{\left(b-a\right)}^2}{12}$
Multivariate Gaussian distribution	$\mathcal{N}(x|\mu ,\Sigma )$	$\frac{1}{{(2\mathsf{\pi })}^{d/2}{|\Sigma |}^{1/2}}\exp \left(-\frac{1}{2}{(x-\mu )}^T{\Sigma }^{-1}(x-\mu )\right),$ $\Sigma$ is a $d\times d$ symmetric, positive definite matrix	$\mathsf{\mu }$	$\Sigma$
Exponential distribution	$\mathrm{Exp}(x|\lambda )$	$\lambda \exp \left(-\lambda x\right),x>0$	$\frac{1}{\lambda }$	$\frac{1}{{\lambda }^2}$
Laplace distribution	$\mathrm{Lap}(x|\mu ,\sigma )$	$\frac{1}{2\sigma }\exp \left(-\frac{|x-\mu |}{\sigma }\right)$	$\mu$	${\sigma }^2$
Gamma distribution	$\mathrm{Gam}(x|a,b)$	$\frac{1}{\Gamma (a)}{b}^a{x}^{a-1}\exp (-bx),x>0$	$\frac{a}{b}$	$\frac{a}{b^2}$
Inverse-Gamma distribution	$\mathrm{IGam}(x|a,b)$	$\frac{1}{\Gamma (a)}{b}^a{x}^{-a-1}\exp \left(-\frac{b}{x}\right),x>0$	$\frac{b}{a-1}$ for $a>1$	$\frac{b^2}{{(a-1)}^2(a-2)}$ for $a>2$
Student’s t-distribution	$\mathrm{St}(x|\mu ,\lambda ,\nu )$	$\frac{\Gamma \left((v+1)/2\right)}{\Gamma (v/2)}{\left(\frac{\lambda }{\mathsf{\pi }v}\right)}^{1/2}{\left(1+\frac{\lambda {(x-\mu )}^2}{v}\right)}^{-(v+1)/2}$	$\mu$	$\frac{v}{\lambda (v-2)}$ for $v>2$
Beta distribution	$\mathrm{Beta}(x|a,b)$	$\frac{\Gamma (a+b)}{\Gamma (a)\Gamma (b)}{x}^{a-1}{(1-x)}^{b-1},x\in [0,1]$	$\frac{a}{a+b}$	$\frac{ab}{{(a+b)}^2(a+b+1)}$
Wishart distribution	$\mathcal{W}(\Lambda |W,v)$	$B(W,v){|\Lambda |}^{(\nu -d-1)/2}\exp \left(-\frac{1}{2}\mathrm{Tr}(W^{-1}\Lambda )\right),$ $\Lambda$ is a $d\times d$ symmetric, positive definite matrix	$vW$	$v(W_{ij}^2+W_{ii}{W}_{jj})$ for ${\Lambda }_{ij}$
Inverse Gaussian distribution	$\mathrm{IG}(x|\mu ,\lambda )$	$\sqrt{\frac{\lambda }{2\mathsf{\pi }{x}^3}}\exp \left(-\frac{\lambda {(x-\mu )}^2}{2{\mu }^{2}x}\right), x>0$	$\mu$	$\frac{{\mu }^3}{\lambda }$
Generalized inverse Gaussian distribution	$\mathrm{GIG}(x|p,a,b)$	$\frac{{(a/b)}^{p/2}{x}^{p-1}}{2K_P\left(\sqrt{ab}\right)}\exp \left(-\frac{1}{2}\left(ax+\frac{b}{x}\right)\right), x>0$	$\frac{\sqrt{b}{K}_{p+1}(c)}{\sqrt{a}{K}_p(c)},c=\sqrt{ab}$	$\frac{b{K}_{p+2}(c)}{a{K}_p(c)}-\frac{b{K}_{p+1}^2(c)}{a{K}_p^2(c)}$

. In this table, we list the notation for each probability distribution, the probability density/mass function, the expectation and the variance/covariance respectively. For the Wishart distribution $\mathcal{W}(\Lambda |W,v)$, the term $B(W,v)$ is given by

$${|W|}^{-v/2}{\left(2^{\nu d/2}{\mathsf{\pi }}^{d(d-1)/4}{\displaystyle \prod _{i=1}^d\Gamma \left((v+1-i)/2\right)}\right)}^{-1}$$

(2)

where the positive integer $v\ge d-1$ is named the degrees of freedom. For the generalized inverse Gaussian distribution $\mathrm{GIG}(x|p,a,b)$, $K_p(c)$ is a modified Bessel function of the second kind. For brevity of notation, we stipulate a few simple representations of a random variable or vector. For instance, if $x$ follows a multivariate Gaussian distribution with mean $\mathsf{\mu }$ and covariance matrix $\mathsf{\Sigma }$, we can write this distribution as $x~\mathcal{N}(x|\mathsf{\mu },\mathsf{\Sigma })$, or $x~\mathcal{N}(\mathsf{\mu },\mathsf{\Sigma })$, or $p(x)=\mathcal{N}(x|\mathsf{\mu },\mathsf{\Sigma })$. In addition, some probability distributions can be extended to the multivariate case under identical conditions. The following cites an example: if random variables $x_i~\mathrm{St}({\mu }_i,{\lambda }_i,v_i)$ and $x_1,x_2,\dots ,x_N$ are independent, then we have $x={(x_1,x_2,\dots ,x_N)}^T~\mathrm{St}(\mathsf{\mu },\mathsf{\lambda },v)$ with the probability density function $p(x)={\displaystyle {\prod }_{i=1}^N\mathrm{St}(x_i|{\mu }_i,{\lambda }_i,v_i)}$, where $\mathsf{\mu }={({\mu }_1,{\mu }_2,\dots ,{\mu }_N)}^T,$ $\mathsf{\lambda }={({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_N)}^T$ and $v={(v_1,v_2,\dots ,v_N)}^T$.

There exist close relationships among probability distributions listed in Table 1, as shown partially in Figure 1. Refer to [54,62,66] for more probability distributions. Moreover, some continuous random variables can be reformulated as Gaussian scale mixtures with other distributions from a hierarchical viewpoint. Now, we give two most frequently used examples as below.

Example 1.

For given parameter $\mu$, if the conditional probability distribution of $x$ is $p(x|\lambda )=\mathcal{N}(x|\mu ,{\lambda }^{-1})$ and the prior distribution of $\lambda$ is $p(\lambda )=\mathrm{IGam}(\lambda |1,1)$, then we have $p(x)=\mathrm{Lap}(x|\mu ,\sqrt{2}/2)$.

The probability density function of $x$ is derived as:

$$\begin{array}{l}p(x)={\displaystyle {\int }_0^{+\infty }p(\lambda )p(x|\lambda )d}\lambda \\ ={\displaystyle {\int }_0^{+\infty }\frac{1}{\sqrt{2\mathsf{\pi }}}{\lambda }^{-3/2}\exp \left(-\frac{1}{\lambda }\right)\exp \left(-\frac{1}{2}\lambda {\left(x-\mu \right)}^2\right)d}\lambda \\ \stackrel{\lambda =2t}{=}{\displaystyle {\int }_0^{+\infty }\frac{1}{2\sqrt{\mathsf{\pi }}}{t}^{-3/2}\exp \left(-\frac{1}{2t}\right)\exp \left(-t{\left(x-\mu \right)}^2\right)d}t\\ =ℒ\left\{\frac{1}{2\sqrt{\mathsf{\pi }}}{t}^{-3/2}\exp \left(-\frac{1}{2t}\right)\right\}{\left(x-\mu \right)}^2.\end{array}$$

(3)

Meanwhile, it holds that

$${ℒ}^{-1}\left\{\frac{\sqrt{2}}{2}\exp \left(-\sqrt{2{(x-\mu )}^2}\right)\right\}(t)=\frac{1}{2\sqrt{\mathsf{\pi }}}{t}^{-3/2}\exp \left(-\frac{1}{2t}\right)$$

(4)

where $ℒ\{\cdot \}$ is the Laplace transform and ${ℒ}^{-1}\{\cdot \}$ is the inverse Laplace transform. Hence, we get $x~\mathrm{Lap}(\mu ,\sqrt{2}/2)$.

Example 2.

For given parameters $\mu$ and $\tau$, if the conditional probability distribution of $x$ is $p(x|u)=\mathcal{N}\left(x|\mu ,{(\tau u)}^{-1}\right)$ and the prior distribution of $u$ is $p(u|v)=\mathrm{Gam}(u|v/2,v/2)$, then it holds that $p(x|v)=\mathrm{St}(x|\mu ,\tau ,v)$, where $v$ is a fixed parameter.

The derivation process for $p(x|v)$ is described as follows:

$$\begin{array}{l}p(x|v)={\displaystyle {\int }_0^{+\infty }p(x,u|v)du}={\displaystyle {\int }_0^{+\infty }p(u|v)p(x|u)du}\\ \propto {\displaystyle {\int }_0^{+\infty }{u}^{(v+1)/2-1}\exp \left(-\frac{u\left(v+{(x-\mu )}^2\tau \right)}{2}\right)du}\\ \propto \frac{1}{{\left(v+{(x-\mu )}^2\tau \right)}^{\left(v+1\right)/2}}.\end{array}$$

(5)

So, $p(x|v)=\mathrm{St}(x|\mu ,\tau ,v)$.

2.2. Conjugate Priors

Let $x$ be a random vector with the parameter vector $z$ and $X=\left\{x_1,x_2,\dots ,x_N\right\}$ a collection of N observed samples. In the presence of latent variables, they are also absorbed into $z$. For given $z$, the conditional probability density/mass function of $x$ is denoted by $p(x|z)$. Thus, we can construct the likelihood function:

$$L(z|X)=p(X|z)={\displaystyle \prod _{i=1}^{N}p(x_i|z)}.$$

(6)

As for variational Bayesian methods, the parameter vector $z$ is usually assumed to be stochastic. Here, the prior distribution of $z$ is expressed as $p(z)$.

To simplify Bayesian analysis, we hope that the posterior distribution $p(z|X)$ is in the same functional form as the prior $p(z)$. Under this circumstance, the prior and the posterior are called conjugate distributions and the prior is also called a conjugate prior for the likelihood function $L(z|X)$ [54,66]. In the following, we provide three most commonly-used examples of conjugate priors.

Example 3.

Assume that random variable $x$ obeys the Bernoulli distribution with parameter $\mu$. We have the likelihood function for $x$:

$$L(\mu |X)={\displaystyle \prod _{i=1}^N\mathrm{Bern}(x_i|\mu )}={\displaystyle \prod _{i=1}^N{\mu }^{x_i}{(1-\mu )}^{1-x_i}}={\mu }^{{\displaystyle {\sum }_{i=1}^N{x}_i}}{(1-\mu )}^{N-{\displaystyle {\sum }_{i=1}^N{x}_i}}$$

(7)

where the observations $x_i\in \{0,1\}$. In consideration of the form of $L(\mu |X)$, we stipulate the prior distribution of $\mu$ as the Beta distribution with parameters $a$ and $b$:

$$p(\mu )=\mathrm{Beta}(\mu |a,b)=\frac{\Gamma (a+b)}{\Gamma (a)\Gamma (b)}{\mu }^{a-1}{(1-\mu )}^{b-1}.$$

(8)

At this moment, we get the posterior distribution of $\mu$ via the Bayes’ rule:

$$p(\mu |X)=\frac{p(\mu )p(X|\mu )}{p(X)}\propto p(\mu )p(X|\mu ).$$

(9)

Because

$$\begin{array}{l}p(\mu )p(X|\mu )=p(\mu )L(\mu |X)\\ \propto {\mu }^{a-1}{(1-\mu )}^{b-1}{\mu }^{{\displaystyle {\sum }_{i=1}^N{x}_i}}{\left(1-\mu \right)}^{N-{\displaystyle {\sum }_{i=1}^N{x}_i}}\\ \propto {\mu }^{a+{\displaystyle {\sum }_{i=1}^N{x}_i}-1}{(1-\mu )}^{b+N-{\displaystyle {\sum }_{i=1}^N{x}_i}-1}\end{array}$$

(10)

we have $p(\mu |X)~\mathrm{Beta}(a+{\displaystyle {\sum }_{i=1}^N{x}_i},b+N-{\displaystyle {\sum }_{i=1}^N{x}_i})$. The conclusion means that the Beta distribution is the conjugate prior for the Bernoulli likelihood.

Example 4.

Assume that random variable $x$ obeys a univariate Gaussian distribution $\mathcal{N}(\mu ,{\lambda }^{-1})$, where $\lambda$ is also named the precision. The likelihood function of $x$ for given $\mu$ is

$$\begin{array}{l}L(\lambda |\mu ,X)={\displaystyle \prod _{i=1}^N\mathcal{N}(x_i|\mu ,{\lambda }^{-1})}\\ ={\displaystyle \prod _{i=1}^N}\frac{\sqrt{\lambda }}{\sqrt{2\mathsf{\pi }}}\exp \left(-\frac{\lambda }{2}{(x_i-\mu )}^2\right)\\ ={\lambda }^{N/2}{(2\mathsf{\pi })}^{-N/2}\exp \left(-\frac{\lambda }{2}{\displaystyle {\sum }_{i=1}^N{(x_i-\mu )}^2}\right)\end{array}$$

(11)

We further suppose that the prior distribution of $\lambda$ is the Gamma distribution with parameters a and b. Let $p(\lambda |\mu ,X)$ be the posterior distribution of $\lambda$. Then we have

$$p(\lambda |\mu ,X)\propto p(\lambda )p(X,\mu |\lambda ).$$

(12)

Because

$$\begin{array}{l}p(\lambda )p(X,\mu |\lambda )\\ =\frac{1}{\Gamma (a)}{b}^a{\lambda }^{a-1}\exp (-b\lambda ){\lambda }^{N/2}{(2\pi )}^{-N/2}\exp \left(-\frac{\lambda }{2}{\displaystyle {\sum }_{i=1}^N{(x_i-\mu )}^2}\right)\\ \propto {\lambda }^{a+N/2-1}\exp \left(-\lambda \left(b+{\displaystyle {\sum }_{i=1}^N{(x_i-\mu )}^2}/2\right)\right)\end{array}$$

(13)

we get $p(\lambda |\mu ,X)~\mathrm{Gam}(a+N/2,b+{\displaystyle {\sum }_{i=1}^N{(x_i-\mu )}^2}/2)$. Therefore, the conjugate prior of the precision for a Gaussian likelihood is a Gamma distribution.

Example 5.

Assume that random vector $x$ obeys a d-dimensional Gaussian distribution $\mathcal{N}(\mathsf{\mu },{\Lambda }^{-1})$, where $\Lambda$, the inverse of covariance matrix, is called the precision matrix. We consider the case that both $\mathsf{\mu }$ and $\Lambda$ are unknown. Thus, the likelihood function for $x$ is

$$\begin{array}{l}L\left(\mathsf{\mu },\mathbf{\Lambda }|X\right)={\displaystyle {\displaystyle \prod }_{i=1}^N}\mathcal{N}\left(x_i|\mathsf{\mu }, {\mathbf{\Lambda }}^{-1}\right)\\ ={\displaystyle {\displaystyle \prod }_{i=1}^N}\frac{|\mathbf{\Lambda }{|}^{1/2}}{{(2\mathsf{\pi })}^{d/2}}\exp \left(-\frac{1}{2}{(x_i-\mathsf{\mu })}^T\mathbf{\Lambda }\left(x_i-\mathsf{\mu }\right)\right)\\ =\frac{|\mathbf{\Lambda }{|}^{N/2}}{{(2\mathsf{\pi })}^{dN/2}}\exp \left(-\frac{1}{2}{{\displaystyle \sum }}_{i=1}^N{(x_i-\mathsf{\mu })}^T\mathbf{\Lambda }\left(x_i-\mathsf{\mu }\right)\right).\end{array}$$

(14)

The prior distribution of $(\mathsf{\mu },\Lambda )$ is given by a Gaussian–Wishart distribution $\mathcal{N}\mathcal{W}(\mathsf{\mu },\Lambda |{\mathsf{\mu }}_0,\beta ,W,v)$ with the joint probability density function:

$$p(\mathsf{\mu },\mathbf{\Lambda }|{\mathsf{\mu }}_0, \beta ,W, v)=\mathcal{N}(\mathsf{\mu }|{\mathsf{\mu }}_0,\left(\beta \mathbf{\Lambda }{)}^{-1} \right)\mathcal{W}(\mathbf{\Lambda }|W, v)$$

(15)

where ${\mathsf{\mu }}_0,\beta ,W$ and $v$ are the fixed hyperparameters. Hence, it holds that

$$\begin{array}{l}p(\mathsf{\mu },\Lambda |{\mathsf{\mu }}_0,\beta ,W,v)p(X|\mathsf{\mu },\Lambda )\\ =\mathcal{N}(\mathsf{\mu }|{\mathsf{\mu }}_0,{(\beta \Lambda )}^{-1})\mathcal{W}(\Lambda |W,v)L(\mathsf{\mu },\Lambda |X)\\ \propto {|\Lambda |}^{(N+\nu -d)/2}\exp \left(-\frac{1}{2}\left({(\mathsf{\mu }-{\mathsf{\mu }}_0)}^T(\beta \Lambda )(\mathsf{\mu }-{\mathsf{\mu }}_0)+{\displaystyle {\sum }_{i=1}^N{(x_i-\mathsf{\mu })}^T\Lambda (x_i-\mathsf{\mu })+\mathrm{Tr}(W^{-1}\Lambda )}\right)\right).\end{array}$$

(16)

Denote $\overline{x}={\displaystyle {\sum }_{i=1}^N{x}_i}/N$ and $X=(x_1,x_2,\dots ,x_N)$. To obtain the above probability density function, we first derive the following formula:

$$\begin{array}{l}{(\mathsf{\mu }-{\mathsf{\mu }}_0)}^T(\beta \Lambda )(\mathsf{\mu }-{\mathsf{\mu }}_0)+{\displaystyle {\sum }_{i=1}^N{(x_i-\mathsf{\mu })}^T\Lambda (x_i-\mathsf{\mu })}\\ ={\mathsf{\mu }}^T(\beta \Lambda )\mathsf{\mu }-2{\mathsf{\mu }}_0^T(\beta \Lambda )\mathsf{\mu }+{\mathsf{\mu }}_0^T(\beta \Lambda ){\mathsf{\mu }}_0+{\displaystyle {\sum }_{i=1}^N\left({\mathsf{\mu }}^T\Lambda \mathsf{\mu }-2x_i{}^T\Lambda \mathsf{\mu }+x_i{}^T\Lambda x_i\right)}\\ =(\beta +N){\mathsf{\mu }}^T\Lambda \mathsf{\mu }-2{\left(\Lambda (\beta {\mathsf{\mu }}_0+N\overline{x})\right)}^T\mathsf{\mu }+\mathrm{Tr}\left((\beta {\mathsf{\mu }}_0{\mathsf{\mu }}_0^T+X{X}^T)\Lambda \right)\\ ={\left(\mathsf{\mu }-(\beta {\mathsf{\mu }}_0+N\overline{x})/(\beta +N)\right)}^T\left((\beta +N)\Lambda \right)\left(\mathsf{\mu }-(\beta {\mathsf{\mu }}_0+N\overline{x})/(\beta +N)\right)\\ +\mathrm{Tr}\left(\left(\beta {\mathsf{\mu }}_0{\mathsf{\mu }}_0^T+X{X}^T-(\beta {\mathsf{\mu }}_0+N\overline{x}){(\beta {\mathsf{\mu }}_0+N\overline{x})}^T/(\beta +N)\right)\Lambda \right).\end{array}$$

(17)

Then we get

$$\begin{array}{l}p(\mathsf{\mu },\Lambda |{\mathsf{\mu }}_0,\beta ,W,v)p(X|\mathsf{\mu },\Lambda )\\ \propto {|\Lambda |}^{1/2}\exp \left(-\frac{1}{2}{\left(\mathsf{\mu }-(\beta {\mathsf{\mu }}_0+N\overline{x})/(\beta +N)\right)}^T\left((\beta +N)\Lambda \right)\left(\mathsf{\mu }-(\beta {\mathsf{\mu }}_0+N\overline{x})/(\beta +N)\right)\right)\\ {|\Lambda |}^{(N+\nu -d-1)/2}\exp \left(-\frac{1}{2}\mathrm{Tr}\left(\left(W^{-1}+\beta {\mathsf{\mu }}_0{\mathsf{\mu }}_0^T+X{X}^T-(\beta {\mathsf{\mu }}_0+N\overline{x}){(\beta {\mathsf{\mu }}_0+N\overline{x})}^T/(\beta +N)\right)\Lambda \right)\right).\end{array}$$

(18)

Because $p(\mathsf{\mu },\Lambda |X)\propto p(\mathsf{\mu },\Lambda |{\mathsf{\mu }}_0,\beta ,W,v)p(X|\mathsf{\mu },\Lambda )$, we have

$$\begin{array}{l}p(\mathsf{\mu },\Lambda |X)~\\ \mathcal{N}\mathcal{W}\left(\frac{\beta {\mathsf{\mu }}_0+N\overline{\mathbf{x}}}{\beta +N},\beta +N,{\left(W^{-1}+\beta {\mathsf{\mu }}_0{\mathsf{\mu }}_0^T+X{X}^T-\frac{(\beta {\mathsf{\mu }}_0+N\overline{x}){(\beta {\mathsf{\mu }}_0+N\overline{x})}^T}{(\beta +N)}\right)}^{-1},N+v\right).\end{array}$$

(19)

This example shows that the conjugate prior of $(\mathsf{\mu },\Lambda )$ for a multivariate Gaussian $\mathcal{N}(x|\mathsf{\mu },\Lambda )$ is a Gaussian–Wishart distribution.

There are some other conclusions on the conjugate priors. For instance, given Gaussian likelihood, the conjugate prior for the mean is another Gaussian distribution, the conjugate prior for the precision matrix is a Wishart distribution. All probability distributions discussed in Section 2.1 belong to a broad class of distributions, that is, the exponential family. It is shown that the exponential family distributions have conjugate priors [54,66,72].

3. Gibbs Sampling and Variational Bayesian Inference

Due to the existence of the latent variables or unknown parameters, computing the posterior distribution is analytically intractable in general. In this section, we will provide two approximation inference methods: Gibbs sampling and variational Bayesian inference.

3.1. Gibbs Sampling

As a powerful framework for sampling from a probability distribution, Markov chain Monte Carlo (MCMC) methods, also called Markov chain simulation, construct a Markov chain such that its equilibrium distribution is the desired distribution [73,74,75,76,77,78,79,80]. Random walk Monte Carlo methods are a large subclass of MCMC and they include mainly Metropolis–Hastings [54,62,66,81], Gibbs sampling, Slice sampling [54,62,66], Multiple-try Metropolis [82,83].

Gibbs sampling or Gibbs sampler, a simple MCMC algorithm, is especially applicable for approximating a sequence of observations by a specified probability distribution when direct sampling is intractable. In addition, the basic version of Gibbs sampling is also a special case of the Metropolis–Hastings algorithm [54,62]. In detailed implementation, Gibbs sampling adopts the strategy of sampling from a conditional probability distribution instead of marginalizing the joint probability distribution by integrating over other variables. In other words, Gibbs sampling generates alternatively an instance from its corresponding conditional probability distribution by fixing other variables.

Let ${\mathbf{z}}_1,{\mathbf{z}}_2,\dots ,{\mathbf{z}}_M$ be M blocks of random variables and set $\mathbf{z}=({\mathbf{z}}_1,{\mathbf{z}}_2,\dots ,{\mathbf{z}}_M)$. The joint probability distribution of $\mathbf{z}$ is written as $p(\mathbf{z})=p({\mathbf{z}}_1,{\mathbf{z}}_2,\dots ,{\mathbf{z}}_M)$. In this subsection, we only consider the case that it is difficult to directly sampling from the joint probability distribution $p(\mathbf{z})$. To this end, we can sample each component ${\mathbf{z}}_i$ from marginal distribution in order. The marginal distribution of ${\mathbf{z}}_i$ can be theoretically obtained by the following formulation:

$$p({\mathbf{z}}_i)={\displaystyle \int {\displaystyle \int \cdots {\displaystyle \int p({\mathbf{z}}_1,{\mathbf{z}}_2,\dots ,{\mathbf{z}}_M)}}}\mathrm{d}{\mathbf{z}}_1\cdots \mathrm{d}{\mathbf{z}}_{i-1}\mathrm{d}{\mathbf{z}}_{i+1}\cdots \mathrm{d}{\mathbf{z}}_M.$$

(20)

Generally speaking, the integrating term in the above formula is not tractable. A wise sampling strategy is that we generate ${\mathbf{z}}_i$ according to the conditional probability distribution $p({\mathbf{z}}_i|\mathbf{z}\backslash {\mathbf{z}}_i)$, where $\mathbf{z}\backslash {\mathbf{z}}_i=({\mathbf{z}}_1,\dots ,{\mathbf{z}}_{i-1},{\mathbf{z}}_{i+1},\dots ,{\mathbf{z}}_M)$. By using Bayes’ rule, the relationship between the conditional probability distribution and the joint probability distribution is given as follows:

$$p({\mathbf{z}}_i|\mathbf{z}\backslash {\mathbf{z}}_i)=\frac{p({\mathbf{z}}_1,{\mathbf{z}}_2,\dots ,{\mathbf{z}}_M)}{p({\mathbf{z}}_1,\dots ,{\mathbf{z}}_{i-1},{\mathbf{z}}_{i+1},\dots ,{\mathbf{z}}_M)}\propto p({\mathbf{z}}_1,{\mathbf{z}}_2,\dots ,{\mathbf{z}}_M).$$

(21)

The sampling procedure is repeated by cycling through all variables. We summarize the outline of Gibbs sampling as below.

1. Initialize ${\mathbf{z}}_i={\mathbf{z}}_i^{(0)},i=2,\dots ,M$.

2. For $k=1,2,\dots ,T$

For $i=1,2,\dots ,M$

Generate ${\mathbf{z}}_i^{(k)}$ from the conditional probability distribution    $p({\mathbf{z}}_i|{\mathbf{z}}_1^{(k)},\dots ,{\mathbf{z}}_{i-1}^{(k)},{\mathbf{z}}_{i+1}^{(k-1)},\dots ,{\mathbf{z}}_M^{(k-1)})$.

End

End

In the above sampling procedure, T is the maximum number of iterations. Therefore, we have T samples ${\mathbf{z}}^{(k)}=({\mathbf{z}}_1^{(k)},{\mathbf{z}}_2^{(k)},\dots ,{\mathbf{z}}_M^{(k)}),k=1,2,\dots ,T$. To alleviate the influence of random initializations, we can ignore the samples within the burn-in period. Although Gibbs sampling is commonly efficient in obtaining marginal distributions from conditional distributions, it can be highly inefficient for some cases such as sampling from Mexican hat like distribution [53].

The Metropolis algorithm is an instance of MCMC algorithms and walks randomly with an acceptance/rejection rule until the convergence is achieved. As the generalization of Metropolis algorithm, the Metropolis–Hastings algorithm modifies the jumping rules and its convergence speed is improved [66]. Broadly speaking, the advanced Gibbs sampling algorithm can be regarded as a special case of the Metropolis–Hastings algorithm. The Gibbs sampler is applicable for models that are conditionally conjugate, while the Metropolis algorithm can be used for not conditionally conjugate models. Hence, we can combine both the Gibbs sampler and the Metropolis algorithm to sample from complex distributions [66].

3.2. Variational Bayesian Inference

Another widely used class of approximating the marginal probability distribution $p({\mathbf{z}}_i)$ is the variational Bayesian (VB) inference [54]. We still use the notation $X=(x_1,x_2,\dots ,x_N)$ to represent a matrix composed of $N$ observed datum and $z=(z_1,z_2,\dots ,z_M)$ to represent a vector constructed by $M$ blocks of latent variables. In VB methods, both $X$ and $z$ are assumed to be stochastic.

Given the prior distribution $p(z)$ and the data likelihood $p(X|z)$, the probability distribution of data matrix $X$, also called the evidence, can be calculated as $p(X)={\displaystyle \int p(z)p(X|z)\mathrm{d}z}$. Then we derive the conditional probability distribution of $z$ via Bayes’ rule:

$$p(z|X)=\frac{p(z)p(X|z)}{{\displaystyle \int p(z)p(X|z)\mathrm{d}z}}.$$

(22)

However, the integrating term ${\displaystyle \int p(z)p(X|z)\mathrm{d}z}$ is analytically intractable under normal conditions. Now, we harness VB methods to approximate the posterior distribution $p(z|X)$ as well as $p(X)$.

Let $q(z)$ be the trial probability distribution of $p(z|X)$. The log of the probability distribution $p(X)$ is decomposed as below:

$$\begin{array}{l}\ln p(X)={\displaystyle \int q(z)\ln p(X)}\mathrm{d}z={\displaystyle \int q(z)\ln \frac{p(X,z)}{p(z|X)}}\mathrm{d}z\\ ={\displaystyle \int q(z)\ln p(X,z)}\mathrm{d}z-{\displaystyle \int q(z)\ln p(z|X)}\mathrm{d}z\\ ={\displaystyle \int q(z)\ln p(X,z)}\mathrm{d}z-\left({\displaystyle \int q(z)\ln \frac{p(z|X)}{q(z)}}\mathrm{d}z+{\displaystyle \int q(z)\ln q(z)}\mathrm{d}z\right)\\ =L(q)+KL\left(q||p\right)\ge L(q)\end{array}$$

(23)

where $L(q)={\displaystyle \int q(z)\ln \frac{p(X,z)}{q(z)}}\mathrm{d}z$ and $KL\left(q||p\right)=-{\displaystyle \int q(z)\ln \frac{p(z|X)}{q(z)}}\mathrm{d}z$. The term $KL\left(q||p\right)$ is called the Kullback–Leibler (KL) divergence [54,84,85,86,87] between $q(z)$ and $p(z|X)$, and $L(q)$ is the lower bound of $\ln p(X)$ that achieves its lower bound if and only if $KL\left(q||p\right)=0$ (or equivalently $q=p$). The above divergence can also be explained as the information gain changing from the prior $q(z)$ to the posterior $p(z|X)$ [85].

Another perspective for Equation (23) is that,

$$\ln p(X)=\ln {\displaystyle \int p(X,z)}dz=\ln {\displaystyle \int q(z)\frac{p(X,z)}{q(z)}}dz\ge {\displaystyle \int q(z)\ln \frac{p(X,z)}{q(z)}}dz=L(q).$$

(24)

The above inequality is obtained by Jensen’s inequality. The negative lower bound $-L(q)$ is called the free energy [88,89].

The KL divergence $KL\left(q||p\right)$ can be regarded as a metric for evaluating the approximation performance of the prior distribution $q(z)$ over the posterior distribution $p(z|X)$ [54]. In the light of Equation (23), minimizing $KL\left(q||p\right)$ is equivalent to maximizing $L(q)$. What’s more, the lower bound can be further derived as below:

$$\begin{array}{l}L(q)={\displaystyle \int q(z)\ln p(X|z)}\mathrm{d}z+{\displaystyle \int q(z)\ln p(z)}\mathrm{d}z-{\displaystyle \int q(z)\ln q(z)}\mathrm{d}z\\ ={\displaystyle \int q(z)\ln p(X|z)}\mathrm{d}z+{\displaystyle \int q(z)\ln \frac{p(z)}{q(z)}}\mathrm{d}z\\ ={\mathbb{E}}_{q(z)}\ln p(X|z)-KL(q||p).\end{array}$$

(25)

The above Equation means that $\ln p(X)$ can also be regarded as the expectation of the log likelihood function $\ln p(X|z)$ with respect to $z$.

To reduce the difficulty of approximating the trial distribution $q(z)$, we partition all latent variables into M disjoint block variables $\left\{{\mathbf{z}}_1,{\mathbf{z}}_2,\dots ,{\mathbf{z}}_M\right\}$ and assume that ${\mathbf{z}}_1,{\mathbf{z}}_2,\dots ,{\mathbf{z}}_M$ are mutually independent. Based on the above assumption, we have

$$q(\mathbf{z})=q({\mathbf{z}}_1)q({\mathbf{z}}_2)\cdots q({\mathbf{z}}_M).$$

(26)

Afterwards, we utilize the mean field theory to obtain the approximation of $q(\mathbf{z})$. Concretely speaking, we update each $q({\mathbf{z}}_{\mathbf{i}})$ in turn. Let $q({\mathbf{z}}_i)$ be unknown and $q({\mathbf{z}}_j)(j\ne i)$ be fixed. Under this circumstance, we can get the optimal $q({\mathbf{z}}_i)$ by solving the following variational maximization problem:

$$q^{*}({\mathbf{z}}_i)=\arg \underset{q({\mathbf{z}}_i)}{\max }L(q(\mathbf{z}))$$

(27)

Plugging Equation (26) into the lower bound of $\ln p(X)$, we have

$$\begin{array}{l}L(q)={\displaystyle \int {\displaystyle \prod _{j=1}^{M}q({\mathbf{z}}_j)}}\ln p(X,z){\displaystyle \prod _{j=1}^M\mathrm{d}{z}_j}-{\displaystyle \int {\displaystyle \prod _{j=1}^{M}q({\mathbf{z}}_j)}}{\displaystyle \sum _{j=1}^M\ln q({\mathbf{z}}_j)}{\displaystyle \prod _{j=1}^M\mathrm{d}{\mathbf{z}}_j}\\ ={\displaystyle \int q({\mathbf{z}}_i)}\left({\displaystyle \int \ln p(X,\mathbf{z}){\displaystyle \prod _{j\ne i}q({\mathbf{z}}_j)}{\displaystyle \prod _{j\ne i}\mathrm{d}{\mathbf{z}}_j}}\right)\mathrm{d}{\mathbf{z}}_i-{\displaystyle \int q({\mathbf{z}}_i)\ln q({\mathbf{z}}_i)}\mathrm{d}{\mathbf{z}}_i+\mathrm{const}\\ ={\displaystyle \int q({\mathbf{z}}_i)\ln \tilde{p}(X,{\mathbf{z}}_i)}\mathrm{d}{\mathbf{z}}_i-{\displaystyle \int q({\mathbf{z}}_i)\ln q({\mathbf{z}}_i)}\mathrm{d}{\mathbf{z}}_i+\mathrm{const}\\ ={\displaystyle \int q({\mathbf{z}}_i)\ln \frac{\tilde{p}(X,{\mathbf{z}}_i)}{q({\mathbf{z}}_i)}}\mathrm{d}{\mathbf{z}}_i+\mathrm{const}\end{array}$$

(28)

where “const” is a term independent of ${\mathbf{z}}_i$ and $\tilde{p}(X,{\mathbf{z}}_i)$ is a new probability distribution whose log is defined by

$$\begin{array}{l}\ln \tilde{p}(X,{\mathbf{z}}_i)={\displaystyle \int \ln p(X,z){\displaystyle \prod _{j\ne i}q({\mathbf{z}}_j)}{\displaystyle \prod _{j\ne i}\mathrm{d}{z}_j}}+\mathrm{const}\\ ={\mathbb{E}}_{z\backslash z_i}\left[p(X,z)\right]+\mathrm{const}.\end{array}$$

(29)

Therefore, $\tilde{p}(X,{\mathbf{z}}_i)$ is the optimal solution of problem (27) obtained by minimizing the KL divergence between $q({\mathbf{z}}_i)$ and $\tilde{p}(X,{\mathbf{z}}_i)$.

An advantage of the VB methods is that they are immune to over fitting. But they also have some shortcomings. For example, the probability distributions derived via VB have always less probability mass in the wings than the true solution and this systematic bias may break applications. A VB method approximates a full posterior distribution by maximizing the corresponding lower bound on the marginal likelihood and it can only handle a smaller dataset. In contrast, the stochastic variational inference optimizes a subset at each iteration. This batch inference algorithm is scalable and outperforms traditional variational inference methods [90,91].

3.3. Comparisons between Gibbs Sampling and Variational Bayesian Inference

Evaluating the posterior distribution $p(z|X)$ is a central task in probabilistic models of low-rank matrix factorizations. In many practical applications, it is often infeasible to compute the posterior distribution or the expectations with respect to this distribution. Gibbs sampling and VB inference are two dominant methods for approximating the posterior distribution. The former is a stochastic approximation and the latter is a deterministic technique.

The fatal shortcoming of Expectation Maximization (EM) algorithm is that the posterior distribution of the latent variables should be given in advance. Both Gibbs sampling and VB inference can ameliorate the shortcoming of EM algorithm. Gibbs sampling is easy to implement due to the fact it adopts a Monte Carlo procedure. Another advantage of Gibbs sampling is that it can generate exact results [54]. But this method is often suitable for small-scale problems because it costs a large amount of computation. Compared with Gibbs sampling, VB inference does not generate exact results and has small computation complexity. Therefore, their strengths and weaknesses are complementary rather than competitive.

4. Probabilistic Models of Principal Component Analysis

Principal Component Analysis (PCA) is a special type of low-rank matrix factorizations. This section first introduces the classical PCA and then reviews its probabilistic models.

4.1. Principal Component Analysis

The classical PCA converts a set of samples with possibly correlated variables into another set of samples with linearly uncorrelated variables via an orthogonal transformation [1]. Based on this, PCA is an effective technique widely used in performing dimensionality reduction and extracting features.

Let $X=\left\{{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_N\right\}$ be a collection of $N$ samples with $d$ dimensions and ${\mathbf{I}}_r(r<d)$ an r-by-r identity matrix. Given a projection transformation matrix $\mathbf{W}\in {ℝ}^{d\times r}$, the PCA model can be expressed as

$${\mathbf{x}}_i=\mathbf{W}{\mathbf{z}}_i+\overline{\mathbf{x}}+{\mathsf{\epsilon }}_i,i=1,2,\dots ,N$$

(30)

where the mean vector $\overline{\mathbf{x}}={\displaystyle {\sum }_{i=1}^N{\mathbf{x}}_i}/N$, $\mathbf{W}$ satisfies ${\mathbf{W}}^T\mathbf{W}={\mathbf{I}}_r$, ${\mathbf{z}}_i$ is a representation coefficient vector with r dimensions and ${\mathsf{\epsilon }}_1,{\mathsf{\epsilon }}_2,\dots ,{\mathsf{\epsilon }}_N~\mathcal{N}(\mathbf{0}\mathbf{,}{\sigma }^2{\mathbf{I}}_d)$ are independent and identically distributed noise vectors. Denote $\mathbf{X}=\left({\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_N\right)$, $\mathbf{Z}=\left({\mathbf{z}}_1,{\mathbf{z}}_2,\dots ,{\mathbf{z}}_N\right)$, $\overline{\mathbf{X}}=\left(\overline{\mathbf{x}},\overline{\mathbf{x}},\dots ,\overline{\mathbf{x}}\right)\in {ℝ}^{d\times N}$ and $\mathbf{E}=\left({\mathsf{\epsilon }}_1,{\mathsf{\epsilon }}_2,\dots ,{\mathsf{\epsilon }}_N\right)$. Then PCA can be rewritten as the following matrix factorization formulation:

$$\mathbf{X}=\mathbf{W}\mathbf{Z}+\overline{\mathbf{X}}+\mathbf{E}.$$

(31)

According to the maximum likelihood estimation, the optimal $\mathbf{W}$ and $\mathbf{Z}$ can be obtained by solving the minimization problem:

$$\underset{\mathbf{W},\mathbf{Z}}{\min }{\Vert \mathbf{X}-\overline{\mathbf{X}}-\mathbf{W}\mathbf{Z}\Vert }_F^2,\mathrm{s}.\mathrm{t}.{\mathbf{W}}^T\mathbf{W}={\mathbf{I}}_r$$

(32)

where ${\Vert \cdot \Vert }_F$ is the Frobenious norm of a matrix. Although this optimization problem is not convex, we can obtain the optimal transform matrix ${\mathbf{W}}^{*}$ by stacking r singular vectors corresponding to the first r largest singular values of the sample covariance matrix $\mathbf{S}={\displaystyle {\sum }_{i=1}^N({\mathbf{x}}_i-\overline{\mathbf{x}})}{({\mathbf{x}}_i-\overline{\mathbf{x}})}^T/(N-1)$. Let ${\mathbf{z}}_i^{*}$ be the optimal low-dimensional representation of ${\mathbf{x}}_i$. Then it holds that ${\mathbf{z}}_i^{*}={\mathbf{W}}^{*T}({\mathbf{x}}_i-\overline{\mathbf{x}})$.

The PCA technique only supposes that the dataset is contaminated by isotropic Gaussian noise. The advantage of PCA is that it is very simple and effective in achieving the point estimations of $\mathbf{W}$ and $\mathbf{Z}$. But we can not obtain their probability distributions. In fact, the probability distributions of parameters are more useful and valuable than the point estimations in exploiting the intrinsic essence and investigating the procedure of data generation. Probabilistic models of PCA are just a class of methods for inferring the probability distributions of parameters.

4.2. Probabilistic Principal Component Analysis

In Equation (30), the low-dimensional representation $z_i$ is an unknown and deterministic parameter vector. In contrast, the original probabilistic PCA [10] regards $z_i$ as a stochastic vector. This probability model provides a general form of decomposing a d-dimensional input sample $x$:

$$x=Wz+\mathsf{\mu }+\mathsf{\epsilon }$$

(33)

where the latent random vector $z~\mathcal{N}(0,I_r)$, the noise $\mathsf{\epsilon }~\mathcal{N}(0,{\sigma }^2{I}_d)$ and $\mathsf{\mu }$ is the mean vector. In the following, a Maximum Likelihood (ML) method is proposed to obtain the point estimations of $W$, $\mathsf{\mu }$ and ${\sigma }^2$.

It is obvious that $p(x|z)=\mathcal{N}(x|Wz+\mathsf{\mu },{\sigma }^2{I}_d)$. Hence we have

$$p(x)={\displaystyle \int p(x|z)p(z)}\mathrm{d}z=\mathcal{N}(x|\mathsf{\mu },{\sigma }^2{I}_d+W{W}^T).$$

(34)

Given the observed sample set $X=\left\{{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_N\right\}$, the log-likelihood function is

$$L(W,\mathsf{\mu },{\sigma }^2|X)={\displaystyle \sum _{i=1}^N\ln p(x_i)}={\displaystyle \sum _{i=1}^N\ln }\mathcal{N}(x_i|\mathsf{\mu },{\sigma }^2{I}_d+W{W}^T).$$

(35)

The optimal $W,\mathsf{\mu },{\sigma }^2$ can be obtained by maximizing $L(W,\mathsf{\mu },{\sigma }^2|X)$. By letting the partial derivative of $L(W,\mathsf{\mu },{\sigma }^2|X)$ with respect to $\mathsf{\mu }$ be the zero vector, we can easily get the maximum likelihood estimation of $\mathsf{\mu }$: ${\mathsf{\mu }}_{\mathrm{ML}}=\overline{x}$. Hence, the optimal $W$ and ${\sigma }^2$ can be achieved by the stationary points of $L(W,{\mathsf{\mu }}_{\mathrm{ML}},{\sigma }^2|X)$ with respect to $W$ and ${\sigma }^2$.

The aforementioned method is slightly complex when computing $W$ and ${\sigma }^2$. For this purpose, an Expectation Maximization (EM) algorithm was also presented to solve probabilistic PCA. If $W$, $\mathsf{\mu }$ and ${\sigma }^2$ are given, then the joint probability distribution of $z$ and $x$ can be derived as follows:

$$\begin{array}{l}p(z,x|W,\mathsf{\mu },{\sigma }^2)=p(z|W,\mathsf{\mu },{\sigma }^2)p(x|z,W,\mathsf{\mu },{\sigma }^2)\\ =\mathcal{N}(z|0,I_r)\mathcal{N}(x|Wz+\mathsf{\mu },{\sigma }^2{I}_d).\end{array}$$

(36)

The posterior distribution of $z$ for given $x$ is

$$p(z|x,W,\mathsf{\mu },{\sigma }^2)=\frac{p(z,x|W,\mathsf{\mu },{\sigma }^2)}{p(x|W,\mathsf{\mu },{\sigma }^2)}\propto \mathcal{N}(z|0,I_r)\mathcal{N}(x|Wz+\mathsf{\mu },{\sigma }^2{I}_d).$$

(37)

Because

$$\begin{array}{l}\mathcal{N}(z|0,I_r)\mathcal{N}(x|Wz+\mathsf{\mu },{\sigma }^2{I}_d)\\ \propto \exp \left(-\frac{1}{2}{z}^{T}z\right)\exp \left(-\frac{1}{2{\sigma }^2}{\left(x-Wz-\mathsf{\mu }\right)}^T\left(x-Wz-\mathsf{\mu }\right)\right)\\ \propto \exp \left(-\frac{1}{2}{\left[z-{\left({\sigma }^{2}I+W^{T}W\right)}^{-1}{W}^T\left(x-\mathsf{\mu }\right)\right]}^T\frac{{\sigma }^{2}I+W^{T}W}{{\sigma }^2}\left[z-{\left({\sigma }^{2}I+W^{T}W\right)}^{-1}{W}^T\left(x-\mathsf{\mu }\right)\right]\right)\end{array}$$

(38)

we have

$$p(z|x,W,\mathsf{\mu },{\sigma }^2)~\mathcal{N}\left(\langle z\rangle ,{\Sigma }_z\right)$$

(39)

where ${\Sigma }_z=\mathrm{Cov}(z,z)={\sigma }^2{\left({\sigma }^{2}I+W^{T}W\right)}^{-1}$ is the covariance matrix of $z$ and $\langle z\rangle ={\Sigma }_z{W}^T\left(x-\mathsf{\mu }\right)/{\sigma }^2$ is the mean of $z$. Based on the above analysis, we give the complete-data log-likelihood function:

$$\begin{array}{l}{L}_C(W,\mathsf{\mu },{\sigma }^2|X)={\displaystyle \sum _{i=1}^N\ln p(x_i,z_i|W,\mathsf{\mu },{\sigma }^2)}\\ ={\displaystyle \sum _{i=1}^N\ln }\mathcal{N}(z_i|0,I_r)+\ln \mathcal{N}(x_i|W{z}_i+\mathsf{\mu },{\sigma }^{2}I)\\ =-\frac{dN}{2}\ln {\sigma }^2-\frac{1}{2}{\displaystyle \sum _{i=1}^N\left(z_i^T{z}_i+{\left(x_i-W{z}_i-\mathsf{\mu }\right)}^T\left(x_i-W{z}_i-\mathsf{\mu }\right)/{\sigma }^2\right)}+\mathrm{const}\\ =-\frac{dN}{2}\ln {\sigma }^2-\frac{1}{2{\sigma }^2}{\displaystyle \sum _{i=1}^N\left(z_i^T({\sigma }^{2}I+W^{T}W)z_i+{\left(x_i-\mathsf{\mu }\right)}^T\left(x_i-\mathsf{\mu }\right)+2{\left(x_i-\mathsf{\mu }\right)}^{T}W{z}_i\right)}+\mathrm{const}\end{array}$$

(40)

where “const” is a term independent of $W,\mathsf{\mu }$ and ${\sigma }^2$.

In the expectation step, we take expectation on $L_C(W,\mathsf{\mu },{\sigma }^2|X)$ with respect to $z$:

$$\begin{array}{l}{\mathbb{E}}_z\left[L_C(W,\mathsf{\mu },{\sigma }^2|X)\right]\\ =-\frac{dN}{2}\ln {\sigma }^2-\frac{1}{2{\sigma }^2}{\displaystyle \sum _{i=1}^N\left(\mathrm{tr}\left(({\sigma }^{2}I+W^{T}W)\mathbb{E}\left[z_i{z}_i^T\right]\right)+{\left(x_i-\mathsf{\mu }\right)}^T\left(x_i-\mathsf{\mu }\right)-2{\left(x_i-\mathsf{\mu }\right)}^{T}W\mathbb{E}\left[z_i\right]\right)}\\ +\mathrm{const}.\end{array}$$

(41)

According to Equation (39), we have $\mathbb{E}\left[z_i\right]={\Sigma }_z{W}^T\left(x_i-\mathsf{\mu }\right)/{\sigma }^2$ and $\mathbb{E}\left[z_i{z}_i^T\right]={\Sigma }_z+\mathbb{E}\left[z_i\right]\mathbb{E}{\left[z_i\right]}^T$.

In the maximization step, we first obtain the maximum likelihood estimation of $\mathsf{\mu }$: ${\mathsf{\mu }}_{\mathrm{ML}}=\overline{x}$ by setting the partial derivative of ${\mathbb{E}}_z\left[L_C(W,\mathsf{\mu },{\sigma }^2|X)\right]$ with respect to $\mathsf{\mu }$ be a zero vector. Similarly, we can also obtain the optimal estimations of $W$ and ${\sigma }^2$ by solving the equation set:

$$\{\begin{array}{l}\frac{\partial }{\partial W}{\mathbb{E}}_z\left[L_C(W,{\mathsf{\mu }}_{\mathrm{ML}},{\sigma }^2|X)\right]=0,\\ \frac{\partial }{\partial {\sigma }^2}{\mathbb{E}}_z\left[L_C(W,{\mathsf{\mu }}_{\mathrm{ML}},{\sigma }^2|X)\right]=0.\end{array}$$

(42)

Moreover, a mixture model for probabilistic PCA was proposed in [92]. Khan et al. replaced the Gaussian noise with Gaussian process and incorporated the information of preference pairs into the collaborative filtering [93].

4.3. Bayesian Principal Component Analysis

In probabilistic PCA, both $z$ and $\mathsf{\epsilon }$ obey Gaussian distributions, and both $W$ and $\mathsf{\mu }$ are non-stochastic parameters. Now, we further treat $w_{\cdot 1},w_{\cdot 2},\dots ,w_{\cdot r}$ and $\mathsf{\mu }$ as independent random variables, that is, $w_{\cdot i}~\mathcal{N}(0,{\alpha }_i^{-1}{I}_d)$ and $\mathsf{\mu }~\mathcal{N}(0,{\beta }^{-1}{I}_d)$. At this moment, the corresponding probabilistic model of PCA is called Bayesian PCA [11]. Set $\alpha ={({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_r)}^T$ and $\tau ={\sigma }^{-2}$. Then we give the prior distributions of $\alpha$ and $\tau$ as follows:

$$p(\alpha )={\displaystyle \prod _{i=1}^r\mathrm{Gam}({\alpha }_i|a_0,b_0)},p(\tau )=\mathrm{Gam}(\tau |c_0,d_0)$$

(43)

where $a_0,b_0,c_0,d_0$ are four given hyperparameters.

The true joint probability distribution of complete data is given by

$$\begin{array}{l}p(X,W,\alpha ,Z,\mathsf{\mu },\tau )\\ =p(W,\alpha ,Z,\mathsf{\mu },\tau )p(X|W,\alpha ,Z,\mathsf{\mu },\tau )\\ =p(\alpha )p(W|\alpha )p(Z)p(\mathsf{\mu })p(\tau )p(X|W,Z,\mathsf{\mu }).\end{array}$$

(44)

We suppose that the trial distribution of $p(X,W,\alpha ,Z,\mathsf{\mu },\tau )$ has the following form:

$$q(W,\alpha ,Z,\mathsf{\mu },\tau )=q(W)q(\alpha )q(Z)q(\mathsf{\mu })q(\tau ).$$

(45)

By making use of VB inference, we can obtain the trial probability distributions of $W,\alpha ,Z,\mathsf{\mu }$ and $\tau$ respectively. In detailed implementation, the hyperparameters $a_0,b_0,c_0,d_0,\beta$ can be set to be small positive numbers to obtain the broad priors. Unlike other approximation methods, the proposed method maximizes a rigorous lower bound.

4.4. Robust L1 Principal Component Analysis

The probabilistic model of robust L1 PCA [12] regards both $\mathsf{\mu }$ and $W$ as deterministic parameters and $\mathsf{\mu }$ is set to a zero vector without loss of the generality. We still suppose that $z$ obeys a spherical multivariate Gaussian distribution, that is, $z~\mathcal{N}(0,I_r)$. To improve the robustness, the noise ${\mathsf{\epsilon }}_i$ is assumed to follow a multivariate Laplace distribution:

$$p({\mathsf{\epsilon }}_i|\sigma )={\left(\frac{1}{2\sigma }\right)}^{d/2}\exp \left(-\frac{1}{\sigma }{\Vert {\mathsf{\epsilon }}_i\Vert }_1\right)$$

(46)

where ${\Vert \cdot \Vert }_1$ is the L1 norm of a vector. Due to the fact that the Laplace distribution has heavy tail, the proposed model in [12] is more robust against data outliers.

We use the hierarchical model to deal with the Laplace distribution. The probability density distribution of a univariate Laplace distribution $\mathrm{Lap}(0,\sigma )$ can be rewritten as

$$p(\epsilon |\sigma )={\displaystyle \int \sqrt{\frac{\beta }{2\mathsf{\pi }{\sigma }^2}}\exp \left(-\frac{\beta }{2{\sigma }^2}{\epsilon }^2\right)\frac{1}{2}{\beta }^{-2}\exp \left(-\frac{1}{2\beta }\right)d\beta }$$

(47)

Hence, we can set ${\mathsf{\epsilon }}_i~\mathcal{N}\left(0,({\sigma }^2/{\beta }_i)I_d\right)$ and ${\beta }_i~\mathrm{IGam}(1,1/2)$. Let $\rho =1/{\sigma }^2$ and give its prior distribution:

$$p(\rho )=\mathrm{Gam}(\rho |a,b)$$

(48)

where $a$ and $b$ are two given hyperparameters.

We take $Z,\mathsf{\beta }$ and $\rho$ as latent variables and $W$ as the hyperparameters. For fixed $W$, the true joint probability distribution of all latent variables is

$$p(Z,\mathsf{\beta },\rho ,X|W)=p(\rho )p(\mathsf{\beta })p(Z)p(X|W,Z).$$

(49)

The trial joint distribution of $Z,\mathsf{\beta }$ and $\rho$ is chosen as $q(Z,\mathsf{\beta },\rho )=q(Z)q(\mathsf{\beta })q(\rho )$. By applying the VB inference, the probability distributions of $Z,\mathsf{\beta }$ and $\rho$ are approximated respectively. What’s more, $W$ is updated by minimizing the robust reconstruction error of all samples.

4.5. Bayesian Robust Factor Analysis

In previous probabilistic models of PCA, the noise $\mathsf{\epsilon }$ obeys the same Gaussian distributions. However, different features maybe have different noise levels in practical applications. Now, we set $\mathsf{\epsilon }~\mathcal{N}(0,\mathrm{diag}(\mathsf{\tau }))$, where $\mathsf{\tau }={({\tau }_1,{\tau }_2,\dots ,{\tau }_d)}^T$ and $\mathrm{diag}(\mathsf{\tau })$ is a diagonal matrix. Probabilistic Factor Analysis (or PCA) [13] further assumes that $W_{ij}~\mathcal{N}(0,{\tau }_i^{-1}{\alpha }_j^{-1})$, ${\mu }_i~\mathcal{N}(0,{\tau }_i^{-1}{\beta }^{-1})$. In other words, different $W_{ij}$ or ${\mu }_i$ have different variances and the variances of $W_{ij}$ and ${\mu }_i$ also have a coupling relationship. Given $W,z,\mathsf{\mu },\mathsf{\tau }$, the conditional probability distribution of $x$ is written as

$$p(x|W,z,\mathsf{\mu },\mathsf{\tau })=\mathcal{N}(x|Wz+\mathsf{\mu },\mathrm{diag}(\mathsf{\tau })).$$

(50)

Meanwhile, the prior distributions for ${\tau }_i,{\alpha }_j$ and $\beta$ are given as follows:

$${\tau }_i~\mathrm{Gam}(a_0,b_0),{\alpha }_j~\mathrm{Gam}(c_0,d_0),\beta ~\mathrm{Gam}(e_0,f_0)$$

(51)

where $\{a_0,b_0,c_0,d_0,e_0,f_0\}$ is the set of hyperparameters.

The robust version of Bayesian factor analysis [13] considers the Student’s t-distributions instead of Gaussian noise corruptions due to the fact that the heavy tail of Student’s t-distributions makes it more robust to outliers or large sparse noise. Let $\mathsf{\epsilon }~\mathrm{St}(0,\mathsf{\tau }{.}^{-1},v)$, where $\mathsf{\tau }{.}^{-1}={(1/{\tau }_1,1/{\tau }_2,\dots ,1/{\tau }_N)}^T$ and $v=(v_1,v_2,\dots ,v_N{)}^T$. In this case, the probability distribution of $x$ for given $W,z,\mathsf{\mu },\mathsf{\tau },v$ is

$$p(x|W,z,\mathsf{\mu },\mathsf{\tau },v)=\mathrm{St}(x|Wz+\mathsf{\mu },\mathsf{\tau }{.}^{-1},v).$$

(52)

We can represent hierarchically the above Student’s t-distributions. Firstly, the conditional probability distribution of $x_k$ can be expressed as:

$$p(x_k|W,z_k,\mathsf{\mu },\mathsf{\tau },u_k)=\mathcal{N}(x_k|W{z}_k+\mathsf{\mu },\mathrm{diag}(\mathsf{\tau }{.}^{-1}.*u_k{.}^{-1}))$$

(53)

where $u_k$ is a d-dimensional column vector, “$.*$” is the Hadamard product (also known as entrywise product). Then we give the prior of $u_k$ for fixed d-dimensional hyper parameter vector $v$ as below:

$$p(u_k|v)=\mathrm{Gam}(u_k|v/2,v/2).$$

(54)

Another Bayesian robust factor analysis is on the basis of the Laplace distribution. At this moment, we assume that $\mathsf{\epsilon }~\mathrm{Lap}(0,\mathsf{\tau }{.}^{-1/2})$, where $\mathsf{\tau }{.}^{-1/2}={\left({\tau }_1^{-1/2},{\tau }_2^{-1/2},\dots ,{\tau }_N^{-1/2}\right)}^T$. In this case, we have

$$p(x|W,z,\mathsf{\mu },\mathsf{\tau })=\mathrm{Lap}(x|Wz+\mathsf{\mu },\mathsf{\tau }{.}^{-1/2})$$

(55)

The Laplace distribution generally leads to adverse effects on inferring the probability distributions of other random variables. Here, we still employ the hierarchical method, that is,

$$p(x|W,z,\mathsf{\mu },\mathsf{\tau },u)=\mathcal{N}(x|Wz+\mathsf{\mu },\mathrm{diag}(\mathsf{\tau }{.}^{-1}.*u{.}^{-1})),p(u)=\mathrm{IGam}(u|1,1/2).$$

(56)

Under this circumstance, we have

$$\begin{array}{l}p(x|W,z,\mathsf{\mu },\mathsf{\tau })={\displaystyle \int p(x,u|W,z,\mathsf{\mu },\mathsf{\tau })}\mathrm{d}u={\displaystyle \int p(x|W,z,\mathsf{\mu },\mathsf{\tau },u)p(u)}\mathrm{d}u\\ ={\displaystyle \int \mathcal{N}(x|Wz+\mathsf{\mu },\mathrm{diag}(\mathsf{\tau }{.}^{-1}.*u{.}^{-1}))\mathrm{IGam}(u|1,1/2)}\mathrm{d}u\\ \propto {\displaystyle \int {\displaystyle \prod _{i=1}^N{\left({\tau }_i{u}_i\right)}^{1/2}\exp \left(-\frac{1}{2}{\tau }_i{u}_i{(x-Wz-\mathsf{\mu })}_i^2\right)u_i{}^{-2}\exp \left(-\frac{1}{2u_i}\right)\mathrm{d}{u}_i}}\\ \propto {\displaystyle \prod _{i=1}^N{\displaystyle \int u_i^{-3/2}\exp \left(-\frac{1}{2u_i}\right)\exp \left(-\frac{1}{2}{\tau }_i{(x-Wz-\mathsf{\mu })}_i^2{u}_i\right)}}\mathrm{d}{u}_i\\ \propto {\displaystyle \prod _{i=1}^Nℒ\left\{u_i^{-3/2}\exp \left(-\frac{1}{2u_i}\right)\right\}\left(\frac{1}{2}{\tau }_i{(x-Wz-\mathsf{\mu })}_i^2\right)}\end{array}$$

(57)

where ${(x-Wz-\mathsf{\mu })}_i$ is the i-th element of vector $x-Wz-\mathsf{\mu }$. Because

$${ℒ}^{-1}\left\{\frac{\sqrt{2}}{2}\exp \left(-\sqrt{{\tau }_i{(x-Wz-\mathsf{\mu })}_i^2}\right)\right\}(t)=\frac{1}{2\sqrt{\pi }}{t}^{-3/2}\exp \left(-\frac{1}{2t}\right)$$

(58)

we have $p(x|W,z,\mathsf{\mu },\mathsf{\tau })\propto {\displaystyle \prod _{i=1}^N\exp \left(-|\sqrt{{\tau }_i}{(x-Wz-\mathsf{\mu })}_i|\right)}$. Hence, Equation (56) holds.

For the aforementioned two probabilistic models of robust factor analysis, the VB inference was proposed to approximate the posterior distributions [13]. In practice, the probability distribution of the noise should be chosen based on the application. The probabilistic models of PCA are compared in

Table 2. Probabilistic models of Principal Component Analysis (PCA) with the factorization $x=Wz+\mathsf{\mu }+\mathsf{\epsilon }$.

Probabilistic Model	Deterministic Variables	Random Variables	Prior Distributions	Solving Strategy
Probabilistic PCA [10]	$W,\mathsf{\mu }$	$\begin{array}{l}\mathsf{\epsilon }~\mathcal{N}(0,{\sigma }^2{I}_d),\hfill \\ z~\mathcal{N}(0,I_r).\hfill \end{array}$	-	ML EM
Bayesian PCA [11]	-	$\begin{array}{l}\mathsf{\epsilon }~\mathcal{N}(0,{\sigma }^2{I}_d),\hfill \\ z~\mathcal{N}(0,I_r),\hfill \\ w_{\cdot i}~\mathcal{N}(0,{\alpha }_i^{-1}{I}_d),\hfill \\ \mathsf{\mu }~\mathcal{N}(0,{\beta }^{-1}I).\hfill \end{array}$	$\begin{array}{l}p(\alpha )={\displaystyle \prod _{i=1}^r\mathrm{Gam}({\alpha }_i|a_0,b_0)},\hfill \\ p(\tau )=\mathrm{Gam}(\tau |c_0,d_0),\hfill \\ \tau ={\sigma }^2.\hfill \end{array}$	VB
Robust L1 PCA [12]	$\begin{array}{l}W,\mathsf{\mu }\hfill \\ (\mu =0)\hfill \end{array}$	$\begin{array}{l}{\mathsf{\epsilon }}_i~\mathcal{N}\left(0,({\sigma }^2/{\beta }_i)I_d\right),\hfill \\ z~\mathcal{N}(0,I_d).\hfill \end{array}$	$\begin{array}{l}{\beta }_i~\mathrm{IGam}(1,1/2),\hfill \\ p(\rho )=\mathrm{Gam}(\rho |a,b),\hfill \\ \rho =1/{\sigma }^2.\hfill \end{array}$	VB
Probabilistic factor analysis [13]	-	$\begin{array}{l}\mathsf{\epsilon }~\mathcal{N}(0,{\sigma }^2{I}_d),\hfill \\ z~\mathcal{N}(0,I_r),\hfill \\ w_{ij}~\mathcal{N}(0,{\tau }_i^{-1}{\alpha }_j^{-1}),\hfill \\ {\mu }_i~\mathcal{N}(0,{\tau }_i^{-1}{\beta }^{-1}).\hfill \end{array}$	$\begin{array}{l}{\tau }_i~\mathrm{Gam}(a_0,b_0),\hfill \\ {\alpha }_j~\mathrm{Gam}(c_0,d_0),\hfill \\ \beta ~\mathrm{Gam}(e_0,f_0).\hfill \end{array}$	VB
Bayesian robust PCA I [13]	-	$$\begin{gathered} {\mathsf{\epsilon }}_k~\mathcal{N}(0,\mathrm{diag}(\mathsf{\tau }{.}^{-1}*u_k{}^{-1})), \\ z~\mathcal{N}(0,I_r), \\ {w}_{ij}~\mathcal{N}(0,{\tau }_i^{-1}{\beta }^{-1}), \\ {\mu }_i~\mathcal{N}(0,{\tau }_i^{-1}{\beta }^{-1}). \end{gathered}$$	$\begin{array}{l}p(u_k|v)=\mathrm{Gam}(u_k|v/2,v/2),\hfill \\ {\tau }_i~\mathrm{Gam}(a_0,b_0),\hfill \\ {\alpha }_j~\mathrm{Gam}(c_0,d_0),\hfill \\ \beta ~\mathrm{Gam}(e_0,f_0).\hfill \end{array}$	VB
Bayesian robust PCA II [13]	-	$\begin{array}{l}\mathsf{\epsilon }~\mathcal{N}(0,\mathrm{diag}(\tau {.}^{-1}*u^{-1})),\hfill \\ z~\mathcal{N}(0,I_r),\hfill \\ w_{ij}~\mathcal{N}(0,{\tau }_i^{-1}{\alpha }_j^{-1}),\hfill \\ {\mu }_i~\mathcal{N}(0,{\tau }_i^{-1}{\beta }^{-1}).\hfill \end{array}$	$\begin{array}{l}p(u)=\mathrm{IGam}(u|1,1/2),\hfill \\ {\tau }_i~\mathrm{Gam}(a_0,b_0),\hfill \\ {\alpha }_j~\mathrm{Gam}(c_0,d_0),\hfill \\ \beta ~\mathrm{Gam}(e_0,f_0).\hfill \end{array}$	VB

.

5. Probabilistic Models of Matrix Factorizations

Matrix factorizations are a type of methods to approximating the data matrix by the product of two low-rank matrices. They can be regarded as a special case of PCA without considering the mean. This section will discuss the existing probabilistic models of matrix factorizations.

5.1. Matrix Factorizations

For given data matrix $X$, its low-rank factorization model is written as

$$X=WZ+E$$

(59)

where $\mathbf{W}\in {ℝ}^{d\times r}$, $\mathsf{Z}\in {ℝ}^{r\times N}$, $E$ is the noise matrix and $r<d$. We assume that $E_{ij}~\mathcal{N}(0,{\sigma }^2)$ are independent and identically distributed. We can get the optimal W and Z according to the maximum likelihood estimation. More specifically, we need to solve the following minimization problem:

$$\underset{W,Z}{\min }{\Vert X-WZ\Vert }_F.$$

(60)

The closed-form solution to problem (60) can be obtained by the Singular Value Decomposition (SVD).

To enhance the robustness to outliers or large sparse noise, we now assume that $E_{ij}~\mathrm{Lap}(0,\sigma )$. For the moment, we solve the following optimization problem:

$$\underset{W,Z}{\min }{\Vert X-WZ\Vert }_1$$

(61)

where ${\Vert \cdot \Vert }_1$ is the L1-norm of a matrix (i.e., the sum of the absolute value of all elements). This problem is also called L1-norm PCA and the corresponding optimization methods were proposed in [3,4]. Srebro and Jaakkola considered the weighted low-rank approximations problems and provided an EM algorithm [94].

5.2. Probabilistic Matrix Factorization

We still consider Gaussian noise corruptions, that is, $E_{ij}~\mathcal{N}(0,{\sigma }^2)$. Furthermore, the zero-mean spherical Gaussian priors are respectively imposed on each row of $W$ and each column of $Z$:

$$p(W|{\sigma }_W^2)={\displaystyle \prod _{i=1}^d\mathcal{N}(w_{i\cdot }|0,{\sigma }_w^2{I}_r)},p(Z|{\sigma }_Z^2)={\displaystyle \prod _{j=1}^N\mathcal{N}(z_{\cdot j}|0,{\sigma }_z^2{I}_r)}.$$

(62)

Probabilistic matrix factorization (PMF) [15] regards both ${\sigma }_w^2$ and ${\sigma }_z^2$ as two deterministic parameters. The point estimations of $W,Z$ can be obtained by maximizing the posterior distribution with the following form:

$$\begin{array}{l}p(W,Z|X,{\sigma }^2,{\sigma }_w^2,{\sigma }_z^2)\propto p(W,Z,X|{\sigma }^2,{\sigma }_w^2,{\sigma }_z^2)\\ \propto p(W|{\sigma }_w^2)p(Z|{\sigma }_z^2)p(X|W,Z,{\sigma }^2).\end{array}$$

(63)

If the priors are respectively placed on ${\sigma }_w^2$ and ${\sigma }_z^2$, we can obtain the generalized model of probabilistic matrix factorization [15]. In this case, the likelihood function is derived as

$$\begin{array}{l}p(W,Z,{\sigma }^2,{\sigma }_w^2,{\sigma }_z^2|X)\propto p(W,Z,{\sigma }^2,{\sigma }_w^2,{\sigma }_z^2,X)\\ \propto p(X|W,Z,{\sigma }^2)p(W|{\sigma }_w^2)p(Z|{\sigma }_z^2)p({\sigma }_w^2)p({\sigma }_z^2).\end{array}$$

(64)

By maximizing the above posterior distribution, we can obtain the point estimations of parameters $\{W,Z\}$ and hyperparameters $\{{\sigma }_w^2,{\sigma }_z^2\}$. Furthermore, two derivatives of PMF are also presented, i.e., PMF with a adaptive prior and PMF with constraining user-specific feature vectors.

In [89], a fully-observed variational Bayesian matrix factorization, an extension of PMF, was discussed. Meanwhile, it is shown that this new probabilistic matrix factorization can weaken the decomposability assumption and obtain the exact global analytic solution for rectangular cases.

5.3. Variational Bayesian Approach to Probabilistic Matrix Factorization

In PMF, $W_{ik}$ are independent and identically distributed and so are $Z_{kj}$. Variational Bayesian PMF [14] assumes the entries from different columns of $W$ or $Z^T$ have different variances, that is, $W_{ik}~\mathcal{N}(0,{\sigma }_k^2),$ $Z_{kj}~\mathcal{N}(0,{\rho }_k^2),$ $k=1,2,\dots ,r$. For given hyperparameters $\{{\sigma }^2,{\mathsf{\sigma }}_W^2,{\mathsf{\rho }}_Z^2\}$, we get the joint probability distribution:

$$p(X,W,Z|{\sigma }^2,{\mathsf{\sigma }}_W^2,{\mathsf{\rho }}_Z^2)=p(X|W,Z,{\sigma }^2)p(W|{\mathsf{\sigma }}_W^2)p(Z|{\mathsf{\rho }}_Z^2)$$

(65)

where ${\mathsf{\sigma }}_W^2={({\sigma }_1^2,{\sigma }_2^2,\dots ,{\sigma }_r^2)}^T$ and ${\mathsf{\rho }}_Z^2={({\rho }_1^2,{\rho }_2^2,\dots ,{\rho }_r^2)}^T$.

We assume that the trial joint distribution of $W$ and $Z$ is decomposable, that is, $q(W,Z)=q(W)q(Z)$. Using VB method, we can update alternatively $q(W)$ and $q(Z)$. The variances $\{{\sigma }^2,{\mathsf{\sigma }}_W^2,{\mathsf{\rho }}_Z^2\}$ can be determined by maximizing the following lower bound:

$$L\left(q\left(W,Z\right)\right)={\mathbb{E}}_{q\left(W\right),q\left(Z\right)}\left[\ln p\left(X,W,Z|{\sigma }^2,{\mathsf{\sigma }}_W^2,{\mathsf{\rho }}_Z^2\right)-\ln q\left(W\right)q\left(Z\right)\right].$$

(66)

The experimental results in Netflix Prize competition show that the Variational Bayesian approach has superiority over MAP and greedy residual fitting.

5.4. Bayesian Probabilistic Matrix Factorizations Using Markov Chain Monte Carlo

Variational Bayesian approach to PMF only discusses the case that $W_{ik}$ (or $Z_{kj}$) are independent and identically distributed and their means are zeros. However, Bayesian PMF [16] assumes that $w_{i\cdot }$ (or $z_{\cdot j}$) are independent and identically distributed and their mean vectors are not zero vectors. Concretely speaking, we stipulate that

$$p(W|{\mathsf{\mu }}_W,{\Lambda }_W)={\displaystyle \prod _{i=1}^d\mathcal{N}\left(w_{i\cdot }|{\mathsf{\mu }}_W,{\Lambda }_W^{-1}\right)},p(Z|{\mathsf{\mu }}_Z,{\Lambda }_Z)={\displaystyle \prod _{j=1}^N\mathcal{N}\left(z_{\cdot j}|{\mathsf{\mu }}_Z,{\Lambda }_Z^{-1}\right)}.$$

(67)

Let ${\Theta }_W=\left\{{\mathsf{\mu }}_W,{\Lambda }_W\right\},{\Theta }_Z=\left\{{\mathsf{\mu }}_Z,{\Lambda }_Z\right\}$. We further suppose the prior distributions of ${\Theta }_W$ and ${\Theta }_Z$ are Gaussian–Wishart priors:

$$\{\begin{array}{l}p\left({\Theta }_W|{\Theta }_0\right)=p\left({\mathsf{\mu }}_W|{\Lambda }_W\right)p\left({\Lambda }_W\right)=\mathcal{N}\left({\mathsf{\mu }}_W|{\mathsf{\mu }}_0,{\left({\beta }_0{\Lambda }_W\right)}^{-1}\right)\mathcal{W}\left({\Lambda }_W|W_0,v_0\right)\\ p\left({\Theta }_Z|{\Theta }_0\right)=p\left({\mathsf{\mu }}_Z|{\Lambda }_Z\right)p\left({\Lambda }_Z\right)=\mathcal{N}\left({\mathsf{\mu }}_Z|{\mathsf{\mu }}_0,{\left({\beta }_0{\Lambda }_Z\right)}^{-1}\right)\mathcal{W}\left({\Lambda }_Z|W_0,v_0\right)\end{array}$$

(68)

where ${\Theta }_0=\left\{{\mathsf{\mu }}_0,v_0,W_0\right\}$, ${\beta }_0$ is a hyper parameter.

We can initialize the parameters ${\Theta }_0$ as follows: ${\mathsf{\mu }}_0=0,v_0=d,W_0=I_d$. In theory, the predictive probability distribution of $X_{ij}^{*}$ can be obtained by marginalizing over model parameters and hyperparameters:

$$\begin{array}{l}p(X_{ij}^{*}|X,{\Theta }_0)={\displaystyle \int {\displaystyle \int {\displaystyle \int {\displaystyle \int p(X_{ij}^{*},W,Z,{\Theta }_W,{\Theta }_Z|X,{\Theta }_0)}}}}\mathrm{d}W\mathrm{d}Z\mathrm{d}{\Theta }_W\mathrm{d}{\Theta }_Z\\ ={\displaystyle \int {\displaystyle \int {\displaystyle \int {\displaystyle \int p(X_{ij}^{*}|w_{i\cdot },z_{\cdot j})p(W,Z|X,{\Theta }_W,{\Theta }_Z)p({\Theta }_W|{\Theta }_0)p({\Theta }_Z|{\Theta }_0)}}}}\mathrm{d}W\mathrm{d}Z\mathrm{d}{\Theta }_W\mathrm{d}{\Theta }_Z.\end{array}$$

(69)

However, the above integral is analytically intractable due to the fact that it is very difficult to determine the posterior distribution. Based on this, Gibbs sampling, one of the simplest Markov chain Monte Carlo, was proposed to approximate the predictive distribution $p(X_{ij}^{*}|X,{\Theta }_0)$. It is noted that MCMC methods for large-scale problems require especial care for efficient proposals and may be very slow if the sample correlation is too long.

5.5. Sparse Bayesian Matrix Completion

We consider the case that some elements of data matrix $X$ are missing and the observed index set is denoted by $\Omega$. Matrix completion assumes that $X$ is approximately low-rank and its goal is to recover all missing elements from observed elements.

For noise matrix $E$, we assume that $E_{ij}~\mathcal{N}(0,{\beta }^{-1})$. In sparse Bayesian matrix completion [19], the Gaussian distributions are imposed on two low-rank matrices, that is,

$$p(W|\mathsf{\gamma })={\displaystyle \prod _{i=1}^r\mathcal{N}\left(w_{\cdot i}|{\gamma }_i^{-1}{I}_d\right)},p(Z|\mathsf{\gamma })={\displaystyle \prod _{i=1}^r\mathcal{N}\left(z_{i\cdot }|{\gamma }_i^{-1}{I}_r\right)}.$$

(70)

Moreover, the priors of ${\gamma }_i$ are given by ${\gamma }_i~\mathrm{Gam}\left(a,b\right)$ and the prior of $\beta$ is assigned the noninformative Jeffrey’s prior: $p(\beta )\propto {\beta }^{-1}$. It is obvious that

$$p(X|W,Z,\beta )={\displaystyle \prod _{(i,j)\in \Omega }\mathcal{N}\left(X_{ij}|w_{i\cdot }{z}_{\cdot j},{\beta }^{-1}\right)}.$$

(71)

Then the joint probability distribution is

$$p(X,W,Z,\mathsf{\gamma },\beta )=p(X|W,Z,\beta )p(W|\mathsf{\gamma })p(Z|\mathsf{\gamma })p(\mathsf{\gamma })p(\beta ).$$

(72)

Let $q(W,Z,\mathsf{\gamma },\beta )$ be the approximated distribution of $p(W,Z,\mathsf{\gamma },\beta )$. The approximated procedure can be achieved by VB inference. It is demonstrated that the proposed method achieves a better prediction error in recommendation systems.

5.6. Robust Bayesian Matrix Factorization

In previous probabilistic models of matrix factorizations, there is no relationship among the variances of $W_{ik},Z_{kj},E_{ij}$. Now, we reconsider the probability distributions of $W_{ik},Z_{kj},E_{ij}$. The noise $E_{ij}$ is chosen as the heteroscedastic Gaussian scale mixture distribution:

$$E_{ij}~\mathcal{N}(0,{(\tau {\alpha }_i{\beta }_j)}^{-1})$$

(73)

and the probability distributions of $w_{i\cdot }$ and $z_{\cdot j}$ are given by:

$$p(w_{i\cdot }|{\alpha }_i)=\mathcal{N}(w_{i\cdot }|0,{({\alpha }_i{\Lambda }_w)}^{-1}),p(z_{\cdot j}|{\beta }_j)=\mathcal{N}(z_{\cdot j}|0,{({\beta }_j{\Lambda }_z)}^{-1}).$$

(74)

We also impose Gamma distribution priors on ${\alpha }_i$ and ${\beta }_j$:

$$p({\alpha }_i)=\mathrm{Gam}\left({\alpha }_i|a_0/2,b_0/2\right),p({\beta }_j)=\mathrm{Gam}\left({\beta }_j|c_0/2,d_0/2\right)$$

(75)

where $\left\{a_0,b_0,c_0,d_0\right\}$ is a given set of hyper-parameters. To reduce this problem’s complexity, we restrict ${\Lambda }_W$ and ${\Lambda }_Z$ to be diagonal matrices, that is,

$${\Lambda }_W^{-1}=\mathrm{diag}\left({\sigma }_1^2,\dots ,{\sigma }_r^2\right),{\Lambda }_Z^{-1}=\mathrm{diag}\left({\rho }_1^2,\dots ,{\rho }_r^2\right).$$

(76)

Let $\theta =\{\tau ,{\Lambda }_W,{\Lambda }_Z,a_0,b_0,c_0,d_0\}$. For the fixed parameters $\theta$, the joint probability distribution is

$$p(X,W,Z,\mathsf{\alpha },\mathsf{\beta }|\theta )=p(X|W,Z,\mathsf{\alpha },\mathsf{\beta },\theta )p(W|\mathsf{\alpha })p(Z|\mathsf{\beta })p(\mathsf{\alpha }|\theta )p(\mathsf{\beta }|\theta ).$$

(77)

We consider two types of approximated posteriors of $p(W,Z,\mathsf{\alpha },\mathsf{\beta }|X,\theta )$, one is

$$q(W,Z,\mathsf{\alpha },\mathsf{\beta })=q(W)q(Z)q(\mathsf{\alpha })q(\mathsf{\beta })$$

(78)

and another has the form

$$q(W,Z,\mathsf{\alpha },\mathsf{\beta })=q(W,\mathsf{\alpha })q(Z,\mathsf{\beta }).$$

(79)

For the above two cases, we can obtain respectively the trial probability distribution $q(W,Z,\mathsf{\alpha },\mathsf{\beta })$ by VB method [17].

In addition, a structured variational approximation was also proposed in [17]. We assume that the variational posteriors of $W$ and $Z$ obey Gaussian distributions:

$$q(W|\mathsf{\alpha })={\displaystyle \prod _{i=1}^d\mathcal{N}(w_{i\cdot }|{\overline{w}}_{i\cdot },{\alpha }_i{\overline{S}}_i)},q(Z|\mathsf{\beta })={\displaystyle \prod _{j=1}^N\mathcal{N}(z_{\cdot j}|{\overline{z}}_{\cdot j},{\beta }_j{\overline{R}}_j)}.$$

(80)

The free energy function is defined by the negative lower bound:

$$-L(q)=-{\displaystyle \int {\displaystyle \int {\displaystyle \int {\displaystyle \int q(W,Z,\mathsf{\alpha },\mathsf{\beta })\ln \frac{p(X,W,Z,\mathsf{\alpha },\mathsf{\beta }|\theta )}{q(W,Z,\mathsf{\alpha },\mathsf{\beta })}}\mathrm{d}W\mathrm{d}Z\mathrm{d}\mathsf{\alpha }\mathrm{d}\mathsf{\beta }}}}$$

(81)

By directly minimizing the free energy function with respect to ${\overline{w}}_{i\cdot },{\overline{S}}_i,{\overline{z}}_{\cdot j},{\overline{R}}_j$, we can obtain the optimal ${\overline{w}}_{i\cdot },{\overline{S}}_i,{\overline{z}}_{\cdot j},{\overline{R}}_j$. The variational posteriors of scale variables $\mathsf{\alpha }$ and $\mathsf{\beta }$ can also be recognized as the generalized inverse Gaussians.

The parameters $\theta$ can be estimated by type II maximum likelihood or empirical Bayes. In other words, we update the parameters by minimizing directly the negative lower bound. Robust Bayesian matrix factorization shows that the heavy-tailed distributions are useful to incorporate robustness information to the probabilistic models.

5.7. Probabilistic Robust Matrix Factorization

The model of probabilistic robust matrix factorization [18] considers the sparse noise corruptions. In this model, the Gaussian distributions are also imposed on $W_{ik}$ and $Z_{kj}$:

$$p(W_{ik}|{\lambda }_W)=\mathcal{N}(W_{ik}|0,{\lambda }_W^{-1}),p(Z_{kj}|{\lambda }_Z)=\mathcal{N}(Z_{kj}|0,{\lambda }_Z^{-1})$$

(82)

and the Laplace noise is placed on $E_{ij}$, that is, $E_{ij}~\mathrm{Lap}(0,\lambda )$.

From Bayes’ rule, we have

$$p(W,Z|X,\lambda ,{\lambda }_W,{\lambda }_Z)\propto p(X|W,Z,\lambda )p(W|{\lambda }_W)p(Z|{\lambda }_Z).$$

(83)

We construct the hierarchical model for the Laplace distribution:

$$p(X_{ij}|W,Z,T)=\mathcal{N}\left(X_{ij}|w_{i\cdot }{z}_{\cdot j},T_{ij}\right),p(T_{ij}|\lambda )=\mathrm{Exp}\left(T_{ij}|\lambda /2\right).$$

(84)

We regard $T$ as a latent variable matrix and denote $\theta =\left\{W,Z\right\},\widehat{\theta }=\left\{\widehat{W},\widehat{Z}\right\}$, where $\widehat{\theta }$ is the current estimation of $\theta$. An EM algorithm was proposed to inferring $W$ and $Z$ [18]. To this end, we construct the so-called Q-function:

$$Q(Z|\widehat{\theta })={\mathbb{E}}_T\left[\log p(Z|\widehat{W},X,T)|X,\widehat{\theta }\right].$$

(85)

The posterior of complete-data is

$$\begin{array}{l}p(Z|\widehat{W},X,T)=\frac{p(Z,X|\widehat{W},T)}{p(X|\widehat{W},T)}\propto p(Z,X|\widehat{W},T)\\ \propto p(X|Z,\widehat{W},T)p(Z)\end{array}$$

(86)

Hence, its log is

$$\log p(Z|\widehat{W},X,T)=-\frac{1}{2}{\displaystyle \sum _{i=1}^d{\displaystyle \sum _{j=1}^N{T}_{ij}^{-1}}}{(X_{ij}-{\widehat{w}}_{i\cdot }{z}_{\cdot j})}^2-\frac{{\lambda }_Z}{2}{\displaystyle \sum _j{z}_{\cdot j}^T}{z}_{\cdot j}+\mathrm{const}$$

(87)

where “const” is a term independent of $Z$.

To compute the expectations of T, we derive the conditional probability distribution of $T$ as follows:

$$p(T|X,\widehat{W},\widehat{Z})=\frac{p(T,X|\widehat{W},\widehat{Z})}{p(X|\widehat{W},\widehat{Z})}\propto p(T,X|\widehat{W},\widehat{Z})\propto p(T)p(X|T,\widehat{W},\widehat{Z}).$$

(88)

Hence, $T_{ij}^{-1}$ follows an inverse Gaussian distribution and its posterior expectation is given by

$$\mathbb{E}\left[T_{ij}^{-1}|X,\widehat{W},\widehat{Z}\right]=\frac{\sqrt{\lambda }}{|X_{ij}-{(\widehat{W}\widehat{Z})}_{ij}|}.$$

(89)

Thus, we get $Q(Z|\widehat{\theta })$.

To obtain the update of $Z$, we maximize the function $Q(Z|\widehat{\theta })$ with respect to $z_{\cdot j}$. By setting $\frac{\partial }{\partial z_{\cdot j}}Q(Z|\widehat{\theta })=0$, we can get the closed-form solution of $z_{\cdot j}$. The update rule for $W$ is similar to that of $Z$. The proposed probabilistic model is robust again outliers and missing data and equivalent to robust PCA under mild conditions [18].

5.8. Bayesian Robust Matrix Factorization

Another robust probabilistic model of matrix factorizations is Bayesian robust matrix factorization [20]. Gaussian distributions are still imposed on $W$ and $Z$, that is,

$$w_{i\cdot }~\mathcal{N}({\mathsf{\mu }}_W,{\Lambda }_W^{-1}),z_{\cdot j}~\mathcal{N}({\mathsf{\mu }}_Z,{\Lambda }_Z^{-1}),i=1,2,\dots ,d,j=1,2,\dots ,N.$$

(90)

We further assume both $({\mathsf{\mu }}_W,{\Lambda }_W)$ and $({\mathsf{\mu }}_Z,{\Lambda }_Z)$ follow Gaussian–Wishart distributions:

$${\Lambda }_W~\mathcal{W}(W_0,v_0),{\mathsf{\mu }}_W~\mathcal{N}({\mathsf{\mu }}_0,{({\beta }_0{\Lambda }_W)}^{-1}),{\Lambda }_Z~\mathcal{W}(W_0,v_0),{\mathsf{\mu }}_Z~\mathcal{N}({\mathsf{\mu }}_0,{({\beta }_0{\Lambda }_Z)}^{-1})$$

(91)

where $W_0,v_0,{\mathsf{\mu }}_0,{\beta }_0$ are the hyperparameters.

To enhance the robustness, we suppose the noise is the mixture of a Laplace distribution and a GIG distribution. Concretely speaking, the noise term $E_{ij}~\mathrm{Lap}(0,{\eta }_{ij})$ and the prior of ${\eta }_{ij}$ is given by $\mathrm{GIG}(p,a,b)$, where $p,a,b$ are three hyperparameters. Hence, $E_{ij}~\mathrm{Lap}(0,{\eta }_{ij})$ can be represented by two steps: $E_{ij}~\mathcal{N}(0,T_{ij}),T_{ij}~\mathrm{Exp}({\eta }_{ij}/2)$. According to the above probability distributions, we can generate ${\Lambda }_W,{\mathsf{\mu }}_W,W,{\mathsf{\mu }}_Z,{\Lambda }_Z,Z,\mathsf{\eta },T$ and $X$ in turn.

Gibbs sampling is proposed to infer the posterior distributions. For this purpose, we need to derive the posterior distributions of all random variables. Due to the fact that the derivation process of $\{{\mathsf{\mu }}_Z,{\Lambda }_Z,Z\}$ is similar to that of $\{{\mathsf{\mu }}_W,{\Lambda }_W,W\}$, the following only considers approximating the posterior distributions of $({\mathsf{\mu }}_W,{\Lambda }_W,W)$ for brevity. Firstly, the posterior distribution of $({\mathsf{\mu }}_W,{\Lambda }_W)$ is a Gaussian–Wishart distribution because that

$$\begin{array}{l}p({\mathsf{\mu }}_W,{\Lambda }_W|W,W_0,{\mathsf{\mu }}_0,v_0,{\beta }_0)\propto p({\mathsf{\mu }}_W,{\Lambda }_W,W,W_0,{\mathsf{\mu }}_0,v_0,{\beta }_0)\\ \propto {\displaystyle \prod _{i=1}^{d}p(w_{i.}|{\mathsf{\mu }}_W,{\Lambda }_W)}p({\Lambda }_W|W_0,v_0)p({\mathsf{\mu }}_W|{\mathsf{\mu }}_0,{\Lambda }_W,{\beta }_0)\end{array}$$

(92)

Then, we compute the conditional probability distribution of $w_{i.}$:

$$\begin{array}{l}p(w_{i\cdot }|X,Z,{\mathsf{\mu }}_W,{\Lambda }_W,T)\propto p(w_{i\cdot },X,Z,{\mathsf{\mu }}_W,{\Lambda }_W,T)\\ \propto p(x_{i\cdot }|w_{i\cdot }{z}_{\cdot j},t_{i\cdot })p(w_{i\cdot }|{\mathsf{\mu }}_W,{\Lambda }_W)\\ \propto {\displaystyle \prod _{j=1}^N\mathcal{N}(X_{ij}|w_{i\cdot }{z}_{\cdot j},T_{ij})}\mathcal{N}(w_{i\cdot }|{\mathsf{\mu }}_W,{\Lambda }_W).\end{array}$$

(93)

According to the above Equation, $p(w_{i\cdot }|X,Z,{\mathsf{\mu }}_W,{\Lambda }_W,T)$ is also Gaussian. Next, for given $\{X_{ij},w_{i\cdot },z_{\cdot j},{\eta }_{ij}\}$, the probability distribution of $T_{ij}$ is derived as below:

$$\begin{array}{l}p(T_{ij}|X_{ij},w_{i\cdot },z_{\cdot j},{\eta }_{ij})\propto p(T_{ij},X_{ij},w_{i\cdot },z_{\cdot j},{\eta }_{ij})\\ \propto p(X_{ij}|w_{i\cdot },z_{\cdot j},t_{ij})p(T_{ij}|{\eta }_{ij}).\end{array}$$

(94)

So, $T_{ij}|X_{ij},w_{i\cdot },z_{\cdot j},{\eta }_{ij}~\mathrm{GIG}(1/2,{\eta }_{ij},r_{ij}^2)$, where $r_{ij}=X_{ij}-w_{i\cdot }{z}_{\cdot j}$.

Finally, the posterior of ${\eta }_{ij}$ satisfies that

$$p({\eta }_{ij}|T_{ij},p,a,b)\propto p({\eta }_{ij},T_{ij},p,a,b)\propto p(T_{ij}|{\eta }_{ij})p({\eta }_{ij}|p,a,b)$$

(95)

Hence, it holds ${\eta }_{ij}|T_{ij},p,a,b~\mathrm{GIG}(p+1,T_{ij}+a,b)$. Bayesian robust matrix factorization incorporates spatial or temporal proximity in computer vision applications and batch algorithms are proposed to infer parameters.

5.9. Bayesian Model for L1-Norm Low-Rank Matrix Factorizations

For low-rank matrix factorizations, L1-norm minimization is more robust than L2-norm minimization in the presence of outliers or non-Gaussian noises. Based on this, we assume that the noise $E_{ij}$ follows the Laplace distribution: $E_{ij}~\mathrm{Lap}(0,\sqrt{\lambda /2})$. Since the Laplace noise is inconvenient for Bayesian inference, a hierarchical Bayesian model was formulated in [21]. Concretely speaking, a two-level hierarchical prior is imposed on the Laplace prior:

$$E_{ij}~\mathcal{N}(0,T_{ij}),T_{ij}~\mathrm{Exp}(\lambda ).$$

(96)

The generative models of $W$ and $Z$ are constructed as follows:

$$w_{i\cdot }~\mathcal{N}(0,{\tau }_W{I}_r),z_{\cdot j}~\mathcal{N}(0,{\tau }_Z{I}_r),i=1,2,\dots ,d,j=1,2,\dots ,N.$$

(97)

In addition, Gamma priors are placed on the precision parameters of the above Gaussian distributions:

$${\tau }_W~\mathrm{Gam}(a_0,b_0),{\tau }_Z~\mathrm{Gam}(c_0,d_0).$$

(98)

The trial posterior distribution for $\{W,Z,{\tau }_W,{\tau }_Z,T\}$ is specified as:

$$q(W,Z,{\tau }_W,{\tau }_Z,T)={\displaystyle \prod _{i=1}^{d}q}(w_{i\cdot }){\displaystyle \prod _{j=1}^{N}q(z_{\cdot j})}q({\tau }_W)q({\tau }_Z).$$

(99)

And the joint probability distribution is expressed as

$$p(W,Z,{\tau }_W,{\tau }_Z,T,X)=p(X|W,Z,T)p(W|{\tau }_W)p(Z|{\tau }_Z)p({\tau }_W)p({\tau }_Z)p(T).$$

(100)

VB inference was adopted to approximate the full posterior distribution [21]. Furthermore, varying precision parameters are also considered for different rows of $W$ or different columns of $Z$. All parameters are automatically tuned to adapt to the data, and the proposed method is applied in computer vision problems to validate its efficiency and robustness.

In

Table 3. Probabilistic models of matrix factorizations with $X=WZ+E$.

Probabilistic Model	Random Variables	Prior Distributions	Solving Strategy
PMF [15]	$\begin{array}{l}{E}_{ij}~\mathcal{N}(0,{\sigma }^2),\hfill \\ P(W|{\sigma }_W^2)={\displaystyle \prod _{i=1}^r\mathcal{N}(w_{\cdot i}|0,{\sigma }_W^2{I}_d)},\hfill \\ P(Z|{\sigma }_Z^2)={\displaystyle \prod _{j=1}^N\mathcal{N}(z_{\cdot j}|0,{\sigma }_Z^2{I}_r)}.\hfill \end{array}$	-	MAP
Variational Bayesian PMF [14]	$\begin{array}{l}{E}_{ij}~\mathcal{N}(0,{\sigma }^2),\hfill \\ W_{ik}~\mathcal{N}(0,{\sigma }_k^2),\hfill \\ Z_{kj}~\mathcal{N}(0,{\rho }_k^2).\hfill \end{array}$	-	VB
Bayesian PMF [16]	$\begin{array}{l}{E}_{ij}~\mathcal{N}(0,{\beta }^{-1}),\hfill \\ p(W|{\mathsf{\mu }}_W,{\Lambda }_W)={\displaystyle \prod _{i=1}^d\mathcal{N}\left(w_{i\cdot }|{\mathsf{\mu }}_W,{\Lambda }_W^{-1}\right)},\hfill \\ p(Z|{\mathsf{\mu }}_Z,{\Lambda }_Z)={\displaystyle \prod _{j=1}^N\mathcal{N}\left(z_{\cdot j}|{\mathsf{\mu }}_Z,{\Lambda }_Z^{-1}\right)}.\hfill \end{array}$	$\begin{array}{l}p\left({\Theta }_W|{\Theta }_0\right)=\hfill \\ \mathcal{N}\left({\mathsf{\mu }}_W|{\mathsf{\mu }}_0,{\left({\beta }_0{\Lambda }_W\right)}^{-1}\right)\mathcal{W}\left({\Lambda }_W|W_0,v_0\right),\hfill \\ p\left({\Theta }_Z|{\Theta }_0\right)=\hfill \\ \mathcal{N}\left({\mathsf{\mu }}_Z|{\mathsf{\mu }}_0,{\left({\beta }_0{\Lambda }_Z\right)}^{-1}\right)\mathcal{W}\left({\Lambda }_Z|W_0,v_0\right).\hfill \end{array}$	Gibbs sampling
Sparse Bayesian matrix completion [19]	$\begin{array}{l}{E}_{ij}~\mathcal{N}(0,{\beta }^{-1}),\hfill \\ p(W|\mathsf{\gamma })={\displaystyle \prod _{i=1}^r\mathcal{N}\left(w_{\cdot i}|{\gamma }_i^{-1}{I}_d\right)},\hfill \\ p(Z|\mathsf{\gamma })={\displaystyle \prod _{i=1}^r\mathcal{N}\left(z_{i\cdot }|{\gamma }_i^{-1}{I}_r\right)}.\hfill \end{array}$	$\begin{array}{l}{\gamma }_i~\mathrm{Gam}\left(a,b\right),\hfill \\ p(\beta )\propto {\beta }^{-1}.\hfill \end{array}$	VB
Robust Bayesian matrix factorization [17]	$\begin{array}{l}{E}_{ij}~\mathcal{N}(0,{(\tau {\alpha }_i{\beta }_j)}^{-1}),\hfill \\ p(w_{i\cdot }|{\alpha }_i)=\mathcal{N}(w_{i\cdot }|0,{({\alpha }_i{\Lambda }_w)}^{-1}),\hfill \\ p(z_{\cdot j}|{\beta }_j)=\mathcal{N}(z_{\cdot j}|0,{({\beta }_j{\Lambda }_z)}^{-1}).\hfill \end{array}$	$\begin{array}{l}p({\alpha }_i)=\mathrm{Gam}\left({\alpha }_i|a_0/2,b_0/2\right),\hfill \\ p({\beta }_j)=\mathrm{Gam}\left({\beta }_j|c_0/2,d_0/2\right).\hfill \end{array}$	VB, type II ML, empirical Bayes
Probabilistic robust matrix factorization [18]	$\begin{array}{l}{E}_{ij}~\mathcal{N}\left(X_{ij}|0,T_{ij}\right),\hfill \\ p(W_{ik}|{\lambda }_W)=\mathcal{N}(W_{ik}|0,{\lambda }_W^{-1}),\hfill \\ p(Z_{kj}|{\lambda }_Z)=\mathcal{N}(Z_{kj}|0,{\lambda }_Z^{-1}).\hfill \end{array}$	$T_{ij}~\mathrm{Exp}\left(\lambda /2\right).$	EM
Bayesian robust matrix factorization [20]	$\begin{array}{l}{E}_{ij}~\mathcal{N}\left(0,T_{ij}\right),\hfill \\ w_{i\cdot }~\mathcal{N}({\mathsf{\mu }}_W,{\Lambda }_W^{-1}),\hfill \\ z_{\cdot j}~\mathcal{N}({\mathsf{\mu }}_Z,{\Lambda }_Z^{-1}).\hfill \end{array}$	$\begin{array}{l}{T}_{ij}~\mathrm{Exp}({\eta }_{ij}/2),\hfill \\ {\Lambda }_W~\mathcal{W}(W_0,v_0),\hfill \\ {\mathsf{\mu }}_W~\mathcal{N}({\mathsf{\mu }}_0,{({\beta }_0{\Lambda }_W)}^{-1}),\hfill \\ {\Lambda }_Z~\mathcal{W}(W_0,v_0),\hfill \\ {\mathsf{\mu }}_Z~\mathcal{N}({\mathsf{\mu }}_0,{({\beta }_0{\Lambda }_Z)}^{-1}).\hfill \end{array}$	Gibbs sampling
Bayesian model for L1-norm low-rank matrix factorizations [21]	$\begin{array}{l}{E}_{ij}~\mathcal{N}(0,T_{ij}),\hfill \\ w_{i\cdot }~\mathcal{N}(0,{\tau }_W{I}_r),\hfill \\ z_{\cdot j}~\mathcal{N}(0,{\tau }_Z{I}_r).\hfill \end{array}$	$\begin{array}{l}{T}_{ij}~\mathrm{Exp}(\lambda ),\hfill \\ {\tau }_W~\mathrm{Gam}(a_0,b_0),\hfill \\ {\tau }_Z~\mathrm{Gam}(c_0,d_0).\hfill \end{array}$	VB

, we list all probabilistic models of matrix factorizations discussed in this section, and compare their probability distributions, priors and solving strategy.

6. Probabilistic Models of Robust PCA

Compared with traditional PCA, robust PCA is more robust to outliers or large sparse noise. Stable robust PCA [95], a stable version of robust PCA, decomposes the data matrix into the sum of a low-rank matrix, a sparse noise matrix and a dense noise matrix. The low-rank matrix obtained by solving a relaxed principal component pursuit is simultaneously stable to small noise and robust to gross sparse errors. This section will review probabilistic models of robust PCA.

6.1. Bayesian Robust PCA

In [22], a stable robust PCA is modeled as:

$$X=W\left(D\Lambda \right)Z+B.*S+E$$

(101)

where $D=\mathrm{diag}(d_{11},d_{22},\dots ,d_{rr})$, $\Lambda =\mathrm{diag}({\lambda }_{11},{\lambda }_{22},\dots ,{\lambda }_{rr})$, $B\in {\{0,1\}}^{d\times N}$. The three terms $W\left(D\Lambda \right)Z$, $B.*S$ and $E$ are the low-rank, the sparse noise and the dense noise terms respectively. If we do not consider the sparse noise term $B.*S$, then Equation (101) is equivalent to Equation (59). If all columns of $B.*S$ are same, then the stable robust PCA becomes to be the PCA model (31).

The following considers the probability distributions of all matrices in the right of Equation (101). We assume that $w_{\cdot i}~\mathcal{N}(0,I_d/d),d_{ii}~\mathrm{Bern}(p_i),{\lambda }_{ii}~\mathcal{N}(0,{\tau }^{-1}),z_{\cdot j}~\mathcal{N}(0,I_r/r)$, $b_{\cdot j}~{\displaystyle \prod _{k=1}^r\mathrm{Bern}}({\pi }_k)$, $s_{\cdot j}~\mathcal{N}(0,v^{-1}{I}_d)$, $E_{ij}~\mathcal{N}(0,{\gamma }^{-1})$, $i=1,2,\dots ,r,j=1,2,\dots ,N$. The priors of $p_i$, $\tau$, ${\pi }_k$, $v$ and $\gamma$ are given by $p_i~\mathrm{Beta}({\alpha }_0,{\beta }_0),$ $\tau ~\mathrm{Gam}(a_0,b_0)$, ${\pi }_k~\mathrm{Beta}({\alpha }_1,{\beta }_1)$, $v~\mathrm{Gam}(c_0,d_0)$ and $\gamma ~\mathrm{Gam}(e_0,f_0)$ respectively.

Two methods were proposed in [22], that is, Gibbs sampling and VB inference. For the second method, the joint probability distribution is

$$\begin{array}{l}p\left(X,W,Z,D,\Lambda ,B,S,p,\tau ,\mathsf{\pi },v,\gamma \right)\\ =p\left(X|W,Z,D,\Lambda ,B,S,\gamma \right)p\left(W\right)p\left(Z\right)p\left(\tau \right)p\left(\Lambda |\tau \right)p\left(D|p\right)p\left(p\right)\\ p\left(\mathsf{\pi }\right)p\left(B|\mathsf{\pi }\right)p\left(v\right)p\left(S|v\right)p\left(\gamma \right).\end{array}$$

(102)

VB inference is employed to approximating the posterior distribution $p\left(W,Z,D,\Lambda ,B,S,p,\tau ,\mathsf{\pi },v,\gamma |X\right)$. Bayesian robust PCA is robust to different noise levels without changing model hyperparameters and exploits additional structure of the matrix in video applications.

6.2. Variational Bayes Approach to Robust PCA

In Equation (101), if both $D$ and $\Lambda$ are set to be identical matrices, then we have another form of stable robust PCA:

$$X=WZ+B.*S+E.$$

(103)

Essentially, both Equations (101) and (103) are equivalent. Formally, Equation (103) can also be transformed into another formula:

$$X=WZ+B.*S+(1-B).*E.$$

(104)

We assume that $w_{\cdot i}~\mathcal{N}({\mathsf{\mu }}_0,{\sigma }_0^2{I}_d)$, $z_{j\cdot }~\mathcal{N}(v_0,{\sigma }_0^2{I}_N)$, $b_i{}_j~\mathrm{Bern}(b_0)$. Let the means of $S_{ij}$ and $E_{ij}$ be zeros and their precision be ${\tau }_S$ and ${\tau }_E$ respectively. The priors of ${\tau }_S$ and ${\tau }_E$ are given by

$${\tau }_S~\mathrm{Gam}({\alpha }_S,{\beta }_S),{\tau }_E~\mathrm{Gam}({\alpha }_E,{\beta }_E).$$

(105)

A naïve VB approach was proposed in [23]. The trial distribution is stipulated as

$$q(W,Z,B,{\tau }_S,{\tau }_E)=q(W)q(Z)q(B)q({\tau }_S)q({\tau }_E).$$

(106)

The likelihood function is constructed as

$$L=p(X|W,Z,B,{\tau }_S,{\tau }_E)={\displaystyle \prod _{i,j}\mathcal{N}{(X_{ij}-w_{i\cdot }{z}_{\cdot j}|{\tau }_S)}^{B_{ij}}}\mathcal{N}{(X_{ij}-w_{i\cdot }{z}_{\cdot j}|{\tau }_E)}^{1-B_{ij}}.$$

(107)

The prior distribution is represented by

$$\pi (W,Z,B,{\tau }_S,{\tau }_E)=p(W|{\mathsf{\mu }}_0,{\sigma }_0^2{I}_d)p(Z|v_0,{\sigma }_0^2{I}_N)p(B)p({\tau }_S)p({\tau }_E).$$

(108)

We first construct a function: $G(q)={\mathbb{E}}_q\left[\log L\right]-D_{KL}(q||\pi )$. To simplify this problem, we let

$$q(w_{i\cdot })~\mathcal{N}({\mathsf{\mu }}_{wi},{\Sigma }_{wi}),q(z_{\cdot i})~\mathcal{N}({\mathsf{\mu }}_{zi},{\Sigma }_{zi}).$$

(109)

To find the updates for $W,Z,$ we can maximize the function $G$ with respect to ${\mathsf{\mu }}_{wi},{\Sigma }_{wi},{\mathsf{\mu }}_{zi},{\Sigma }_{zi}$ respectively. The main advantage of the proposed approach is that it can incorporate additional prior information and cope with missing data.

6.3. Sparse Bayesian Robust PCA

In Equation (101), we replace $B.*S$ by $S$ for the sake of simplicity. Thus, we have a model of sparse Bayesian robust PCA [19]. Assume that $w_{\cdot k}~\mathcal{N}(0,{\gamma }_k^{-1}{I}_d),z_{k\cdot }~\mathcal{N}(0,{\gamma }_k^{-1}{I}_N)$, $S_{ij}~\mathcal{N}(0,{\alpha }_{ij}^{-1})$ and $E_{ij}~\mathcal{N}(0,{\beta }^{-1})$, where $k=1,2,\dots ,r,i=1,2,\dots ,d,j=1,2,\dots ,N$. The priors of ${\gamma }_k$ are given by ${\gamma }_k~\mathrm{Gam}(a,b)$. What’s more, we assign Jeffrey’s priors to ${\alpha }_{ij}$ and $\beta$:

$$p(\beta )={\beta }^{-1},p({\alpha }_{ij})={\alpha }_{ij}^{-1},i=1,2,\dots ,d,j=1,2,\dots ,N.$$

(110)

The joint distribution is expressed as

$$p(X,W,Z,S,\mathsf{\gamma },\mathsf{\alpha },\beta )=p(X|W,Z,S,\beta )p(W|\mathsf{\gamma })p(Z|\mathsf{\gamma })p(S|\mathsf{\alpha })p(\mathsf{\gamma })p(\mathsf{\alpha })p(\beta ).$$

(111)

VB inference was used to approximate the posterior distributions of all variables matrices [19]. Experimental results in video background/foreground separation show that the proposed method is more effective than MAP and Gibbs sampling.

7. Probabilistic Models of Non-Negative Matrix Factorization

Non-negative matrix factorization (NMF) decomposes a non-negative data matrix into the product of two non-negative low-rank matrices. Mathematically, we can formulate NMF as follows:

$$X=WZ+E$$

(112)

where $X_{ij},W_{ik},Z_{kj}$ are non-negative. Multiplicative algorithms [96] are often used to obtain the point estimations of both $W$ and $Z$. This section will introduce probabilistic models of NMF.

7.1. Probabilistic Non-Negative Matrix Factorization

Equation (112) can be rewritten as $X\approx WZ={\displaystyle {\sum }_{k=1}^r{w}_{\cdot k}{z}_{k\cdot }}$. Let ${\theta }_k=\left\{w_{\cdot k},z_{k\cdot }\right\}$ and $\theta =\{{\theta }_1,{\theta }_2,\dots ,{\theta }_r\}$. Probabilistic non-negative matrix factorization [24] introduces a generative model:

$$X_{ij}={\displaystyle \sum _{k=1}^r{C}_{k,ij}},C_{k,ij}~p(C_{k,ij}|{\theta }_k).$$

(113)

The probability distributions of $C_{k,ij}$ can be assumed to follow Gaussian or Poisson distributions because they are closed under summation. This assumption means that we can get easily the probability distributions of $X_{ij}$. Four algorithms were proposed in [24], that is, multiplicative, EM, Gibbs sampling and VB algorithms.

7.2. Bayesian Inference for Nonnegative Matrix Factorization

For arbitrary $k\in \{1,2,\dots ,r\}$, Bayesian non-negative matrix factorization [25] introduces variables $S_k=\{S_{ikj}|i=1,2,\dots ,d,j=1,2,\dots ,N\}$ as latent sources. The hierarchical model of $X_{ij}$ is given by

$$S_{ikj}~\mathrm{Poiss}(W_{ik}{Z}_{kj}),X_{ij}={\displaystyle \sum _{k=1}^r{S}_{ikj}}.$$

(114)

In view of the fact that a Gamma distribution is the conjugate prior to Poisson distribution, the hierarchical priors of $W_{ik}$ and $Z_{kj}$ are proposed as follows

$$W_{ik}~\mathrm{Gam}(a_{ik}^W,b_{ik}^W),Z_{kj}~\mathrm{Gam}(a_{ik}^Z,b_{ik}^Z).$$

(115)

Let $\Theta =\left\{a_{ik}^W,b_{ik}^W,a_{ik}^Z,b_{ik}^Z\right\}$ and $S=\{S_1,S_2,\dots ,S_r\}$. For given $\Theta$ and $X$, the posterior distribution is expressed as $p(W,Z,S|X,\Theta )$. We assume the trial distribution of $p(W,Z,S|X,\Theta )$ is factorable: $q(W,Z,S)=q(W)q(Z)q(S)$. Both VB inference and Gibbs sampling were proposed to infer the probability distributions of all variables [25].

The above Bayesian nonnegative matrix factorization is not a matrix factorization approach to latent Dirichlet allocation. To this end, another Bayesian extension of the nonnegative matrix factorization algorithm was proposed in [27]. What’s more, Paisley et al. also provided a correlated nonnegative matrix factorization based on the correlated topic model [27]. The stochastic variational inference algorithms were presented to solve the proposed two models.

7.3. Bayesian Nonparametric Matrix Factorization

Gamma process nonnegative matrix factorization (GaP-NMF) was developed in [26]. This Bayesian nonparametric matrix factorization considers the case that the number of sources $r$ is unknown. Let non-negative hidden variable ${\theta }_k$ be the overall gain of the k-th source and $L$ a large number of sources. We assume that

$$W_{ik}~\mathrm{Gam}(a,a),Z_{kj}~\mathrm{Gam}(b,b),X_{ij}~\mathrm{Exp}({\displaystyle \sum _{k=1}^r{\theta }_k{W}_{ik}{Z}_{kj}}),{\theta }_k~\mathrm{Gam}(\alpha /L,\alpha c).$$

(116)

The posterior distribution is expressed as $p(W,Z,\mathsf{\theta }|X,a,b,\alpha ,c)$ for given hyperparameters $a,b,\alpha ,c$. The trial distribution of $\mathsf{\theta },W,Z$ is assumed to be factorable, that is, $q(W,Z,\mathsf{\theta })=q(W)q(Z)q(\mathsf{\theta })$. The flexible generalized inverse-Gaussian distributions are imposed on $W_{ik},Z_{kj}$ and ${\theta }_k$ respectively,

$$q(W_{ik})~\mathrm{GIG}({\gamma }_{ik}^W,{\rho }_{ik}^W,{\tau }_{ik}^W),q(Z_{kj})~\mathrm{GIG}({\gamma }_{kj}^Z,{\rho }_{kj}^Z,{\tau }_{kj}^Z),q({\theta }_k)~\mathrm{GIG}({\gamma }_k^{\mathsf{\theta }},{\rho }_k^{\mathsf{\theta }},{\tau }_k^{\mathsf{\theta }}).$$

(117)

The lower bound of the marginal likelihood is computed as

$$\begin{array}{l}\log p\left(X,a,b,\alpha ,c\right)\ge {\mathbb{E}}_q\left[\log p(X|W,Z,\mathsf{\theta })\right]\\ +{\mathbb{E}}_q\left[\log p(W|a)\right]-{\mathbb{E}}_q\left[\log q\left(W\right)\right]\\ +{\mathbb{E}}_q\left[\log p(Z|b)\right]-{\mathbb{E}}_q\left[\log q\left(Z\right)\right]\\ +{\mathbb{E}}_q\left[\log p(\mathsf{\theta }|\alpha ,c)\right]-{\mathbb{E}}_q\left[\log q\left(\mathsf{\theta }\right)\right]\end{array}$$

(118)

By maximizing the lower bound of $\log p(X,a,b,\alpha ,c)$, we can yield an approximation distribution of $q(W,Z,\mathsf{\theta })$. GaP-NMF is applied in recorded music and the number of latent sources is discovered automatically.

7.4. Beta-Gamma Non-Negative Matrix Factorization

For the moment, we do not factorize the data matrix $X$ into the product of two low-rank matrices. Different from previous probabilistic models of NMF, we assume that $X_{ij}$ is generated from a Beta distribution: $\mathrm{Beta}(A_{ij},B_{ij})$. For two matrices $A$ and $B$, we jointly factorize them as

$$A\approx CH,B\approx DH$$

(119)

where $C\in {ℝ}_{+}^{d\times r},D\in {ℝ}_{+}^{d\times r},H\in {ℝ}_{+}^{r\times N}$.

In Beta-Gamma non-negative matrix factorization [28], the generative model is given by

$$X_{ij}~\mathrm{Beta}({\displaystyle \sum _{k=1}^r{C}_{ik}{H}_{kj}},{\displaystyle \sum _{k=1}^r{D}_{ik}{H}_{kj}}),C_{ik}~\mathrm{Gam}({\mu }_0,{\alpha }_0),D_{ik}~\mathrm{Gam}({\upsilon }_0,{\beta }_0),H_{kj}~\mathrm{Gam}({\rho }_0,{\varsigma }_0)$$

(120)

Variational inference framework was adopted and a new lower-bound was proposed to approximate the objective function to derive an analytically tractable approximate solution of the posterior distribution. Beta-Gamma non-negative matrix factorization is used in source separation, collaborative filtering and cancer epigenomics analysis.

8. Other Probabilistic Models of Low-Rank Matrix/Tensor Factorizations

Besides the probabilistic models discussed in foregoing sections, there are many other types of probabilistic low-rank matrix factorization models. A successful application of probabilistic low-rank matrix factorization is the collaborative filtering in recommendation systems. In collaborative filtering, there are several other Poisson models in which the observations are usually modeled with a Poisson distribution, and these models mainly include [97,98,99,100,101,102,103,104,105]. As a matter of fact, the Poisson factorization roots in the nonnegative matrix factorization and takes advantage of the sparse essence of user behavior data and scales [103].For some probabilistic models with respect to collaborative filtering, the Poisson distribution is changed into other probability distributions and this change deals with logistic function [106,107,108], Heaviside step function [107,109], Gaussian cumulative density function [110] and so on. In addition, side information on the a low-dimensional latent presentations is integrated into probabilistic low-rank matrix factorization models [111,112,113], and the case that the data is missing not at random is taken into consideration [109,114,115].

It is worthy to pay attention to other applications of probabilistic low-rank matrix factorization models. For instance, [116] developed a probabilistic model for low-rank subspace clustering. In [88], a sparse additive matrix factorization was proposed by a Bayesian regularization effect and the corresponding model was applied into a foreground/background video separation problem.

Recently, probabilistic low-rank matrix factorizations have been extended into the case of tensor decompositions (factorizations). Tucker decomposition and CP decomposition are two popular tensor decomposition approaches. The probabilistic Tucker decomposition models mainly include probabilistic Tucker decomposition [34], exponential family tensor factorization [35] and InfTucker model [36]. Probabilistic Tucker decomposition was closely related to probabilistic PCA. In [35], an integration method was proposed to model heterogeneously attributed data tensors. InfTucker, a tensor-variate latent nonparametric Bayesian model, conducted Tucker decomposition in an infinite feature space.

More probabilistic models of tensor factorizations focus on CP tensor decomposition model. For example, Ermis and Cemgil investigated variational inference for probabilistic latent tensor factorization [37]. Based on hierarchical dirichlet process, a Bayesian probabilistic model for unsupervised tensor factorization was proposed [38]. In [39], a novel probabilistic tensor factorization was proposed by extending probabilistic matrix factorization. A probabilistic latent tensor factorization was proposed in [40] to address the task of link pattern prediction. Based on the Polya-Gamma augmentation strategy and online expectation maximization algorithm, [41] proposed a scalable probabilistic tensor factorization framework. As the generalization of Poisson matrix factorization, Poisson tensor factorization was presented in [42]. In [43], a Bayesian tensor factorization models was proposed to infer the latent group structures from dynamic pairwise interaction patterns. In [44], a Bayesian non-negative tensor factorization model was presented for count-valued tensor data and scalable inference algorithms were developed. A scalable Bayesian framework for low-rank CP decomposition was presented and it can analyses both continuous and binary datasets [45]. A zero-truncated Poisson tensor factorization for binary tensors was proposed in [46]. A Bayesian robust tensor factorization [47] was proposed and it is the extension of probabilistic stable robust PCA. And in [48], the CP factorization was formulated by a hierarchical probabilistic model.

9. Conclusions and Future Work

In this paper, we have made a survey on probabilistic models of low-rank matrix factorizations and the related works. To classify the main probabilistic models, we divide low-rank matrix factorizations into several groups such as PCA, matrix factorizations, robust PCA, non-negative matrix factorization and so on. For each category, we list representative probabilistic models, describe the probability distributions of all random matrices or latent variables, present the corresponding inference methods and compare their similarity and difference. Besides, we further provide an overview of probabilistic models of low-rank tensor factorizations and discuss other probabilistic matrix factorizations models.

Although probabilistic low-rank matrix/tensor factorizations have made some progresses, we still face some challenges in theories and applications. Future research may concern the following aspects:

• Scalable algorithms to infer the probability distributions and parameters. Although both Gibbs sampling and variational Bayesian inference have their own advantages, they need large computation cost for real large-scale problems. A promising future direction is to design scalable algorithms.
• Constructing new probabilistic models of low-rank matrix factorizations. It is necessary to develop other probabilistic models according to the actual situation. For example, we can consider different types of sparse noise and different probability distributions (including the prior distributions) of low-rank components or latent variables.
• Probabilistic models of non-negative tensor factorizations. There is not much research on this type of probabilistic models. Compared with probabilistic models of tensor factorizations, the probabilistic non-negative tensor factorizations models are more complex and difficult in inferring the posterior distributions.
• Probabilistic TT format. In contrast to both CP and Tucker decompositions, the TT format provides stable representations and is formally free from the curse of dimensionality. Hence, probabilistic model of the TT format would be an interesting research issue.