A brief introduction to Bayesian inference, expectation maximization and variational Bayesian inference

• Maximizing log likelihood and maximizing a posterior

In signal processing field, parameter estimation is frequently discussed. Classical method is to use the maximum likelihood (ML) estimation, which can be expressed as following:

$\hat{\mathcal{x}}=arg~~ max ~~log(p(\mathcal{y}; \mathcal{x}))$

where $p(\mathcal{y};\mathcal{x})$ is the probability of the observation $\mathcal{y}$ given the variable $\mathcal{x}$.

Often, there is prior information about the variable to estimate, which can be expressed as the prior probability distribution  $p(\mathcal{x})$. According to the Bayesian rule, posterior of variable $\mathcal{x}$ can be written as

$p(\mathcal{x}|\mathcal{y})=p(\mathcal{x})p(\mathcal{y}|\mathcal{x})/p(\mathcal{y})$ 

With the prior probability distribution, the Maximum A Posterior (MAP) can be formulated as:

$\hat{\mathcal{x}}=arg~~max~~log(p(\mathcal{x}|\mathcal{y}))$

Note that here $p(\mathcal{y};\mathcal{x})$ implicitly means that the variable $\mathcal{x}$ is a non-random variable, while $p(\mathcal{y}|\mathcal{x})$ denotes that the variable $\mathcal{x}$ is a random variable. Although in papers these two expressions are interchangely used, it is better to make them clear here.

Moreover, when $\mathcal{x}$ is not a random variable, the problem is termed parameter estimation problem; if the variable is a random variable, the problem is usually termed parameter inference problem.

• Problems with ML and MAP

It is known that ML is an unbiased. However, some times it is hard to maximize the likelihood. In MAP, it is uneasy to make sure that the prior distribution could make the estimation unbiased. 

If one models the probability distribution of $\mathcal{x}$ as $p(\mathcal{x})$, usually the distribution is determined by some parameters $\theta$. More specifically, the distribution of $\mathcal{x}$ should be written as $p(\mathcal{x};\theta)$. In this way, the MAP estimation can be rewritten as

$\hat{\mathcal{x}}=arg~~max~log( p(\mathcal{x};\theta)p(\mathcal{y}|\mathcal{x}))$

From above equation, one can see that in MAP estimation, the hyperparameter $\theta$ (since we term $\mathcal{x}$ as parameter) is assumed to be fixed beforehand. A natural question is that: why make the hyperparameter $\theta$ fixed? Why not estimate the value of $\theta$ from the observation? 

• Expectation Maximization method

One simple method to tackle with $\theta$ is to integrate out the variable $\mathcal{x}$, which can be written as

$p(y;\theta)=\int{p(\mathcal{y}|\mathcal{x})p(\mathcal{x};\theta)}d\mathcal{x}$ (empirical Bayesian estimation)

Maximizing the above equation is often termed empirical Bayesian estimation because the hyperparameter $\theta$ is estimated from the observation. The problem with the maximization of $p(y;\mathcal{\theta})$ is the difficulty in calculation. 

The way to solve the problem empirical Bayesian estimation problem is to use Expectation Maximization (EM) algorithm. In EM algorithm, it is assumed that $\mathcal{y}$ and $\mathcal{x}$ are together complete observations. As we have no access to $\mathcal{x}$, $\mathcal{x}$ is regarded as missing data or missing observation. Since the complete observation is $\mathcal{x}$ and $\mathcal{y}$, the maximum likelihood should be

$\hat{\mathcal{x}}=arg~~max~~log(p(\mathcal{x},\mathcal{y};\theta))$

However, we do not know the distribution $p(\mathcal{x},\mathcal{y})$. The way to deal with this problem is as follows:

 

$log(p(\mathcal{y};\theta))=log(\int{p(\mathcal{x},\mathcal{y};\theta)d\mathcal{x}})$

$~~~~~~~~~~~~~~~~~~=log(\int{p(\mathcal{x},\mathcal{y};\theta)q(\mathcal{x})/q(\mathcal{x})}d\mathcal{x})$

$~~~~~~~~~~~~~~~~~~\geq \int{q(\mathcal{x})log(p(\mathcal{x},\mathcal{y};\theta)/q(\mathcal{x}))}d\mathcal{x}$

 

where the upper inequality is according to Jenson's inequality. The inequality becomes equality only when $q(\mathcal{x})=p(\mathcal{x}|\mathcal{y})$. Instead of maximizing $log(p(\mathcal{y};\theta))$, one focus on its lower bound. Define

$g(\theta, \hat{\theta})=\int{p(\mathcal{x}|\mathcal{y};\hat{\theta})log(p(\mathcal{x},\mathcal{y};\theta)/p(\mathcal{x}|\mathcal{y};\hat{\theta}))}$

The EM algorithm is as follows:

1. Expectation step: assign $\hat{\theta}$ with the lastest value; calculate $g(\theta,\hat{\theta}^{(t)})$ (t is index of current iteration)

