Introduction of the variational inference method

Visualize VI solution to bivariate Gaussian

• Variational inference minimizes KL$(q\parallel p)$: mean of the approximation is correct, but variance (along the orthogonal direction) is significantly under-estimated

• Expectation propagation minimizes KL$(p\parallel q)$: solution equals marginal distributions

Figure 1: Left: variational inference. Right: expectation propagation

Example: univariate Gaussian

Model selection

Mixture of Gaussians

• For each observation $x_n\in R^D$, we have a corresponding latent variable $z_n$, a 1-of-$K$ binary group indicator vector

• Mixture of Gasussians joint likelihood, based on $N$ observations $$\begin{array}{cc}p(Z\mid \pi )& =\prod _{n=1}^N\prod _{k=1}^K{\pi }_k^{z_{nk}}\\ p(X\mid Z,\mu ,\Lambda )& =\prod _{n=1}^N\prod _{k=1}^K\text{N}{(x_n\mid {\mu }_k,{\Lambda }_k^{-1})}^{z_{nk}}\end{array}$$

Figure 2: Graph representation of mixture of Gaussians

Conjugate priors

• Dirichlet for $\pi$ $$p\left(\pi \right)=\text{Dir}(\pi \mid {\alpha }_0)\propto \prod _{k=1}^K{\pi }_k^{{\alpha }_{0k}-1}$$

• Independent Gaussian-Wishart for $\mu ,\Lambda$ $$\begin{array}{cc}p(\mu ,\Lambda )& =\prod _{k=1}^{K}p({\mu }_k\mid {\Lambda }_k)p\left({\Lambda }_k\right)\\ & =\prod _{k=1}^K\text{N}({\mu }_k\mid m_0,({\beta }_0{\Lambda }_k{)}^{-1})\text{W}({\Lambda }_k\mid W_0,{\nu }_0)\end{array}$$

  • Usually, the prior mean $m_0=0$

Variational distribution

• Joint posterior $$p(X,Z,\pi ,\mu ,\Lambda )=p(X\mid Z,\mu ,\Lambda )p(Z\mid \pi )p\left(\pi \right)p(\mu \mid \Lambda )p\left(\Lambda \right)$$

• Variational distribution factorizes between the latent variables and the parameters $$\begin{array}{cc}q(Z,\pi ,\mu ,\Lambda )& =q\left(Z\right)q(\pi ,\mu ,\Lambda )\\ & =q\left(Z\right)q\left(\pi \right)\prod _{k=1}^{K}q({\mu }_k,{\Lambda }_k)\end{array}$$

Variational solution for $Z$

• Optimized factor $$\begin{array}{cc}\log q^{\ast }\left(Z\right)& =E_{\pi ,\mu ,\Lambda }\left[\log p\right(X,Z,\pi ,\mu ,\Lambda \left)\right]\\ & =E_{\pi }\left[\log p\right(Z\mid \pi \left)\right]+E_{\mu ,\Lambda }\left[\log p\right(X\mid Z,\mu ,\Lambda \left)\right]\\ & =\sum _{n=1}^N\sum _{k=1}^K{z}_{nk}\log {\rho }_{nk}+\text{const}\\ \log {\rho }_{nk}=& E\left[\log {\pi }_k\right]+\frac{1}{2}E\left[\log \left|{\Lambda }_k\right|\right]-\frac{D}{2}\log \left(2\pi \right)\\ & -\frac{1}{2}{E}_{\mu ,\Lambda }\left[{(x_n-{\mu }_k)}^{\prime }{\Lambda }_k(x_n-{\mu }_k)\right]\end{array}$$

• Thus, the factor $q^{\ast }\left(Z\right)$ takes the same functional form as the prior $p(Z\mid \pi )$ $$q^{\ast }\left(Z\right)=\prod _{n=1}^N\prod _{k=1}^K{r}_{nk}^{z_{nk}},r_{n}k=\frac{{\rho }_{nk}}{{\sum }_{j=1}^K{\rho }_{nj}}$$

  • By $q^{\ast }\left(Z\right)$, the posterior mean (i.e., responsibility) $E\left[z_{nk}\right]=r_{nk}$

Define three statistics wrt the responsibilities

• For each of group $k=1,\dots ,K$, denote $$\begin{array}{cc}{N}_k& =\sum _{n=1}^N{r}_{nk}\\ {\overline{x}}_k& =\frac{1}{N_k}\sum _{n=1}^N{r}_{nk}{x}_n\\ S_k& =\frac{1}{N_k}\sum _{n=1}^N{r}_{nk}(x_n-{\overline{x}}_k){(x_n-{\overline{x}}_k)}^{\prime }\end{array}$$

Variational solution for $\pi$