2. Maximization step: update the value of $\theta$ by maximizing $g(\theta,\hat{\theta}^{(t)})$

The essence of EM algorithm is Majorization Maximization technique. This can be viewed in the figure below.

From the above figure, one can see that using $p(\mathcal{y};\theta)$ one can construct a surrogate function $g(\theta,\hat{\theta})$. The E-step and M-step in EM algorithm are actually iteratively maximizing $p(\mathcal{y};\theta)$.

According to the EM algorithm, one could easily varify that:

 

$p(\mathcal{y};\hat{\theta}^{(t+1)})\geq g(\hat{\theta}^{(t+1)},\hat{\theta}^{(t)})=max~g(\theta,\hat{\theta}^{(t)})\geq g(\hat{\theta}^{(t)},\hat{\theta}^{(t)})=p(\mathcal{y};\hat{\theta}^{(t)})$

 

which means that EM algorithm indeed can increase the value of $p(\mathcal{y};\theta)$.

EM algorithm terminates either the number of iterations pass the pre-defined value or the difference of $\theta$ between subsequent iterations are below a pre-defined threshold. After get the estimation of $\theta$, which is $\hat{\theta}$, by EM algorithm, one can easily calculate $p(\mathcal{x};\hat{\theta})$. The value of $\mathcal{x}$ can be estimated by

$\hat{x}=arg~~max~log(p(\mathcal{y}|\mathcal{x})p(\mathcal{x};\hat{\theta}))$

There is another way to derive EM algorithm as follows.

 

$log p(\mathcal{y};\theta)=log(\frac{p(\mathcal{x},\mathcal{y};\theta)}{p(\mathcal{x}|\mathcal{y};\theta)})$

$~~~~~~~~~~~~~~~~=\int{q(\mathcal{x})log(\frac{p(\mathcal{x},\mathcal{y};\theta)q(\mathcal{x})}{p(\mathcal{x}|\mathcal{y};\theta)q(\mathcal{x})})}d\mathcal{x}$

$~~~~~~~~~~~~~~~~=\int{q(\mathcal{x})log(\frac{p(\mathcal{x},\mathcal{y};\theta)}{q(\mathcal{x})})}d\mathcal{x}-\int{q(\mathcal{x})log(\frac{p(\mathcal{x}|\mathcal{y};\theta)}{q(\mathcal{x})})}d\mathcal{x}$

$~~~~~~~~~~~~~~~~=F(q,\theta)+KL(q||p)$

 

where $F(q,\theta)=\int{q(\mathcal{x})log(\frac{p(\mathcal{x},\mathcal{y};\theta)}{q(\mathcal{x})})}d\mathcal{x}$ and $KL(p||q)$ is the Kullback-Leibler divergence between $p(\mathcal{x}|\mathcal{y};\theta)$ and $q(\mathcal{x})$. Since $KL(p||q)$ is nonnegative, $F(q,\theta)$ is the lower bound of the log likelihood $log(p(\mathcal{y};\theta))$. When $q(\mathcal{x})=p(\mathcal{x}|\mathcal{y};\theta)$, one has $KL(P||q)=0$ and therefore $F(q,\theta)=log(p(\mathcal{y};\theta))$.

Our original objective is to maximize the log likelihood. Here, we try to maximize the lower bound of the log likelihood. Assume that we have an estimate of $\theta$ as $\hat{\theta}$. Let $q(\mathcal{x})=p(\mathcal{x}|\mathcal{y};\hat{\theta})$. The lower bound can be written as

$F(q,\theta)=\int{p(\mathcal{x}|\mathcal{y};\hat{\theta})log(\frac{p(\mathcal{x},\mathcal{y};\theta)}{p(\mathcal{x}|\mathcal{y};\hat{\theta})})}d\mathcal{x}$

$~~~~~~~~~~~=\int{p(\mathcal{x}|\mathcal{y};\hat{\theta})log(p(\mathcal{x},\mathcal{y};\theta))}d\mathcal{x}-\int{p(\mathcal{x}|\mathcal{y};\hat{\theta})log(p(\mathcal{x}|\mathcal{y};\hat{\theta}))}d\mathcal{x}$

$~~~~~~~~~~~=g(\theta,\hat{theta})+C$

 

where $C$ is the entropy of the distribution $p(\mathcal{x}|\mathcal{y};\hat{\theta})$.

One can see that to maximize the lower bound $F(q, \theta)$, $g(\theta, \hat{\theta})$ needs to be maximized. Therefore, EM algorithm is employed to update the value of $\theta$.

• Variational Bayesian Inference

In the aforementioned parts, it is assumed that the expression of $p(\mathcal{x}|\mathcal{y};\theta)$ can be given. In the last part, we make an assumption that $q(\mathcal{x})=p(\mathcal{x}|\mathcal{y};\theta)$. However, in practice, what if the expression of $p(\mathcal{x}|\mathcal{y};\theta)$ is not available? In this case, variation Bayesian inference is ready to sovle this problem.

In the aforementioned parts, $q(\mathcal{x})$ is an unknown function. To get its expression, functional derivative is needed. Here, variational method is employed to find the function. This is done by assuming that the function over which optimization is performed has specific forms. In Bayesian inference, the widely assumed form is the factorized form, which means that $q(\mathcal{x})$ can be factorized as 

$q(\mathcal{x})=\prod q_{i}(\mathcal{x}_{i})$

where $\mathcal{x}_{i}$ is the $i$th element of $\mathcal{x}$.

In this case, the lower bound $F(q,\theta)$ can be expressed as

$F(q,\theta)=\int{\prod_{i}q_{i}[log(p(\mathcal{x},\mathcal{y};\theta))-\sum_{i}log(q_{i})]}d\mathcal{x}$

$~~~~~~~~~~~=\int{\prod_{i} q_{i}log(p(\mathcal{x},\mathcal{y};\theta))\prod_{i}d\mathcal{x}_{i}}-\sum_{i}\int{\prod_{j}q_{j}log(q_{j})d\mathcal{x}_{i}}$

$~~~~~~~~~~~=\int{q_{j}[\prod_{i\neq j}p(\mathcal{x}_{i})log(p(\mathcal{x},\mathcal{y}|\theta))d\mathcal{x}_{i}]}d\mathcal{x}_{j}-\int{q_{j}log(q_{j})d\mathcal{x}_{j}}-\sum_{i\neq j}\int{q_{i}log(q_{i})}d\mathcal{x}_{i}$

$~~~~~~~~~~~=\int{q_{j}f(\mathcal{x}_{j},\mathcal{y},\theta)}d\mathcal{x}_{j}-\int{q_{j}log(q_{j})}d\mathcal{x}_{j}-\sum_{i\neq j}\int{q_{i}log(q_{i})}d\mathcal{x}$

$~~~~~~~~~~~=-KL(q_{j}|| f(\mathcal{x}_{j},\mathcal{y},\theta))-\sum_{i\neq j}\int{q_{i}log(q_{i})}d\mathcal{x}$

where $f(\mathcal{x}_{j},\mathcal{y};\theta)=\int{log(p(\mathcal{x},\mathcal{y};\theta))}\prod_{i\neq j}q_{i}d\mathcal{x}_{i}$.

Clearly, the bound in the upper equation is maximized when the Kullback-Leibler becomes zeos, which means that $q_{j}(\mathcal{x}_{j}=f(\mathcal{x}_{j},\mathcal{y};\theta))$. 

• Graphic models

In real applicaitons, when one wishes to model the problem as a Bayesian inference problem, it is better to draw a graphic model. For instance, the linear measurement equation is 

$\mathcal{t}=\mathcal{\Phi w}+\mathcal{\epsilon}$

where $\mathcal{w}$ is the vector we are to estimate, $\mathcal{\epsilon}$ is the measurement noise and $\mathcal{t}$ is the measurement signal. The problem is to recovery $\mathcal{w}$ with given $\mathcal{\Phi}$ and $\mathcal{t}$.

Let $\beta$ be the variance of the i.i.d Gaussian noise $\mathcal{\epsilon}$. There are many ways to model this problem in Bayesian inference way, such as:

1. Model $\mathcal{w}$ and $\mathcal{\epsilon}$ as non-random variable. This problem is equal to least square regression problem.
2. Model $\mathcal{w}$ as i.i.d Gaussian distribution with non-random variance. In this case complete observations should be $\mathcal{t}$ and $\mathcal{w}$. Variatioanl Bayesian inference can be employed to estimate the variance of $\mathcal{w}$ along with $\mathcal{w}$.
3. Model $\mathcal{w}$ as i.i.d Gaussian distribution with variance being Gamma distribution. Model $\beta$ as Gamma distribution. In this case the complete observation should be $\mathcal{t}$, $\mathcal{t}$, $\beta$ and the variance of $\mathcal{w}$.

All the three models are illustrated in the following figure, where sqaure represents non-random variable and circle represents random variable. For complete procedure for estimation, please refer to the following paper:

Dimitris G. Tzikas, etc. The variational approximation for Bayesian inference. IEEE Signal Processing Magazine, 2008.

The rule of thumb is that: after you draw the graphic model, always treat the random variable as observations or hidden observations. Write down the likelihood function and use variational Bayesian inference to estimate the hyper-parameter. 

• Reference

1. Tzikas, D.G., etc. The variational approximation for Bayesian inference, IEEE Signal Processing Magazine, vol. 25, no. 6, pp: 131-146, 2008.
2. Chuong B. Do, etc. What is the expectation maximization algorithm? Nature Biotechnology, no. 26, pp: 897-899, 2008
3. Wei Wu, etc. Bayesian machine learning: EEG/MEG signal processing measurement. IEEE Signal Processing Magazine, vol. 33, no. 1, pp:14-36, 2016.