• Optimized factor $$\begin{array}{cc}\log q^{\ast }\left(\pi \right)& =\log p\left(\pi \right)+E_Z\left[p\right(Z\mid \pi \left)\right]\\ & =({\alpha }_0-1)\sum _{k=1}^K\log {\pi }_k+\sum _{k=1}^K\sum _{n=1}^N{r}_{nk}\log {\pi }_{nk}+\text{const}\end{array}$$

• Thus, $q^{\ast }\left(\pi \right)$ is a Dirichlet distribution $$q^{\ast }\left(\pi \right)=\text{Dir}\left(\alpha \right),{\alpha }_k={\alpha }_0+N_k$$

Variational solution for ${\mu }_k,{\Lambda }_k$

• Optimized factor for $({\mu }_k,{\Lambda }_k)$ $$\begin{array}{cc}\log q^{\ast }({\mu }_k,{\Lambda }_k)=& E_Z\left[\sum _{n=1}^N{z}_{nk}\log \text{N}(x_n\mid {\mu }_k,{\Lambda }_k^{-1})\right]\\ & +\log p({\mu }_k\mid {\Lambda }_k)+\log p\left({\Lambda }_k\right)\end{array}$$

• Thus, $q^{\ast }({\mu }_k,{\Lambda }_k)$ is Gaussian-Wishart $$\begin{array}{ccc}{q}^{\ast }({\mu }_k\mid {\Lambda }_k)& =\text{N}(m_k,{\left({\beta }_k{\Lambda }_k\right)}^{-1})q^{\ast }\left({\Lambda }_k\right)& =\text{W}({\Lambda }_k\mid W_k,{\nu }_k)\end{array}$$

• Parameters are updated by the data $$\begin{array}{cc}{\beta }_k& ={\beta }_0+N_k,m_k=\frac{1}{{\beta }_k}({\beta }_0{m}_0+N_k{\overline{x}}_k),{\nu }_k={\nu }_0+N_k\\ W_k^{-1}& =W_0^{-1}+N_k{S}_k+\frac{{\beta }_0{N}_k}{{\beta }_0+N_k}({\overline{x}}_k-m_0){({\overline{x}}_k-m_0)}^{\prime }\end{array}$$

Similarity between VI and EM solutions

• Optimization of the variational posterior distribution involves cycling between two stages analogous to the E and M steps of the maximum likelihood EM algorithm

  • Finding $q^{\ast }\left(Z\right)$: analogous to the E step, both need to compute the responsibilities
  • Finding $q^{\ast }(\pi ,\mu ,\Lambda )$: analogous to the M step

• The VI solution (Bayesian approach) has little computational overhead, comparing with the EM solution (maximum likelihood approach). The dominant computational cost for VI are

  • Evaluation of the responsibilities
  • Evaluation and inversion of the weighted data covariance matrices

Advantage of the VI solution over the EM solution:

• Since our priors are conjugate, the variational posterior distributions have the same functional form as the priors

1. No singularity arises in maximum likelihood when a Gassuain component “collapses” onto a specific data point

   • This is actually the advantage of Bayesian solutions (with priors) over frequentist ones

2. No overfitting if we choose a large number $K$. This is helpful in determining the optimal number of components without performing cross validation

   • For ${\alpha }_0<1$, the prior favors soutions where some of the mixing coefficients $\pi$ are zero, thus can result in some less than $K$ number components having nonzero mixing coefficients

Computing variational lower bound

• To test for convergence, it is useful to monitor the bound during the re-estimation.
• At each step of the iterative re-estimation, the value of the lower bound should not decrease $$\begin{array}{cc}L=& \sum _Z\iiint q^{\ast }(Z,\pi ,\mu ,\Lambda )\log \left\{\frac{p(X,Z,\pi ,\mu ,\Lambda )}{q^{\ast }(Z,\pi ,\mu ,\Lambda )}\right\}d\pi d\mu d\Lambda \\ =& E\left[\log p\right(X,Z,\pi ,\mu ,\Lambda \left)\right]-E\left[\log q^{\ast }\right(Z,\pi ,\mu ,\Lambda \left)\right]\\ =& E\left[\log p\right(X\mid Z,\mu ,\Lambda \left)\right]+E\left[\log p\right(Z\mid \pi \left)\right]\\ & +E\left[\log p\right(\pi \left)\right]+E\left[\log p\right(\mu ,\Lambda \left)\right]\\ & -E\left[\log q^{\ast }\right(Z\left)\right]-E\left[\log q^{\ast }\right(\pi \left)\right]-E\left[\log q^{\ast }\right(\mu ,\Lambda \left)\right]\end{array}$$

Label switching problem

• EM solution of maximum likelihood does not have label switching problem, because the initialization will lead to just one of the solutions

• In a Bayesian setting, label switching problem can be an issue, because the marginal posterior is multi-modal.

• Recall that for multi-modal posteriors, variational inference usually approximate the distribution in the neighborhood of one of the modes and ignore the others

Induced factorizations

• Induced factorizations: the additional factorizations that are a consequence of the interaction between

  • the assumed factorization, and
  • the conditional independence properties of the true distribution

• For example, suppose we have three variation groups $A,B,C$

  • We assume the following factorization $$q(A,B,C)=q(A,B)q\left(C\right)$$
  • If $A$ and $B$ are conditional independent $$A\perp B\mid X,C⟺p(A,B\mid X,C)=p(A\mid X,C)p(B\mid X,C)$$ then we have induced factorization $q^{\ast }(A,B)=q^{\ast }\left(A\right)q^{\ast }\left(B\right)$ $$\begin{array}{cc}\log q^{\ast }(A,B)& =E_C\left[\log p\right(A,B\mid X,C\left)\right]+\text{const}\\ & =E_C\left[\log p\right(A\mid X,C\left)\right]+E_C\left[\log p\right(B\mid X,C\left)\right]+\text{const}\end{array}$$

Bayesian linear regression

• Here, I use a denotion system commonly used in statistics textbooks. So its different from the one used in this book.

• Likelihood function $$p(y\mid \beta )=\prod _{n=1}^N\text{N}(y_n\mid x_n\beta ,{\varphi }^{-1})$$

  • $\varphi =1/{\sigma }^2$ is the precision parameter. We assume that it is known.
  • $\beta \in R^p$ includes the intercept

• Prior distributions: Normal Gamma $$\begin{array}{cc}p(\beta \mid \kappa )& =\text{N}(\beta \mid 0,{\kappa }^{-1}I)\\ p\left(\kappa \right)& =\text{Gam}(\kappa \mid a_0,b_0)\end{array}$$

Variational solution for $\kappa$

• Variational posterior factorization $$q(\beta ,\kappa )=q\left(\beta \right)q\left(\kappa \right)$$

• Varitional solution for $\kappa$ $$\begin{array}{cc}\log q^{\ast }\left(\kappa \right)& =\log p\left(\kappa \right)+E_{\beta }\left[\log p\right(\beta \mid \kappa \left)\right]\\ & =(a_0-1)\log \kappa -b_0\kappa +\frac{p}{2}\log \kappa -\frac{\kappa }{2}E\left[{\beta }^{\prime }\beta \right]\end{array}$$

• Varitional posterior is a Gamma $$\begin{array}{cc}\kappa & \sim \text{Gam}(a_N,b_N)\\ a_N& =a_0+\frac{p}{2}\\ b_N& =b_0+\frac{E\left[{\beta }^{\prime }\beta \right]}{2}\end{array}$$

Variational solution for $\beta$

• Variational solution for $\beta$ $$\begin{array}{cc}\log q^{\ast }\left(\beta \right)& =\log p(y\mid \beta )+E_{\kappa }\left[\log p\right(\beta \mid \kappa \left)\right]\\ & =-\frac{\varphi }{2}{(y-X\beta )}^2-\frac{E\left[\kappa \right]}{2}{\beta }^{\prime }\beta \\ & =-\frac{1}{2}{\beta }^{\prime }(E\left[\kappa \right]I+\varphi X^{\prime }X)\beta +\varphi {\beta }^{\prime }{X}^{\prime }y\end{array}$$

• Variational posterior is a Normal $$\begin{array}{cc}\beta & \sim \text{N}(m_N,S_N)\\ S_N& ={(E\left[\kappa \right]I+\varphi X^{\prime }X)}^{-1}\\ m_N& =\varphi S_N{X}^{\prime }y\end{array}$$

Iteratively re-estimate the variational solutions

• Means of the variational posteriors $$\begin{array}{cc}E\left[\kappa \right]& =\frac{a_N}{b_N}\\ E\left[{\beta }^{\prime }\beta \right]& =m_N{m}_N^{\prime }+S_N\end{array}$$

• Lower bound of $\log p\left(y\right)$ can be used in convergence monitoring, and also model selection $$\begin{array}{cc}L=& E\left[\log p\right(\beta ,\kappa ,y\left)\right]-E\left[\log q^{\ast }\right(\beta ,\kappa \left)\right]\\ =& E_{\beta }\left[\log p\right(y\mid \beta \left)\right]+E_{\beta ,\kappa }\left[\log p\right(\beta \mid \kappa \left)\right]+E_{\kappa }\left[\log p\right(\kappa \left)\right]\\ & -E_{\beta }\left[\log q^{\ast }\right(\beta \left)\right]-E_{\kappa }\left[\log q^{\ast }\right(\kappa \left)\right]\end{array}$$

TO BE CONTINUED

References

• Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Springer.