It turns out that we can perform (approximate) Bayesian inference by minimizing the Variational Free Energy (VFE) with respect to the variational distribution $q$ .
To explain inference by VFE minimization (abbreviated as VFEM), we first rewrite the VFE in terms of "inference bound" minus "log-evidence" terms (the bound-evidence (BE) decomposition):
$$\begin{align*}
F[q] &= \underbrace{ \int q(z) \log \frac{q(z)}{p(z)} \mathrm{d}z }_{\text{complexity}} - \underbrace{\int q(z) \log p(x|z) \mathrm{d}z}_{\text{accuracy}} \\
&= \underbrace{\int q(z) \log \frac{q(z)}{p(z|x)}\mathrm{d}z}_{\text{inference bound}\geq 0} - \underbrace{\log p(x)}_{\text{log-evidence}}
\end{align*}$$
persist_js_state$60a50f12-063c-4247-9c9c-bb41d5dc9811running§runtime ݸdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$60a50f12-063c-4247-9c9c-bb41d5dc9811depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&Dghas_pluto_hook_features¤body:
persist_js_state$26c75b5e-d294-11ef-173e-b3f46a1df536running§runtime tdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$26c75b5e-d294-11ef-173e-b3f46a1df536depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&ƥӷhas_pluto_hook_features¤body
Exam guide: Working out VFE minimization for the GMM to these update equations (eqs B-10.58 through B-10.63) is not something that you need to reproduce without assistance at the exam. Rather, the essence is that it is possible to arrive at closed-form variational update equations for the GMM. You should understand though how FEM works conceptually and in principle be able to derive variational update equations for very simple models that do not involve clever mathematical tricks.
persist_js_state$fd233604-120a-4838-a660-d5021bccecd0running§runtimeڀdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$fd233604-120a-4838-a660-d5021bccecd0depends_on_disabled_cells¦queued¤logsoutputmimeimage/svg+xmlrootassigneelast_run_timestampAڹ&has_pluto_hook_features¤body ;
persist_js_state$c18b7c1b-8011-469b-ad92-7d50c23c46e3running§runtime :9depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$c18b7c1b-8011-469b-ad92-7d50c23c46e3depends_on_disabled_cells¦queued¤logsoutputmime!application/vnd.pluto.tree+objectrootassigneeconst plot_limslast_run_timestampAڹ&Yhas_pluto_hook_features¤bodyobjectidae2ec98d0122aedctypeNamedTupleelementsxlimobjectid7b3b6e4d0ce70e0atypeTupleelements1.3text/plain5.5text/plain!application/vnd.pluto.tree+objectylimobjectidfc6c5f93bd4de540typeTupleelements41text/plain99text/plain!application/vnd.pluto.tree+objectpersist_js_state$bbdca8c2-022f-42be-bcf7-80d86f7f269crunning§runtime 7depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$bbdca8c2-022f-42be-bcf7-80d86f7f269cdepends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&ihas_pluto_hook_features¤body
Let us reconsider the Bound-Evidence decomposition of the VFE for a model $p(x,z)$ with variational posterior $q(z)$ ,
$$\begin{align}
\mathrm{F}[q] = \underbrace{\sum_z q(z) \log \frac{q(z)}{p(z|x)}}_{\text{inference bound}\geq 0} \underbrace{- \log p(x)}_{\text{surprise}} \tag{BE}
\end{align}$$
The VFE comprises two cost terms:
The surprise (or negative log-evidence), $-\log p(x)$ , reflects the cost of predicting the data $x$ using a model $p(x, z)$ , assuming that (ideal) Bayesian inference can be performed. Specifically, the evidence $p(x)$ is obtained from the joint model $p(x, z)$ by marginalizing over the latent variables:
$$p(x) = \sum_z p(x,z) \,.$$
$$\sum_z q(z) \log \frac{q(z)}{p(z |x)} \geq 0 \,,$$
quantifies the cost of imperfect inference, i.e., the discrepancy between the variational posterior $q(z)$ and the true Bayesian posterior $p(z | x)$ .
In any practical setting, using a model implies performing inference within that model. Therefore, the effective cost of applying a model is not merely the surprise but also must include the cost of inference.
Put more bluntly: a model with very high Bayesian evidence $p(x)$ may still be practically unusable due to exorbitant inference costs.
In the literature, the VFE is typically interpreted as an approximation (more precisely, an upper-bound) to the surprise, $-\log p(x)$ , which is often regarded as the “true” measure of model performance. However, we argue that this perspective should be reversed: the VFE should be considered the true performance metric in practice, as it accounts for both model fit and the tractability of inference. The surprise can be viewed as a special case of the VFE, corresponding to a zero inference bound, that only applies when ideal Bayesian inference is computationally feasible.
persist_js_state$cc547bfa-a130-4382-af47-73de56e4741brunning§runtimeqdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$cc547bfa-a130-4382-af47-73de56e4741bdepends_on_disabled_cells¦queued¤logsoutputmimetext/plainrootassigneeold_faithfullast_run_timestampAڹ&Mbhas_pluto_hook_features¤body272×2 Matrix{Float64}:
3.6 79.0
1.8 54.0
3.333 74.0
2.283 62.0
4.533 85.0
2.883 55.0
4.7 88.0
⋮
4.75 75.0
4.117 81.0
2.15 46.0
4.417 90.0
1.817 46.0
4.467 74.0persist_js_state$62868f61-95c7-4e07-854f-a171aadc667brunning§runtimeͽѸdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$62868f61-95c7-4e07-854f-a171aadc667bdepends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&xUhas_pluto_hook_features¤bodyH
Code Example: EM-algorithm for the GMM on the Old-Faithful data setpersist_js_state$26c5d734-d294-11ef-20a3-afd2c3324323running§runtime Fdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$26c5d734-d294-11ef-20a3-afd2c3324323depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&±Bhas_pluto_hook_features¤bodyj
GMM is a Flexible Model
GMMs are very popular models. They have decent computational properties and are universal approximators of densities (as long as there are enough Gaussians, of course).
(In the above figure, the Gaussian components are shown in red and the pdf of the mixture models in blue ).
persist_js_state$26c9121e-d294-11ef-18e6-ed8105503adcrunning§runtime vdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$26c9121e-d294-11ef-18e6-ed8105503adcdepends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&˩has_pluto_hook_features¤body \
The Bethe Free Energy and Belief Propagation
We showed that, under mean field assumptions, the FE can be decomposed into a sum of local FE contributions for the nodes ($a$ ) and edges ($i$ ):
$$\begin{align*}
F[q] = \sum_{a=1}^M \underbrace{\sum_{x_a} \left( \prod_{j\in N(a)} q_j(x_j)\cdot \left(-\log p_a(x_a)\right) \right) }_{\text{node energy }U[p_a]} - \sum_{i=1}^N \underbrace{\sum_{x_i} q_i(x_i) \log \frac{1}{q_i(x_i)}}_{\text{edge entropy }H[q_i]}
\end{align*}$$
The mean field assumption is very strong and may lead to large inference costs ($\mathrm{KL}(q(x),p(x|\text{data}))$ ). A more relaxed assumption is to allow joint posterior beliefs over the variables that attach to a node. This idea is expressed by the Bethe Free Energy:
$$\begin{align*}
F_B[q] = \sum_{a=1}^M \left( \sum_{x_a} q_a(x_a) \log \frac{q_a(x_a)}{p_a(x_a)} \right) - \sum_{i=1}^N (d_i - 1) \sum_{x_i} q_i(x_i) \log {q_i(x_i)}
\end{align*}$$
where $q_a(x_a)$ is the posterior joint belief over the variables $x_a$ (i.e., the set of variables that attach to node $a$ ), $q_i(x_i)$ is the posterior marginal belief over the variable $x_i$ and $d_i$ is the number of factor nodes that link to edge $i$ . Moreover, $q_a(x_a)$ and $q_i(x_i)$ are constrained to obey the following equalities:
$$\begin{align*}
\sum_{x_a \backslash x_i} q_a(x_a) &= q_i(x_i), ~~~ \forall i, \forall a \\
\sum_{x_i} q_i(x_i) &= 1, ~~~ \forall i \\
\sum_{x_a} q_a(x_a) &= 1, ~~~ \forall a \\
\end{align*}$$
We form the Lagrangian by augmenting the Bethe Free Energy functional with the constraints:
$$\begin{align*}
L[q] = F_B[q] + \sum_i\sum_{a \in N(i)} \lambda_{ai}(x_i) \left(q_i(x_i) - \sum_{x_a\backslash x_i} q(x_a) \right) + \sum_{i} \gamma_i \left( \sum_{x_i}q_i(x_i) - 1\right) + \sum_{a}\gamma_a \left( \sum_{x_a}q_a(x_a) -1\right)
\end{align*}$$
The stationary solutions for this Lagrangian are given by
$$\begin{align*}
q_a(x_a) &= f_a(x_a) \exp\left(\gamma_a -1 + \sum_{i \in N(a)} \lambda_{ai}(x_i)\right) \\
q_i(x_i) &= \exp\left(1- \gamma_i + \sum_{a \in N(i)} \lambda_{ai}(x_i)\right) ^{\frac{1}{d_i - 1}}
\end{align*}$$
where $N(i)$ denotes the factor nodes that have $x_i$ in their arguments and $N(a)$ denotes the set of variables in the argument of $f_a$ .
Stationary solutions are functions of Lagrange multipliers. This means that Lagrange multipliers need to be determined. Lagrange multipliers can be determined by plugging the stationary solutions back into the constraint specification and solving for the multipliers which ensure that the constraint is satisfied. The first constraint we consider is normalization, which yields the following identification:
$$\begin{align*}
\gamma_a &= 1 - \log \Bigg(\sum_{x_a}f_a(x_a)\exp\left(\sum_{i \in N(a)}\lambda_{ai}(x_i)\right)\Bigg)\\
\gamma_i &= 1 + (d_i-1) \log\Bigg(\sum_{x_i}\exp\left( \frac{1}{d_i-1}\sum_{a \in N(i)} \lambda_{ai}(x_i)\right)\Bigg).
\end{align*}$$
The functional form of the Lagrange multipliers that corresponds to the normalization constraint enforces us to obtain the Lagrange multipliers that correspond to the marginalization constraint. To do so we solve for
$$\begin{align*} \sum_{x_a \backslash x_i} f_a(x_a) \exp\left(\sum_{i \in N(a)} \lambda_{ai}(x_i)\right) &= \exp\left(\sum_{a \in N(i)} \lambda_{ai}(x_i)\right) ^{\frac{1}{d_i - 1}} \exp\left(\lambda_{ai}(x_i)\right)\sum_{x_a \backslash x_i} f_a(x_a) \exp\Bigg(\sum_{\substack{{j \in N(a)} j \neq i}}\lambda_{aj}(x_j)\Bigg) \\
&= \exp\left(\sum_{a \in N(i)} \lambda_{ai}(x_i)\right) ^{\frac{1}{d_i - 1}} \exp\left(\lambda_{ai}(x_i) + \lambda_{ia}(x_i)\right) \\
&= \exp\left(\sum_{a \in N(i)} \lambda_{ai}(x_i)\right) ^{\frac{1}{d_i - 1}}\, ,
\end{align*}$$
where we defined an auxilary function
$$\begin{align*}
\exp(\lambda_{ia}(x_i)) \triangleq \sum_{x_a \backslash x_i} f_a(x_a) \exp\Bigg(\sum_{\substack{{j \in N(a)} j \neq i}}\lambda_{aj}(x_j)\Bigg) \,.
\end{align*}$$
This definition is valid since it can be inverted by the relation
$$\begin{align*}
\lambda_{ia}(x_i) = \frac{2-d_i}{d_i - 1}\lambda_{ai}(x_i) + \frac{1}{d_i -1}\sum_{\substack{c \in N(i)\\c \neq a}}\lambda_{ci}(x_i)
\end{align*}$$
In general it is not possible to solve for the Lagrange multipliers analytically and we resort to iteratively obtaining the solutions. This leads to the Belief Propagation algorithm where the exponentiated Lagrange multipliers (messages) are updated iteratively via
$$\begin{align*}
\mu_{ia}^{(k+1)}(x_i) &= \sum_{x_a \backslash x_i} f_a(x_a) \prod_{\substack{{j \in N(a)} j \neq i}}\mu^{(k)}_{aj}(x_j) \mu_{ai}^{(k)}(x_i) \\
&= \prod_{\substack{c \in N(i) c \neq a}}\mu^{(k)}_{ic}(x_i)\,,
\end{align*}$$
where $k$ denotes iteration number and the messages are defined as
$$\begin{align*}
\mu_{ia}(x_i) &\triangleq \exp(\lambda_{ia}(x_i))\\
\mu_{ai}(x_i) &\triangleq \exp(\lambda_{ai}(x_i))\,.
\end{align*}$$
For a more complete overview of message passing as Bethe Free Energy minimization, see Senoz et al. (2021) .
persist_js_state$26c67f04-d294-11ef-03a4-838ae255689drunning§runtime Fdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$26c67f04-d294-11ef-03a4-838ae255689ddepends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&ùhas_pluto_hook_features¤body
The Variational Free Energy Functional
We'll start from scratch. Consider a model $p(x,z) = p(x|z) p(z)$ , where $x$ and $z$ are observed and latent variables, respectively. $z$ may include parameters but also observation-dependent latent variables.
The goal of Bayesian inference is to transform the (known) likelihood-times-prior factorization of the full model to a posterior-times-evidence decomposition:
$$ \underbrace{p(x|z) p(z)}_{\text{what we know}} \rightarrow \underbrace{p(z|x) p(x)}_{\text{what we want}} $$
Remember from the Bayesian machine learning lesson that negative log-evidence can be decomposed as "complexity" minus "accuracy" terms (the CA decomposition):
$$ -\log p(x) = \underbrace{ \int p(z|x) \log \frac{p(z|x)}{p(z)} \mathrm{d}z }_{\text{complexity}} - \underbrace{\int p(z|x) \log p(x|z) \mathrm{d}z}_{\text{accuracy}}
$$
The CA decomposition cannot be evaluated because it depends on the posterior $p(z|x)$ , which cannot be evaluated since it is the objective of the inference process.
Let's now introduce a distribution $q(z)$ (called the "variational" or "approximate" posterior distribution) that we will use to approximate the posterior $p(z|x)$ . Since we propose the distribution $q(z)$ ourselves, we will assume that $q(z)$ can be evaluated.
If we substitute $q(z)$ for $p(z|x)$ in the CA decomposition, then we obtain
$$ F[q] \triangleq \underbrace{ \int q(z) \log \frac{q(z)}{p(z)} \mathrm{d}z }_{\text{complexity}} - \underbrace{\int q(z) \log p(x|z) \mathrm{d}z}_{\text{accuracy}}
$$
This expression is called the Variational Free Energy (VFE), represented by the symbol $F$ . We treat $F$ as a function of the posterior $q(z)$ . Technically, a function of a function is called a functional, and we write square brackets (e.g., $F[q]$ ) to differentiate functionals from functions (e.g., $q(z)$ ).
Note that all factors in the CA decomposition of VFE (i.e., $q(z)$ , $p(z)$ , and $p(x|z)$ ) can be evaluated as a function of $z$ (and $x$ is observed), and therefore the VFE can be evaluated. This is important: log-evidence $\log p(x)$ cannot be evaluated, but $F[q]$ can be evaluated!
persist_js_state$4ee377c2-a126-4c40-8053-517d40c5ef9drunning§runtimeξ\depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$4ee377c2-a126-4c40-8053-517d40c5ef9ddepends_on_disabled_cells¦queued¤logsoutputmime!application/vnd.pluto.tree+objectrootassigneevfem_resultlast_run_timestampAڹ&has_pluto_hook_features¤bodyobjectidf34ebda5250321fctypeNamedTupleelementsclusters_vb6×120 Matrix{Distribution}:
FullNormal(
dim: 2
μ: [2.33575, 69.4662]
Σ: [0.333333 0.0; 0.0 33.3333]
)
… FullNormal(
dim: 2
μ: [1.98036, 54.6905]
Σ: [0.0396588 0.36343; 0.36343 33.5136]
)
FullNormal(
dim: 2
μ: [2.06093, 90.0417]
Σ: [0.333333 0.0; 0.0 33.3333]
)
FullNormal(
dim: 2
μ: [3.80446, 86.3386]
Σ: [0.409737 0.2216; 0.2216 20.2179]
)
FullNormal(
dim: 2
μ: [3.31476, 72.6624]
Σ: [0.333333 0.0; 0.0 33.3333]
)
FullNormal(
dim: 2
μ: [3.95635, 79.2482]
Σ: [0.163805 1.07272; 1.07272 45.0006]
)
FullNormal(
dim: 2
μ: [3.82633, 60.1778]
Σ: [0.333333 0.0; 0.0 33.3333]
)
FullNormal(
dim: 2
μ: [2.70389, 57.3163]
Σ: [0.423094 4.68636; 4.68636 69.6156]
)
FullNormal(
dim: 2
μ: [4.51578, 72.9886]
Σ: [0.333333 0.0; 0.0 33.3333]
)
FullNormal(
dim: 2
μ: [4.42087, 80.5662]
Σ: [0.101807 0.45467; 0.45467 28.6922]
)
FullNormal(
dim: 2
μ: [3.57291, 62.0741]
Σ: [0.333333 0.0; 0.0 33.3333]
)
… FullNormal(
dim: 2
μ: [2.47015, 55.678]
Σ: [0.265131 2.95173; 2.95173 55.8681]
)
text/plainR 6×272×120 Array{Float64, 3}:
[:, :, 1] =
0.166667 0.166667 0.166667 0.166667 … 0.166667 0.166667 0.166667 0.166667
0.166667 0.166667 0.166667 0.166667 0.166667 0.166667 0.166667 0.166667
0.166667 0.166667 0.166667 0.166667 0.166667 0.166667 0.166667 0.166667
0.166667 0.166667 0.166667 0.166667 0.166667 0.166667 0.166667 0.166667
0.166667 0.166667 0.166667 0.166667 0.166667 0.166667 0.166667 0.166667
0.166667 0.166667 0.166667 0.166667 … 0.166667 0.166667 0.166667 0.166667
[:, :, 2] =
0.0334305 0.791507 0.117296 … 0.000181688 0.35124 0.000683434
0.00660918 1.37119e-7 0.00131983 0.0162787 4.30092e-10 3.01916e-6
0.696345 0.00758745 0.690529 0.119728 0.00164355 0.112473
0.00654938 0.0525124 0.028048 6.40983e-5 0.23365 0.0260373
0.237536 3.09111e-6 0.085711 0.863556 6.58269e-7 0.830583
0.0195301 0.14839 0.0770954 … 0.000192033 0.413466 0.0302208
[:, :, 3] =
4.79769e-5 0.833547 0.00260193 … 2.86942e-11 0.507932 2.60733e-11
0.00240698 1.74887e-13 0.000398605 0.00348038 3.81342e-18 4.02216e-5
0.759349 6.16739e-5 0.888971 0.141414 5.20313e-6 0.0813542
0.00269402 0.0444204 0.0204209 3.27016e-5 0.154463 0.00708219
0.227581 2.50573e-16 0.0286477 0.854988 1.88366e-18 0.905215
0.00792119 0.121971 0.0589594 … 8.45028e-5 0.3376 0.0063085
;;; …
[:, :, 118] =
5.57557e-15 0.971682 5.87457e-10 … 7.95452e-34 0.879584 2.47811e-34
0.000747297 1.70009e-16 0.000140195 0.000987611 2.9377e-22 1.15159e-5
0.795695 4.55769e-8 0.894567 0.210516 4.33408e-9 0.0789169
0.00136513 0.00511513 0.0183056 1.91638e-5 0.0342483 0.00131846
0.199161 5.66805e-17 0.040229 0.78847 5.56927e-19 0.919692
0.00303156 0.0232029 0.0467582 … 7.22977e-6 0.0861679 6.14603e-5
[:, :, 119] =
5.52427e-15 0.971747 5.83772e-10 … 7.79499e-34 0.879796 2.43091e-34
0.000747632 1.66959e-16 0.000140074 0.000987866 2.8557e-22 1.14753e-5
0.795849 4.53433e-8 0.894583 0.210577 4.31487e-9 0.0788944
0.0013654 0.00509487 0.018306 1.92023e-5 0.0341879 0.00131701
0.199007 5.65938e-17 0.0402015 0.788409 5.56709e-19 0.919716
0.00303171 0.0231584 0.0467695 … 7.20931e-6 0.0860159 6.09146e-5
[:, :, 120] =
5.4742e-15 0.971811 5.80166e-10 … 7.64098e-34 0.880006 2.38533e-34
0.000747971 1.63981e-16 0.000139955 0.000988127 2.77643e-22 1.1435e-5
0.796001 4.51133e-8 0.894599 0.210637 4.29596e-9 0.0788722
0.00136566 0.00507485 0.0183063 1.92405e-5 0.0341281 0.00131557
0.198853 5.65073e-17 0.0401742 0.788349 5.56486e-19 0.91974
0.00303184 0.0231145 0.0467807 … 7.18899e-6 0.0858656 6.03781e-5text/plainpersist_js_state$26c5cfb4-d294-11ef-05bb-59d5e27cf37crunning§runtime depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$26c5cfb4-d294-11ef-05bb-59d5e27cf37cdepends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&lhas_pluto_hook_features¤body<
The Marginal Distribution for the GMM
In the literature, the GMM is often introduced by the marginal distribution for an observed data point $x_n$ , given by
$$\begin{align*}{}
p(x_n) &= \sum_{z_n} p(x_n,z_n) \\
&= \sum_{k=1}^K \pi_k \cdot \mathcal{N}\left( x_n | \mu_k, \Sigma_k \right) \tag{B-9.12}
\end{align*}$$
persist_js_state$a7fb83cb-1f40-4c8a-9fda-2165f91e413erunning§runtime ?depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$a7fb83cb-1f40-4c8a-9fda-2165f91e413edepends_on_disabled_cells¦queued¤logsoutputmimeimage/svg+xmlrootassigneelast_run_timestampAڹ&has_pluto_hook_features¤body 7
persist_js_state$64819124-865e-48b8-a916-2ce08dba0accrunning§runtime Hdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$64819124-865e-48b8-a916-2ce08dba0accdepends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&x)has_pluto_hook_features¤body
Show class colors:
persist_js_state$16d90f11-5933-4145-b219-19774eba25d6running§runtime.WSdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$16d90f11-5933-4145-b219-19774eba25d6depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&۷has_pluto_hook_features¤bodyً
Show synthetic label colors:
persist_js_state$06170e31-e865-4178-8af0-41d82df95d71running§runtime 8depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$06170e31-e865-4178-8af0-41d82df95d71depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ& 2has_pluto_hook_features¤body
🎯 Key concept
Constrained VFE minimization is consistent with the Maximum Entropy Principle , which prescribes how to rationally update beliefs when new information becomes available. In this framework, the updated posterior is the distribution that minimizes VFE (or equivalently, KL divergence to the prior) subject to the imposed constraints.
persist_js_state$f1f7407d-86a1-4f24-b78a-61a411d1f371running§runtime ٸdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$f1f7407d-86a1-4f24-b78a-61a411d1f371depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&| has_pluto_hook_features¤bodyV
Variational Inference
persist_js_state$5af3ff1b-1655-4dd9-a089-91544fc85a0erunning§runtime/depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$5af3ff1b-1655-4dd9-a089-91544fc85a0edepends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&ahas_pluto_hook_features¤body<
iteration_em = 1-E persist_js_state$8555aec9-4e80-49e7-8514-ef4a2236801brunning§runtime36depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$8555aec9-4e80-49e7-8514-ef4a2236801bdepends_on_disabled_cells¦queued¤logsoutputmimetext/plainrootassigneeNlast_run_timestampAڹ&,qmhas_pluto_hook_features¤body272persist_js_state$26c65092-d294-11ef-39cc-1953a725f285running§runtime depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$26c65092-d294-11ef-39cc-1953a725f285depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&dhas_pluto_hook_features¤body
Can we still compute an approximate posterior? In this lesson, we introduce an approximate Bayesian inference method known as Variational Bayes (VB) (also known as Variational Inference ) that can be used for Bayesian inference in models with latent variables. Later in this lesson, we will use VB to do inference in the GMM.
persist_js_state$3948225e-90d8-4b78-ba9b-5f98e228285arunning§runtime-(qdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$3948225e-90d8-4b78-ba9b-5f98e228285adepends_on_disabled_cells¦queued¤logsoutputmimeimage/svg+xmlrootassigneelast_run_timestampAڹ&&ohas_pluto_hook_features¤body4
persist_js_state$95d47a10-f3f8-479b-afe0-21241104b758running§runtimee*depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$95d47a10-f3f8-479b-afe0-21241104b758depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&/has_pluto_hook_features¤body*
Code Example: VFEM for GMM on Old Faithfull data setpersist_js_state$26c62ebe-d294-11ef-0cfb-ef186203e890running§runtime 4depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$26c62ebe-d294-11ef-0cfb-ef186203e890depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&ahas_pluto_hook_features¤body$
Since the class labels $y_{nk} \in \{0,1\}$ were assumed to be given by the data set, maximization of this expression decomposed into a set of simple update rules for the Gaussian and multinomial distributions.
persist_js_state$c4b23b39-e6e4-44ec-b204-0c7d7d5a4026running§runtime depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$c4b23b39-e6e4-44ec-b204-0c7d7d5a4026depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&'has_pluto_hook_features¤body
Since the second term (log-evidence) does not involve $q(z)$ , VFEM over $q$ will bring $q(z)$ closer to the Bayesian posterior $p(z|x)$ .
Since $D_{\text{KL}}[q(z),p(z|x)]\geq 0$ for any $q(z)$ , and $D_{\text{KL}}[q(z),p(z|x)]= 0$ only if $q(z) = p(z|x)$ , the VFE is always an upper bound on (minus) log-evidence, i.e.,
$$F[q] \geq -\log p(x) \,.$$
As a result, global minimization of VFE leads to
$$q^*(z) = \arg\min_q F[q]$$
where
$$\begin{align}
q^*(z) &= p(z|x) \tag{= posterior}\\
F[q^*] &= -\log p(x) \tag{= -log evidence}
\end{align}$$
persist_js_state$26c5b896-d294-11ef-1d8e-0feb99d2d45brunning§runtime depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$26c5b896-d294-11ef-1d8e-0feb99d2d45bdepends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&#has_pluto_hook_features¤body
We associate a one-hot coded hidden class label $z_n$ with each observation $x_n$ :
$$\begin{equation*}
z_{nk} = \begin{cases} 1 & \text{if } x_n \in \mathcal{C}_k \text{ (the $k$-th class)}\\
0 & \text{otherwise} \end{cases}
\end{equation*}$$
persist_js_state$26c623f6-d294-11ef-13c0-19edd43592c0running§runtime 8depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$26c623f6-d294-11ef-13c0-19edd43592c0depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&has_pluto_hook_features¤body
Inference for GMM is Difficult
Indeed, the fact that the observation-dependent class labels are unobserved for the GMM, leads to a problem for processing new data by Bayes rule in a GMM.
Consider a given data set $D = \{x_1,x_2,\ldots,x_N\}$ . We recall here the log-likelihood for the Gaussian-Categorial Model, see the generative classification lesson :
$$\log\, p(D|\theta) = \sum_{n,k} y_{nk} \underbrace{ \log\mathcal{N}(x_n|\mu_k,\Sigma) }_{ \text{Gaussian} } + \underbrace{ \sum_{n,k} y_{nk} \log \pi_k }_{ \text{multinomial} } \,.$$
persist_js_state$55570464-89c8-4d9b-b667-dfa64ac62294running§runtime depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$55570464-89c8-4d9b-b667-dfa64ac62294depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&]has_pluto_hook_features¤body4
Code
persist_js_state$26c59b52-d294-11ef-1eba-d3f235f85eeerunning§runtime -depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$26c59b52-d294-11ef-1eba-d3f235f85eeedepends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&Xhas_pluto_hook_features¤body
Unobserved Classes
Consider again a set of observed data $D=\{x_1,\dotsc,x_N\}$ .
This time, we suspect that there are unobserved class labels that would help explain (or predict) the data, e.g.,
the observed data are the color of living things; the unobserved classes are animals and plants.
observed are wheel sizes; unobserved categories are trucks and personal cars.
observed is an audio signal; unobserved classes include speech, music, traffic noise, etc.
persist_js_state$f17c9d8a-9291-4110-bcf4-c582d23f986brunning§runtime depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$f17c9d8a-9291-4110-bcf4-c582d23f986bdepends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&ʐChas_pluto_hook_features¤body
Variational Message Passing
Let us now consider the local free energy that is associated with edge corresponding to $x_j$ .
persist_js_state$dd1242db-fb20-4732-ac55-a3e021bbd2b7running§runtime ;depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$dd1242db-fb20-4732-ac55-a3e021bbd2b7depends_on_disabled_cells¦queued¤logsoutputmime!application/vnd.pluto.tree+objectrootassigneeconst data_plot_kwargslast_run_timestampAڹ&thas_pluto_hook_features¤bodyobjectidcc4be3e3968086c0typeNamedTupleelementsmarkersize4text/plainmarkerstrokewidth0text/plaincolor:red4text/plainopacity0.8text/plainpersist_js_state$26c73cf0-d294-11ef-297b-354eb9c71f57running§runtime Kdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$26c73cf0-d294-11ef-297b-354eb9c71f57depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ& has_pluto_hook_features¤body(
persist_js_state$2d4adbf6-6de8-4e3a-ad6f-fa8bbfa5999erunning§runtime !depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$2d4adbf6-6de8-4e3a-ad6f-fa8bbfa5999edepends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&Ɂkhas_pluto_hook_features¤body
persist_js_state$26c56fd8-d294-11ef-236d-81deef63f37crunning§runtimeͪdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$26c56fd8-d294-11ef-236d-81deef63f37cdepends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&TKhas_pluto_hook_features¤body
Latent Variable Models and Variational Bayespersist_js_state$666680b2-315a-4d95-8f7f-3ae50018e112running§runtime&depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$666680b2-315a-4d95-8f7f-3ae50018e112depends_on_disabled_cells¦queued¤logsoutputmimetext/plainrootassigneelast_run_timestampAڹ&>has_pluto_hook_features¤bodypersist_js_state$0349720e-5de4-4b39-babd-c0881588f1derunning§runtimedepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$0349720e-5de4-4b39-babd-c0881588f1dedepends_on_disabled_cells¦queued¤logsoutputmimetext/plainrootassigneeXlast_run_timestampAڹ&d{̷has_pluto_hook_features¤body2×272 Matrix{Float64}:
3.6 1.8 3.333 2.283 4.533 2.883 … 4.117 2.15 4.417 1.817 4.467
79.0 54.0 74.0 62.0 85.0 55.0 81.0 46.0 90.0 46.0 74.0persist_js_state$cd7cf2d0-aff7-49dc-bb12-dc73c6433768running§runtime%0ĸdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$cd7cf2d0-aff7-49dc-bb12-dc73c6433768depends_on_disabled_cells¦queued¤logsoutputmimeimage/svg+xmlrootassigneelast_run_timestampAڹ&has_pluto_hook_features¤body i
persist_js_state$b4d965d1-91d8-43f8-84a5-b37ad5c6cafarunning§runtime depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$b4d965d1-91d8-43f8-84a5-b37ad5c6cafadepends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&has_pluto_hook_features¤bodyNote that the inference bound is a Kullback-Leibler (KL) divergence between an (approximate) posterior $q(z)$ and the (perfect) Bayesian posterior $p(z|x)$ . To learn more:
Apparently (see previous slide), there are three contributions to the free energy for $x_j$ :
The local free energy for $x_j$ can be written as (exercise)
$$ F[q_j] \propto \sum_{x_j} q(x_j) \log \frac{q_j(x_j)}{\nu_a(x_j)\cdot \nu_b(x_j)}
$$
where
$$\begin{align*}
\nu_a(x_j) &\propto \exp\left( \mathbb{E}_{q_{k}}\left[ \log p_a(x_a)\right]\right) \\
\nu_b(x_j) &\propto \exp\left( \mathbb{E}_{q_{l}}\left[ \log p_b(x_b)\right]\right)
\end{align*}$$
and $\mathbb{E}_{q_{k}}\left[\cdot\right]$ is an expectation w.r.t. all $q(x_k)$ with $k \in N(a)\setminus {j}$ .
$\nu_a(x_j)$ and $\nu_b(x_j)$ can be locally computed in nodes $a$ and $b$ respectively and can be interpreted as colliding messages over edge $x_j$ .
Local free energy minimization is achieved by setting
$$ q_j(x_j) \propto \nu_a(x_j) \cdot \nu_b(x_j)
$$
Note that message $\nu_a(x_j)$ depends on posterior beliefs over incoming edges ($k$ ) for node $a$ , and in turn, the message from node $a$ towards edge $x_k$ depends on the belief $q_j(x_j)$ . I.o.w., direct mutual dependencies exist between posterior beliefs over edges that attach to the same node.
These considerations lead to the Variational Message Passing procedure, which is an iterative free energy minimization procedure that can be executed completely through locally computable messages.
Procedure VMP, see Dauwels (2007), section 3
Initialize all messages $q$ and $ν$ , e.g., $q(\cdot) \propto 1$ and $\nu(\cdot) \propto 1$ . <br/>
Select an edge $z_k$ in the factor graph of $f(z_1,\ldots,z_m)$ .<br/>
Compute the two messages $\overrightarrow{\nu}(z_k)$ and $\overleftarrow{\nu}(z_k)$ by applying the following generic rule:
$$ \overrightarrow{\nu}(y) \propto \exp\left( \mathbb{E}_{q}\left[ \log g(x_1,\dots,x_n,y)\right] \right) $$
Compute the marginal $q(z_k)$
$$ q(z_k) \propto \overrightarrow{\nu}(z_k) \overleftarrow{\nu}(z_k) $$
and send it to the two nodes connected to the edge $x_k$ .
Iterate 2–4 until convergence.
FE Minimization by the Expectation-Maximization (EM) Algorithm The EM algorithm is a special case of VFE minimization that focuses on Maximum-Likelihood estimation for models with latent variables.
Consider a model
$$p(x,z,\theta)$$
with observations $x = \{x_n\}$ , latent variables $z=\{z_n\}$ and parameters $\theta$ .
We can write the following VFE functional for this model:
$$\begin{align*}
F[q] = \sum_z \sum_\theta q(z,\theta) \log \frac{q(z,\theta)}{p(x,z,\theta)}
\end{align*}$$
The EM algorithm makes the following simplifying assumptions:
The prior for the parameters is uninformative (uniform). This implies that
$$p(x,z,\theta) = p(x,z|\theta) p(\theta) \propto p(x,z|\theta)$$
A factorization constraint
$$q(z,\theta) = q(z) q(\theta)$$
The posterior for the parameters is a delta function:
$$q(\theta) = \delta(\theta - \hat{\theta})$$
Basically, these three assumptions turn VFE minimization into maximum likelihood estimation for the parameters $\theta$ and the VFE simplifies to
$$\begin{align*}
F[q,\theta] = \sum_z q(z) \log \frac{q(z)}{p(x,z|\theta)}
\end{align*}$$
The EM algorithm minimizes this FE by iterating (iteration counter: $i$ ) over
$$\begin{align} \mathcal{L}^{(i)}(\theta) &= \sum_z \overbrace{p(z|x,\theta^{(i-1)})}^{q^{(i)}(z)} \log p(x,z|\theta) \tag{the E-step} \\
\theta^{(i)} &= \arg\max_\theta \mathcal{L}^{(i)}(\theta) \tag{the M-step} \end{align}
$$
These choices are optimal for the given FE functional. In order to see this, consider the two decompositions
$$\begin{align*}
F[q,\theta] &= \underbrace{-\sum_z q(z) \log p(x,z|\theta)}_{\text{energy}} - \underbrace{\sum_z q(z) \log \frac{1}{q(z)}}_{\text{entropy}} \qquad &&\text{(EE)}\\
&= \underbrace{\sum_z q(z) \log \frac{q(z)}{p(z|x,\theta)}}_{\text{divergence}} - \underbrace{\log p(x|\theta)}_{\text{log-likelihood}} \qquad &&\text{(DE)}
\end{align*}$$
The DE decomposition shows that the FE is minimized for the choice $q(z) := p(z|x,\theta)$ . Also, for this choice, the FE equals the (negative) log-evidence (, which is this case simplifies to the log-likelihood).
The EE decomposition shows that the FE is minimized wrt $\theta$ by minimizing the energy term. The energy term is computed in the E-step and optimized in the M-step.
In order to execute the EM algorithm, it is assumed that we can analytically execute the E- and M-steps. For a large set of models (including models whose distributions belong to the exponential family of distributions), this is indeed the case and hence the large popularity of the EM algorithm.
The EM algorihm imposes rather severe assumptions on the FE (basically approximating Bayesian inference by maximum likelihood estimation). Over the past few years, the rise of Probabilistic Programming languages has dramatically increased the range of models for which the parameters can by estimated autmatically by (approximate) Bayesian inference, so the popularity of EM is slowly waning. (More on this in the Probabilistic Programming lessons).
Bishop (2006) works out EM for the GMM in section 9.2.
persist_js_state$26c7696e-d294-11ef-25f2-dbc0946c0858running§runtime depends_on_skipped_cellsµpublished_object_keyserrored§cell_id$26c7696e-d294-11ef-25f2-dbc0946c0858depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&&has_pluto_hook_features¤bodyٷBelow we exemplify training of a Gaussian Mixture Model on the Old Faithful data set by VFE minimization, with the constraints as specified above.
Click for the solution
Let's start with a prior estimate $\theta^{(i-1)}$ and we want to minimize the free energy functional wrt $q$ . This leads to
$$\begin{align*}
q^{(i)}(z) &= \arg\min_q F[q](\theta^{(i-1)}) \\
&= \arg\min_q \sum_z q(z) \log \frac{q(z)}{p(D,z|\theta^{(i-1)})} \\
&= \arg\min_q \sum_z q(z) \log \frac{q(z)}{p(z|D,\theta^{(i-1)}) \cdot p(D|\theta^{(i-1)})} \\
&= p(z|D,\theta^{(i-1)})
\end{align*}$$
Next, we use $q^{(i)}(z)=p(z|D,\theta^{(i-1)})$ and minimize the free energy w.r.t. $\theta$ , leading to
$$\begin{align*}
\theta^{(i)} &= \arg\min_\theta F[q^{(i)}](\theta) \\
&= \arg\min_\theta \sum_z p(z|D,\theta^{(i-1)}) \log \frac{p(z|D,\theta^{(i-1)})}{p(D,z|\theta)} \\
&= \arg\max_\theta \sum_z \underbrace{p(z|D,\theta^{(i-1)})}_{q^{(i)}(z)} \log p(D,z|\theta)
\end{align*}$$
persist_js_state$87d94630-c90b-4379-91bd-88641ee7b508running§runtime Odepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$87d94630-c90b-4379-91bd-88641ee7b508depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&$phas_pluto_hook_features¤bodyImplementation
Reading this lecture online? Click "View code" in the top right to read the implementation of this visualisation.
Exercises iteration_vfem = 2 Problem
You are asked to build a density model for a data set (Old Faithful , Bishop pg. 681) that clearly is not distributed as a single Gaussian:
Classification problems with unobserved classes are called Clustering problems. In clustering problems, the learning algorithm needs to discover the underlying classes from the observed data .
The Expectation-Maximization (EM) algorithm (**)
Consider a model $p(x,z|\theta)$ , where $D=\{x_1,x_2,\ldots,x_N\}$ is observed, $z$ are unobserved variables, and $\theta$ are parameters. The Expectation-Maximization (EM) algorithm estimates the parameters by iterating over the following two equations ($i$ is the iteration index):
$$\begin{align*}
q^{(i)}(z) &= p(z|D,\theta^{(i-1)}) \\
\theta^{(i)} &= \arg\max_\theta \sum_z q^{(i)}(z) \cdot \log p(D,z|\theta)
\end{align*}$$
Prove that this algorithm minimizes the Free Energy functional
$$\begin{align*}
F[q](\theta) = \sum_z q(z) \log \frac{q(z)}{p(D,z|\theta)}
\end{align*}$$
The generated figure resembles Figure 10.6 in Bishop. The plots show VFEM results for a GMM of $K = 6$ Gaussians applied to the Old Faithful data set. The ellipses denote the two standard-deviation density contours for each of the components, and the color coding of the data points reflects the "soft" class label assignments.
We'll perform clustering on the data set from the challenge by fitting a GMM consisting of two Gaussians using the EM algorithm.
Optional Slides
persist_js_state$26c64174-d294-11ef-2bbc-ab1a84532311running§runtime jdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$26c64174-d294-11ef-2bbc-ab1a84532311depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&Ihas_pluto_hook_features¤bodyThere is no known conjugate prior for the latent variables in the GMM likelihood. Therefore, Bayes rule does not yield a closed-form expression for the posterior over the latent variables:
$$ \underbrace{p(\{z_{nk}\},\{\mu_k,\Sigma_k,\pi_k\} | D)}_{\text{posterior (no analytical solution)}} \propto \underbrace{p(D\,|\,\{z_{nk}\},\{\mu_k,\Sigma_k,\pi_k\})}_{\text{likelihood}} \cdot \underbrace{p( \{z_{nk}\},\{\mu_k,\Sigma_k,\pi_k\} )}_{\text{prior (no known conjugate)}} $$
Click for the solution
We see here an example of "precision-weighted means add" when two sources of information are fused, just like precision-weighted means add when two Gaussians are multiplied, eg a prior and likelihood. In this case, the prior is $m_0$ and the likelihood estimate is $\bar{x}$ . $\beta_0$ can be interpreted as the number of pseudo-observations in the prior.
persist_js_state$26c6f63c-d294-11ef-1090-e9238dd6ad3frunning§runtime ydepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$26c6f63c-d294-11ef-1090-e9238dd6ad3fdepends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&ݷhas_pluto_hook_features¤bodyConstrained VFE Minimization
It is common to add simplifying constraints to an optimization problem to make a difficult optimization task tractible. This is also standard practice when approximating Bayesian inference by FE minimization.
There are three important cases of adding constraints to the VFE functional that often alleviate the VFE minimization task:
form constraints on $q(z)$
factorization constraints on $q(z)$
other (ad hoc) constraints
We will shortly discuss these simplifications below.
For almost every practical setting, we constrain the posterior $q(z)$ to belong to a specific parameterized family of probability distributions, e.g.,
$$q(z) = \mathcal{N}\left( z | \mu, \Sigma \right)\,.$$
In this case, the functional minimization problem for $F[q]$ simplifies to the minimization of an ordinary function
$$F(\mu,\Sigma) = \int \mathcal{N}\left( z | \mu, \Sigma \right) \log \frac{\mathcal{N}\left( z | \mu, \Sigma \right)}{p(x,z)}\mathrm{d}z$$
with respect to the parameters $\mu$ and $\Sigma$ .
We can often use standard gradient-based optimization methods to minimize the function $F(\mu,\Sigma)\,.$
persist_js_state$26c82f16-d294-11ef-0fe1-07326b56282frunning§runtime "Edepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$26c82f16-d294-11ef-0fe1-07326b56282fdepends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&uM has_pluto_hook_features¤bodyVFE Minimization with Mean-field Factorization Constraints: the CAVI Approach
Let's work out VFE minimization with additional mean-field constraints (=full factorization) constraints:
$$q(z) = \prod_{j=1}^m q_j(z_j)\,.$$
In other words, the posteriors for $z_j$ are all considered independent. This is a strong constraint but often leads to good solutions.
Given the mean-field constraints, it is possible to derive the following expression for the optimal solutions $q_j^*(z_j)$ , for $j=1,\ldots,m$ :
$$\begin{align}
\log q_j^*(z_j) &\propto \mathrm{E}_{q_{-j}^*}\left[ \log p(x,z) \right] \\
&= \underbrace{\sum_{z_{-j}} q_{-j}^*(z_{-j}) \underbrace{\log p(x,z)}_{\text{"field"}}}_{\text{"mean field"}}
\end{align} $$
where we defined $q_{-j}^*(z_{-j}) \triangleq q_1^*(z_1)q_2^*(z_2)\cdots q_{j-1}^*(z_{j-1})q_{j+1}^*(z_{j+1})\cdots q_m^*(z_m)$ .
Proof (from Blei, 2017 ): We first rewrite the FE as a function of $q_j(z_j)$ only:
$$ F[q_j] = \mathbb{E}_{q_{j}}\left[ \mathbb{E}_{q_{-j}}\left[ \log p(x,z_j,z_{-j})\right]\right] - \mathbb{E}_{q_j}\left[ \log q_j(z_j)\right] + \mathtt{const.}\,,$$
where the constant holds all terms that do not depend on $z_j$ . This expression can be written as
$$ F[q_j] = \sum_{z_j} q_j(z_j) \log \frac{q_j(z_j)}{\exp\left( \mathbb{E}_{q_{-j}}\left[ \log p(x,z_j,z_{-j})\right]\right)}$$
which is a KL-divergence that is minimized by Eq. B-10.9. (end proof)
This is not yet a full solution to the FE minimization task since the solution $q_j^*(z_j)$ depends on expectations that involve other solutions $q_{i\neq j}^*(z_{i \neq j})$ , and each of these other solutions $q_{i\neq j}^*(z_{i \neq j})$ depends on an expection that involves $q_j^*(z_j)$ .
In practice, we solve this chicken-and-egg problem by an iterative approach: we first initialize all $q_j(z_j)$ (for $j=1,\ldots,m$ ) to an appropriate initial distribution and then cycle through the factors in turn by solving eq.B-10.9 and update $q_{-j}^*(z_{-j})$ with the latest estimates. (See Blei, 2017 , Algorithm 1, p864).
This algorithm for approximating Bayesian inference is known Coordinate Ascent Variational Inference (CAVI).
🎯 Key concept
For a model $p(x,z) = p(x|z) p(x)$ , the variational free energy functional
$$ F[q] = \underbrace{ \int q(z) \log \frac{q(z)}{p(z)} \mathrm{d}z }_{\text{complexity}} - \underbrace{\int q(z) \log p(x|z) \mathrm{d}z}_{\text{accuracy}}$$
is considered a functional of a variational posterior $q(z)$ over the latent variables $z$ . Minimizing $F[q]$ accomplishes approximate Bayesian inference by
$$\begin{align*}
\hat{q}(z) &\approx p(z|x) \\
F[\hat{q}] &\approx -\log p(x)
\end{align*}$$
rec(),
]
# Execute the algorithm: iteratively update parameters and responsibilities
for i in 1:10
updateResponsibilities!(γ, X, clusters, π_hat)
push!(em_result, (i, :E) => rec())
updateParameters!(clusters, π_hat, X, γ)
push!(em_result, (i, :M) => rec())
end
Text("em_result")
end$26c5a93a-d294-11ef-23a1-cbcf0c370fc9metadatadisabled©show_logsîskip_as_script§cell_id$26c5a93a-d294-11ef-23a1-cbcf0c370fc9code_foldedäcodeٸmd"""
## The Gaussian Mixture Model
The spread of the data in the Old Faithful data set looks like it could be modeled by two Gaussians. Let's develop a model for this data set.
"""$3e897a59-e7b5-492c-8a8a-724248513a72metadatadisabled©show_logsîskip_as_script§cell_id$3e897a59-e7b5-492c-8a8a-724248513a72code_foldedäcodemd"""
##### Model Specification
We consider a Gaussian Mixture Model, specified by
```math
\begin{align*}
p(x,z|\theta) &= p(x|z,\mu,\Lambda)p(z|\pi) \\
&= \prod_{n=1}^N \prod_{k=1}^K \mathcal{N}\left( x_n | \mu_k, \Lambda_k^{-1}\right)^{z_{nk}} \cdot \prod_{n=1}^N \prod_{k=1}^K \pi_k^{z_{nk}} \\
&= \prod_{n=1}^N \prod_{k=1}^K \left(\pi_k \cdot \mathcal{N}\left( x_n | \mu_k, \Lambda_k^{-1}\right)\right)^{z_{nk}} \tag{B-10.37,38}
\end{align*}
```
Let us introduce some priors for the parameters ``\pi``, ``\mu``, and ``\Lambda``. We factorize the prior and choose conjugate distributions by
```math
p(\pi,\mu,\Lambda) = p(\pi) p(\mu|\Lambda) p(\Lambda)
```
with
```math
\begin{align}
p(\pi) &= \mathrm{Dir}(\pi|\alpha_0) = C(\alpha_0) \prod_k \pi_k^{\alpha_0-1} \tag{B-10.39} \\
p(\mu|\Lambda) &= \prod_k \mathcal{N}\left(\mu_k | m_0, \left( \beta_0 \Lambda_k\right)^{-1} \right) \tag{B-10.40} \\
p(\Lambda) &= \prod_k \mathcal{W}\left( \Lambda_k | W_0, \nu_0 \right) \tag{B-10.40}
\end{align}
```
where ``\mathcal{W}\left( \cdot \right)`` is a [Wishart distribution](https://en.wikipedia.org/wiki/Wishart_distribution) (i.e., a multi-dimensional Gamma distribution).
The full generative model is now specified by
```math
p(x,z,\pi,\mu,\Lambda) = \underbrace{p(x|z,\mu,\Lambda) p(z|\pi)}_{\text{B-10.37-38}} \underbrace{p(\pi) p(\mu|\Lambda) p(\Lambda)}_{\text{B-10.39-40}} \tag{B-10.41}
```
with hyperparameters ``\{ \alpha_0, m_0, \beta_0, W_0, \nu_0\}``.
"""$26c8068a-d294-11ef-3983-a1be55128b3fmetadatadisabled©show_logsîskip_as_script§cell_id$26c8068a-d294-11ef-3983-a1be55128b3fcode_foldedäcodemd"""
## Variational Inference in Practice
For most realistic models of complex real-world problems, Bayes rule is not tractable in closed form. As a result, the use of approximate variational Bayesian inference has seen rapid growth in practical applications.
Toolboxes such as [RxInfer](http://rxinfer.com) enable users to define sophisticated probabilistic models and automate the inference process via constrained VFE minimization. Remarkably, specifying even complex models typically requires no more than a single page of code.
In contrast to traditional algorithm design, where solving a problem might require implementing a custom solution in, say, ``40`` pages of code, automated inference in a probabilistic model offers a radically more efficient and modular approach. This shift has the potential to fundamentally change how we design and deploy information processing systems in the future.
"""$ae7ed1fc-fc36-4327-be55-a142477ca0admetadatadisabled©show_logsîskip_as_script§cell_id$ae7ed1fc-fc36-4327-be55-a142477ca0adcode_foldedäcode hide_proof(
md"""
```math
\begin{align*}
F[q] &= \underbrace{ \int q(z) \log \frac{q(z)}{p(z)} \mathrm{d}z }_{\text{complexity}} - \underbrace{\int q(z) \log p(x|z) \mathrm{d}z}_{\text{accuracy}} \\
&= \int q(z) \log \frac{q(z)}{ p(x|z) p(z) }\mathrm{d}z \\
&= \int q(z) \log \frac{q(z)}{ p(z|x) p(x)}\mathrm{d}z \quad \text{( since } p(x|z) p(z) = p(z|x) p(x)\text{ )} \\
&= \underbrace{\int q(z) \log \frac{q(z)}{p(z|x)}\mathrm{d}z}_{\text{inference bound}\geq 0} - \underbrace{\log p(x)}_{\text{log-evidence}}
\end{align*}
```
""")$86c33a7c-135a-461f-a17e-b50bca418e13metadatadisabled©show_logsîskip_as_script§cell_id$86c33a7c-135a-461f-a17e-b50bca418e13code_foldedäcodefunction sufficientStatistics(X,r,k::Int) #function to compute sufficient statistics
N_k = sum(r[k,:])
hat_x_k = sum([r[k,n]*X[:,n] for n in 1:N]) ./ N_k
S_k = sum([r[k,n]*(X[:,n]-hat_x_k)*(X[:,n]-hat_x_k)' for n in 1:N]) ./ N_k
return N_k, hat_x_k, S_k
end$26c8b682-d294-11ef-1331-2bcf8baec73fmetadatadisabled©show_logsîskip_as_script§cell_id$26c8b682-d294-11ef-1331-2bcf8baec73fcode_foldedäcodeUmd"""
## Message Passing for Free Energy Minimization
The Sum-Product (SP) update rule implements perfect Bayesian inference.
Sometimes, the SP update rule is not analytically solvable.
Fortunately, for many well-known Bayesian approximation methods, a message passing update rule can be created, e.g. [Variational Message Passing](https://en.wikipedia.org/wiki/Variational_message_passing) (VMP) for variational inference.
In general, all of these message passing algorithms can be interpreted as minimization of a constrained free energy (e.g., see [Senoz et al. (2021)](https://research.tue.nl/nl/publications/variational-message-passing-and-local-constraint-manipulation-in-), and hence these message passing schemes comply with [Caticha's Method of Maximum Relative Entropy](https://arxiv.org/abs/1011.0723), which, as discussed in the [variational Bayes lesson](https://bmlip.github.io/course/lectures/Latent%20Variable%20Models%20and%20VB.html) is the proper way for updating beliefs.
Different message passing updates rules can be combined to get a hybrid inference method in one model.
"""$f42a1a65-20ce-452f-9974-bc8146943574metadatadisabled©show_logsîskip_as_script§cell_id$f42a1a65-20ce-452f-9974-bc8146943574code_foldedäcode8md"""
# Theoretical Underpinning of VFE Minimization
"""$edb179df-5cff-4e7b-8645-6da4818dceeemetadatadisabled©show_logsîskip_as_script§cell_id$edb179df-5cff-4e7b-8645-6da4818dceeecode_foldedäcodemd"""
##### 3. Other constraints, e.g., the Expectation-Maximization (EM) algorithm
Beyond form and factorization constraints on ``q(z)``, several ad hoc algorithms have been developed to simplify VFE minimization for particular models.
In particular, the [Expectation-Maximization (EM) algorithm](https://en.wikipedia.org/wiki/Expectation%E2%80%93maximization_algorithm) is a famous special case of constrained VFE minimization. The EM algorithm places some constraints on both the posterior ``q(z)`` and the prior ``p(z)`` (see the [OPTIONAL SLIDE](#EM-Algorithm) for more info) that essentially reduce VFE minimization to maximum likelihood estimation.
"""$06512595-bdb7-4adf-88ae-62af20210891metadatadisabled©show_logsîskip_as_script§cell_id$06512595-bdb7-4adf-88ae-62af20210891code_foldedäcodeKchallenge_solution("Modeling of the Old Faithful Data Set"; header_level=1)$3deadfd0-9fbb-476a-a7de-5dd694e55a65metadatadisabled©show_logsîskip_as_script§cell_id$3deadfd0-9fbb-476a-a7de-5dd694e55a65code_foldedäcodeemd"""
Eq. B-9.12 reveals the link to the name Gaussian *mixture model*. The priors ``\pi_k`` for the ``k``-th class are also called **mixture coefficients**.
Be aware that Eq. B-9.12 is not the generative model for the GMM! The generative model is the joint distribution ``p(x,z,\pi,\mu,\Sigma)`` over all variables, including the latent variables.
"""$3654551d-5d08-4bb0-8a0d-c7d42225bc69metadatadisabled©show_logsîskip_as_script§cell_id$3654551d-5d08-4bb0-8a0d-c7d42225bc69code_foldedäcodemd"""
##### 2. Factorization constraints
In addition to form constraints, it is also common to constrain the posterior ``q(z)`` by a specific factorization. For instance, in the *mean-field factorization* constraints, we constrain the variational posterior to factorize fully into a set of independent factors, i.e.,
```math
q(z) = \prod_{j=1}^m q_j(z_j)\,, \tag{B-10.5}
```
Variational inference with mean-field factorization has been worked out in detail as the **Coordinate Ascent Variational Inference** (CAVI) algorithm. See the [Optional Slide on CAVI](#CAVI) for details.
Mean-field factorization is just an example of various _factorization constraints_ that have been successfully applied to VFEM.
"""$670f83cb-66d6-40ab-9c42-50e022a9bd4cmetadatadisabled©show_logsîskip_as_script§cell_id$670f83cb-66d6-40ab-9c42-50e022a9bd4ccode_foldedäcodeَmd"""
In other words, global minimization of ``F[q]`` recovers Bayes rule _exactly_, namely, both the posterior and evidence are revealed.
"""$c563de07-0036-4dc4-87ee-e8363e9b6b1bmetadatadisabled©show_logsîskip_as_script§cell_id$c563de07-0036-4dc4-87ee-e8363e9b6b1bcode_foldedäcode٨navigate_prev_next(
"https://bmlip.github.io/course/lectures/Discriminative%20Classification.html",
"https://bmlip.github.io/course/lectures/Dynamic%20Models.html",
)$90fbe618-fc81-480b-b685-69cd97e5b8edmetadatadisabled©show_logsîskip_as_script§cell_id$90fbe618-fc81-480b-b685-69cd97e5b8edcode_foldedäcodeّResource("https://github.com/bmlip/course/blob/v2/assets/figures/VMP-two-nodes.png?raw=true", :style => "background: white; border-radius: 1em;")$f8c8013a-3e87-4d01-a3ae-86b39cf1f002metadatadisabled©show_logsîskip_as_script§cell_id$f8c8013a-3e87-4d01-a3ae-86b39cf1f002code_foldedäcode&md"""
# The Gaussian Mixture Model
"""$4e0c025d-fa39-462d-8e7d-e66d220e9595metadatadisabled©show_logsîskip_as_script§cell_id$4e0c025d-fa39-462d-8e7d-e66d220e9595code_folded¤codemax_iterations = 120;$ce7d086b-ff20-4da1-a4e8-52b5b7dc9e2bmetadatadisabled©show_logsîskip_as_script§cell_id$ce7d086b-ff20-4da1-a4e8-52b5b7dc9e2bcode_foldedäcodePlutoUI.TableOfContents()$26c58298-d294-11ef-2a53-2b42b48e0725metadatadisabled©show_logsîskip_as_script§cell_id$26c58298-d294-11ef-2a53-2b42b48e0725code_foldedäcode{md"""
## Preliminaries
##### Goal
* Introduction to latent variable models and variational inference by Free energy minimization
##### Materials
* Mandatory
* These lecture notes
* Optional
* Bishop (2016), [PRML book](https://www.microsoft.com/en-us/research/wp-content/uploads/2006/01/Bishop-Pattern-Recognition-and-Machine-Learning-2006.pdf), pp. 461-486 (sections 10.1, 10.2 and 10.3)
* Ariel Caticha (2010), [Entropic Inference](https://arxiv.org/abs/1011.0723)
* tutorial on entropic inference, which is a generalization to Bayes rule and provides a foundation for variational inference.
* references $(HTML("red ") and the pdf of the mixture models in $(html"blue ")).
"""$26c9121e-d294-11ef-18e6-ed8105503adcmetadatadisabled©show_logsîskip_as_script§cell_id$26c9121e-d294-11ef-18e6-ed8105503adccode_foldedäcodekmd"""
## The Bethe Free Energy and Belief Propagation
We showed that, under mean field assumptions, the FE can be decomposed into a sum of local FE contributions for the nodes (``a``) and edges (``i``):
```math
\begin{align*}
F[q] = \sum_{a=1}^M \underbrace{\sum_{x_a} \left( \prod_{j\in N(a)} q_j(x_j)\cdot \left(-\log p_a(x_a)\right) \right) }_{\text{node energy }U[p_a]} - \sum_{i=1}^N \underbrace{\sum_{x_i} q_i(x_i) \log \frac{1}{q_i(x_i)}}_{\text{edge entropy }H[q_i]}
\end{align*}
```
The mean field assumption is very strong and may lead to large inference costs (``\mathrm{KL}(q(x),p(x|\text{data}))``). A more relaxed assumption is to allow joint posterior beliefs over the variables that attach to a node. This idea is expressed by the Bethe Free Energy:
```math
\begin{align*}
F_B[q] = \sum_{a=1}^M \left( \sum_{x_a} q_a(x_a) \log \frac{q_a(x_a)}{p_a(x_a)} \right) - \sum_{i=1}^N (d_i - 1) \sum_{x_i} q_i(x_i) \log {q_i(x_i)}
\end{align*}
```
where ``q_a(x_a)`` is the posterior joint belief over the variables ``x_a`` (i.e., the set of variables that attach to node ``a``), ``q_i(x_i)`` is the posterior marginal belief over the variable ``x_i`` and ``d_i`` is the number of factor nodes that link to edge ``i``. Moreover, ``q_a(x_a)`` and ``q_i(x_i)`` are constrained to obey the following equalities:
```math
\begin{align*}
\sum_{x_a \backslash x_i} q_a(x_a) &= q_i(x_i), ~~~ \forall i, \forall a \\
\sum_{x_i} q_i(x_i) &= 1, ~~~ \forall i \\
\sum_{x_a} q_a(x_a) &= 1, ~~~ \forall a \\
\end{align*}
```
We form the Lagrangian by augmenting the Bethe Free Energy functional with the constraints:
```math
\begin{align*}
L[q] = F_B[q] + \sum_i\sum_{a \in N(i)} \lambda_{ai}(x_i) \left(q_i(x_i) - \sum_{x_a\backslash x_i} q(x_a) \right) + \sum_{i} \gamma_i \left( \sum_{x_i}q_i(x_i) - 1\right) + \sum_{a}\gamma_a \left( \sum_{x_a}q_a(x_a) -1\right)
\end{align*}
```
The stationary solutions for this Lagrangian are given by
```math
\begin{align*}
q_a(x_a) &= f_a(x_a) \exp\left(\gamma_a -1 + \sum_{i \in N(a)} \lambda_{ai}(x_i)\right) \\
q_i(x_i) &= \exp\left(1- \gamma_i + \sum_{a \in N(i)} \lambda_{ai}(x_i)\right) ^{\frac{1}{d_i - 1}}
\end{align*}
```
where ``N(i)`` denotes the factor nodes that have ``x_i`` in their arguments and ``N(a)`` denotes the set of variables in the argument of ``f_a``.
Stationary solutions are functions of Lagrange multipliers. This means that Lagrange multipliers need to be determined. Lagrange multipliers can be determined by plugging the stationary solutions back into the constraint specification and solving for the multipliers which ensure that the constraint is satisfied. The first constraint we consider is normalization, which yields the following identification:
```math
\begin{align*}
\gamma_a &= 1 - \log \Bigg(\sum_{x_a}f_a(x_a)\exp\left(\sum_{i \in N(a)}\lambda_{ai}(x_i)\right)\Bigg)\\
\gamma_i &= 1 + (d_i-1) \log\Bigg(\sum_{x_i}\exp\left( \frac{1}{d_i-1}\sum_{a \in N(i)} \lambda_{ai}(x_i)\right)\Bigg).
\end{align*}
```
The functional form of the Lagrange multipliers that corresponds to the normalization constraint enforces us to obtain the Lagrange multipliers that correspond to the marginalization constraint. To do so we solve for
```math
\begin{align*} \sum_{x_a \backslash x_i} f_a(x_a) \exp\left(\sum_{i \in N(a)} \lambda_{ai}(x_i)\right) &= \exp\left(\sum_{a \in N(i)} \lambda_{ai}(x_i)\right) ^{\frac{1}{d_i - 1}} \exp\left(\lambda_{ai}(x_i)\right)\sum_{x_a \backslash x_i} f_a(x_a) \exp\Bigg(\sum_{\substack{{j \in N(a)} j \neq i}}\lambda_{aj}(x_j)\Bigg) \\
&= \exp\left(\sum_{a \in N(i)} \lambda_{ai}(x_i)\right) ^{\frac{1}{d_i - 1}} \exp\left(\lambda_{ai}(x_i) + \lambda_{ia}(x_i)\right) \\
&= \exp\left(\sum_{a \in N(i)} \lambda_{ai}(x_i)\right) ^{\frac{1}{d_i - 1}}\, ,
\end{align*}
```
where we defined an auxilary function
```math
\begin{align*}
\exp(\lambda_{ia}(x_i)) \triangleq \sum_{x_a \backslash x_i} f_a(x_a) \exp\Bigg(\sum_{\substack{{j \in N(a)} j \neq i}}\lambda_{aj}(x_j)\Bigg) \,.
\end{align*}
```
This definition is valid since it can be inverted by the relation
```math
\begin{align*}
\lambda_{ia}(x_i) = \frac{2-d_i}{d_i - 1}\lambda_{ai}(x_i) + \frac{1}{d_i -1}\sum_{\substack{c \in N(i)\\c \neq a}}\lambda_{ci}(x_i)
\end{align*}
```
In general it is not possible to solve for the Lagrange multipliers analytically and we resort to iteratively obtaining the solutions. This leads to the **Belief Propagation algorithm** where the exponentiated Lagrange multipliers (messages) are updated iteratively via
```math
\begin{align*}
\mu_{ia}^{(k+1)}(x_i) &= \sum_{x_a \backslash x_i} f_a(x_a) \prod_{\substack{{j \in N(a)} j \neq i}}\mu^{(k)}_{aj}(x_j) \mu_{ai}^{(k)}(x_i) \\
&= \prod_{\substack{c \in N(i) c \neq a}}\mu^{(k)}_{ic}(x_i)\,,
\end{align*}
```
where ``k`` denotes iteration number and the messages are defined as
```math
\begin{align*}
\mu_{ia}(x_i) &\triangleq \exp(\lambda_{ia}(x_i))\\
\mu_{ai}(x_i) &\triangleq \exp(\lambda_{ai}(x_i))\,.
\end{align*}
```
For a more complete overview of message passing as Bethe Free Energy minimization, see [Senoz et al. (2021)](https://research.tue.nl/nl/publications/variational-message-passing-and-local-constraint-manipulation-in-).
"""$26c67f04-d294-11ef-03a4-838ae255689dmetadatadisabled©show_logsîskip_as_script§cell_id$26c67f04-d294-11ef-03a4-838ae255689dcode_foldedäcode "md"""
## The Variational Free Energy Functional
We'll start from scratch. Consider a model ``p(x,z) = p(x|z) p(z)``, where ``x`` and ``z`` are observed and latent variables, respectively. ``z`` may include parameters but also observation-dependent latent variables.
The goal of Bayesian inference is to transform the (known) *likelihood-times-prior* factorization of the full model to a *posterior-times-evidence* decomposition:
```math
\underbrace{p(x|z) p(z)}_{\text{what we know}} \rightarrow \underbrace{p(z|x) p(x)}_{\text{what we want}}
```
Remember from the [Bayesian machine learning lesson](https://bmlip.github.io/course/lectures/Bayesian%20Machine%20Learning.html#Bayesian-model-evidence) that negative log-evidence can be decomposed as "complexity" minus "accuracy" terms (the CA decomposition):
```math
-\log p(x) = \underbrace{ \int p(z|x) \log \frac{p(z|x)}{p(z)} \mathrm{d}z }_{\text{complexity}} - \underbrace{\int p(z|x) \log p(x|z) \mathrm{d}z}_{\text{accuracy}}
```
The CA decomposition cannot be evaluated because it depends on the posterior ``p(z|x)``, which cannot be evaluated since it is the objective of the inference process.
Let's now introduce a distribution ``q(z)`` (called the "variational" or "approximate" posterior distribution) that we will use to *approximate* the posterior ``p(z|x)``. Since we propose the distribution ``q(z)`` ourselves, we will assume that ``q(z)`` can be evaluated.
If we substitute ``q(z)`` for ``p(z|x)`` in the CA decomposition, then we obtain
```math
F[q] \triangleq \underbrace{ \int q(z) \log \frac{q(z)}{p(z)} \mathrm{d}z }_{\text{complexity}} - \underbrace{\int q(z) \log p(x|z) \mathrm{d}z}_{\text{accuracy}}
```
This expression is called the **Variational Free Energy** (VFE), represented by the symbol ``F``. We treat ``F`` as a function of the posterior ``q(z)``. Technically, a function of a function is called a functional, and we write square brackets (e.g., ``F[q]``) to differentiate functionals from functions (e.g., ``q(z)``).
Note that all factors in the CA decomposition of VFE (i.e., ``q(z)``, ``p(z)``, and ``p(x|z)``) can be evaluated as a function of ``z`` (and ``x`` is observed), and therefore the VFE can be evaluated. This is important: log-evidence ``\log p(x)`` cannot be evaluated, but ``F[q]`` *can* be evaluated!
"""$4ee377c2-a126-4c40-8053-517d40c5ef9dmetadatadisabled©show_logsîskip_as_script§cell_id$4ee377c2-a126-4c40-8053-517d40c5ef9dcode_foldedäcodevfem_result = let
#store the inference results in these vectors
ν = fill(3.0, K, max_iterations)
β = fill(1.0, K, max_iterations)
α = fill(0.01, K, max_iterations)
R = Array{Float64}(undef,K,N,max_iterations)
M = Array{Float64}(undef,2,K,max_iterations)
Λ = Array{Float64}(undef,2,2,K,max_iterations)
clusters_vb = Array{Distribution}(undef,K,max_iterations) #clusters to be plotted
#initialize prior distribution parameters
M[:,:,1] = rand(MersenneTwister(43), 2, K) .* [4, 50] .+ [1, 50]
for k in 1:K
Λ[:,:,k,1] = [1.0 0;0 0.01]
R[k,:,1] = 1/(K)*ones(N)
clusters_vb[k,1] = MvNormal(M[:,k,1],PDMats.PDMat(convert(Matrix,Hermitian(inv(ν[1,1].*Λ[:,:,k,1])))))
end
#variational inference
for i in 1:max_iterations-1
#variational expectation
R[:,:,i+1] = updateR(Λ[:,:,:,i],M[:,:,i],α[:,i],ν[:,i],β[:,i])
#variational minimisation
M[:,:,i+1],β[:,i+1],Λ[:,:,:,i+1],ν[:,i+1],α[:,i+1] = updateMeanPrecisionPi(M[:,:,i],β[:,i],Λ[:,:,:,i],ν[:,i],α[:,i],R[:,:,i+1])
for k in 1:K
clusters_vb[k,i+1] = MvNormal(M[:,k,i+1],PDMats.PDMat(convert(Matrix,Hermitian(inv(ν[k,i+1].*Λ[:,:,k,i+1])))))
end
end
(; clusters_vb, R)
end$26c5cfb4-d294-11ef-05bb-59d5e27cf37cmetadatadisabled©show_logsîskip_as_script§cell_id$26c5cfb4-d294-11ef-05bb-59d5e27cf37ccode_foldedäcodeZmd"""
## The Marginal Distribution for the GMM
In the literature, the GMM is often introduced by the marginal distribution for an *observed* data point ``x_n``, given by
```math
\begin{align*}{}
p(x_n) &= \sum_{z_n} p(x_n,z_n) \\
&= \sum_{k=1}^K \pi_k \cdot \mathcal{N}\left( x_n | \mu_k, \Sigma_k \right) \tag{B-9.12}
\end{align*}
```
"""$a7fb83cb-1f40-4c8a-9fda-2165f91e413emetadatadisabled©show_logsîskip_as_script§cell_id$a7fb83cb-1f40-4c8a-9fda-2165f91e413ecode_foldedäcodelet
(i, step), res = em_result[iteration_em]
# Get the current iteration's parameters
current_clusters = res.clusters
current_responsibilities = res.γ
# Compute mixture weights (π) from responsibilities
π = vec(mean(current_responsibilities, dims=2))
if any(isnan, π)
π = ones(size(π))
end
π = π ./ sum(π) # Normalize to ensure they sum to 1
# Generate synthetic data of the same size as original
N_synthetic = size(X, 2) # Same number of points as original data
X_synthetic = Matrix{Float64}(undef, 2, N_synthetic)
# Sample from the mixture model
Random.seed!(42) # For reproducibility
components = rand(Categorical(π), N_synthetic)
for n in 1:N_synthetic
# Sample which component to use
# Sample from that component
X_synthetic[:, n] = rand(current_clusters[components[n]])
end
# Create the plot
scatter(
X_synthetic[1,:], X_synthetic[2,:];
title="Synthetic data from this model",
size=(600,300),
legend=nothing,
plot_lims...,
data_plot_kwargs...,
color=show_em_synthetic_labels ? cluster_colors[components] : :purple,
)
end$64819124-865e-48b8-a916-2ce08dba0accmetadatadisabled©show_logsîskip_as_script§cell_id$64819124-865e-48b8-a916-2ce08dba0acccode_foldedäcodeh@htl """
Show class colors: $(@bind show_em_synthetic_labels CheckBox(default=false))
"""$16d90f11-5933-4145-b219-19774eba25d6metadatadisabled©show_logsîskip_as_script§cell_id$16d90f11-5933-4145-b219-19774eba25d6code_foldedäcodet@htl """
Show synthetic label colors: $(@bind show_vfem_synthetic_labels CheckBox(default=false))
"""$06170e31-e865-4178-8af0-41d82df95d71metadatadisabled©show_logsîskip_as_script§cell_id$06170e31-e865-4178-8af0-41d82df95d71code_foldedäcode`keyconcept("","Constrained VFE minimization is consistent with the **Maximum Entropy Principle**, which prescribes how to rationally update beliefs when new information becomes available. In this framework, the updated posterior is the distribution that minimizes VFE (or equivalently, KL divergence to the prior) subject to the imposed constraints. ")$f1f7407d-86a1-4f24-b78a-61a411d1f371metadatadisabled©show_logsîskip_as_script§cell_id$f1f7407d-86a1-4f24-b78a-61a411d1f371code_foldedäcode!md"""
# Variational Inference
"""$5af3ff1b-1655-4dd9-a089-91544fc85a0emetadatadisabled©show_logsîskip_as_script§cell_id$5af3ff1b-1655-4dd9-a089-91544fc85a0ecode_foldedäcodeh@bindname iteration_em Slider(eachindex(em_result); default=2, show_value=i->join(em_result[i][1], "-"))$8555aec9-4e80-49e7-8514-ef4a2236801bmetadatadisabled©show_logsîskip_as_script§cell_id$8555aec9-4e80-49e7-8514-ef4a2236801bcode_folded¤codeN = size(X, 2)$26c65092-d294-11ef-39cc-1953a725f285metadatadisabled©show_logsîskip_as_script§cell_id$26c65092-d294-11ef-39cc-1953a725f285code_foldedäcodeZmd"""
Can we still compute an approximate posterior? In this lesson, we introduce an approximate Bayesian inference method known as **Variational Bayes** (VB) (also known as **Variational Inference**) that can be used for Bayesian inference in models with latent variables. Later in this lesson, we will use VB to do inference in the GMM.
"""$3948225e-90d8-4b78-ba9b-5f98e228285ametadatadisabled©show_logsîskip_as_script§cell_id$3948225e-90d8-4b78-ba9b-5f98e228285acode_foldedäcodeBscatter(X[1,:], X[2,:]; label="observations", data_plot_kwargs...)$95d47a10-f3f8-479b-afe0-21241104b758metadatadisabled©show_logsîskip_as_script§cell_id$95d47a10-f3f8-479b-afe0-21241104b758code_foldedäcode7code_example("VFEM for GMM on Old Faithfull data set")
$26c62ebe-d294-11ef-0cfb-ef186203e890metadatadisabled©show_logsîskip_as_script§cell_id$26c62ebe-d294-11ef-0cfb-ef186203e890code_foldedäcodemd"""
Since the class labels ``y_{nk} \in \{0,1\}`` were assumed to be given by the data set, maximization of this expression decomposed into a set of simple update rules for the Gaussian and multinomial distributions.
"""$c4b23b39-e6e4-44ec-b204-0c7d7d5a4026metadatadisabled©show_logsîskip_as_script§cell_id$c4b23b39-e6e4-44ec-b204-0c7d7d5a4026code_foldedäcode]md"""
Since the second term (log-evidence) does not involve ``q(z)``, VFEM over ``q`` will bring ``q(z)`` closer to the Bayesian posterior ``p(z|x)``.
Since ``D_{\text{KL}}[q(z),p(z|x)]\geq 0`` for any ``q(z)``, and ``D_{\text{KL}}[q(z),p(z|x)]= 0`` only if ``q(z) = p(z|x)``, the VFE is always an upper bound on (minus) log-evidence, i.e.,
```math
F[q] \geq -\log p(x) \,.
```
As a result, global minimization of VFE leads to
```math
q^*(z) = \arg\min_q F[q]
```
where
```math
\begin{align}
q^*(z) &= p(z|x) \tag{= posterior}\\
F[q^*] &= -\log p(x) \tag{= -log evidence}
\end{align}
```
"""$26c5b896-d294-11ef-1d8e-0feb99d2d45bmetadatadisabled©show_logsîskip_as_script§cell_id$26c5b896-d294-11ef-1d8e-0feb99d2d45bcode_foldedäcode"md"""
We associate a one-hot coded hidden class label ``z_n`` with each observation ``x_n``:
```math
\begin{equation*}
z_{nk} = \begin{cases} 1 & \text{if } x_n \in \mathcal{C}_k \text{ (the $k$-th class)}\\
0 & \text{otherwise} \end{cases}
\end{equation*}
```
"""$26c623f6-d294-11ef-13c0-19edd43592c0metadatadisabled©show_logsîskip_as_script§cell_id$26c623f6-d294-11ef-13c0-19edd43592c0code_foldedäcodemd"""
## Inference for GMM is Difficult
Indeed, the fact that the observation-dependent class labels are *unobserved* for the GMM, leads to a problem for processing new data by Bayes rule in a GMM.
Consider a given data set ``D = \{x_1,x_2,\ldots,x_N\}``. We recall here the log-likelihood for the Gaussian-Categorial Model, see the [generative classification lesson](https://bmlip.github.io/course/lectures/Generative%20Classification.html):
```math
\log\, p(D|\theta) = \sum_{n,k} y_{nk} \underbrace{ \log\mathcal{N}(x_n|\mu_k,\Sigma) }_{ \text{Gaussian} } + \underbrace{ \sum_{n,k} y_{nk} \log \pi_k }_{ \text{multinomial} } \,.
```
"""$55570464-89c8-4d9b-b667-dfa64ac62294metadatadisabled©show_logsîskip_as_script§cell_id$55570464-89c8-4d9b-b667-dfa64ac62294code_foldedäcodemd"""
# Code
"""$26c59b52-d294-11ef-1eba-d3f235f85eeemetadatadisabled©show_logsîskip_as_script§cell_id$26c59b52-d294-11ef-1eba-d3f235f85eeecode_foldedäcodemd"""
## Unobserved Classes
Consider again a set of observed data ``D=\{x_1,\dotsc,x_N\}``.
This time, we suspect that there are *unobserved* class labels that would help explain (or predict) the data, e.g.,
* the observed data are the color of living things; the unobserved classes are animals and plants.
* observed are wheel sizes; unobserved categories are trucks and personal cars.
* observed is an audio signal; unobserved classes include speech, music, traffic noise, etc.
"""$f17c9d8a-9291-4110-bcf4-c582d23f986bmetadatadisabled©show_logsîskip_as_script§cell_id$f17c9d8a-9291-4110-bcf4-c582d23f986bcode_foldedäcodeٌmd"""
## Variational Message Passing
Let us now consider the local free energy that is associated with edge corresponding to ``x_j``.
"""$dd1242db-fb20-4732-ac55-a3e021bbd2b7metadatadisabled©show_logsîskip_as_script§cell_id$dd1242db-fb20-4732-ac55-a3e021bbd2b7code_folded¤code^const data_plot_kwargs = (
markersize=4,
markerstrokewidth=0,
color=:red4,
opacity=.8,
)$26c73cf0-d294-11ef-297b-354eb9c71f57metadatadisabled©show_logsîskip_as_script§cell_id$26c73cf0-d294-11ef-297b-354eb9c71f57code_foldedäcodemd"""
## Derivation of VFEM Update Equations
Let's get back to the illustrative challenge at the beginning of this lesson: we want to do [density modeling for the Old Faithful data set](#challenge_hello).
"""$2d4adbf6-6de8-4e3a-ad6f-fa8bbfa5999emetadatadisabled©show_logsîskip_as_script§cell_id$2d4adbf6-6de8-4e3a-ad6f-fa8bbfa5999ecode_foldedäcodemd"""
#### Mean Updating (*)
Explain the following update rule for the [mean of the Gaussian cluster-conditional data distribution](##Update-Equations):
```math
m_k = \frac{1}{\beta_k} \left( \beta_0 m_0 + N_k \bar{x}_k \right) \tag{B-10.61}
```
"""$26c56fd8-d294-11ef-236d-81deef63f37cmetadatadisabled©show_logsîskip_as_script§cell_id$26c56fd8-d294-11ef-236d-81deef63f37ccode_foldedäcode5title("Latent Variable Models and Variational Bayes")$666680b2-315a-4d95-8f7f-3ae50018e112metadatadisabled©show_logsîskip_as_script§cell_id$666680b2-315a-4d95-8f7f-3ae50018e112code_folded¤codeK = 6;$0349720e-5de4-4b39-babd-c0881588f1demetadatadisabled©show_logsîskip_as_script§cell_id$0349720e-5de4-4b39-babd-c0881588f1decode_foldedäcode)X = Array(Matrix{Float64}(old_faithful)')$cd7cf2d0-aff7-49dc-bb12-dc73c6433768metadatadisabled©show_logsîskip_as_script§cell_id$cd7cf2d0-aff7-49dc-bb12-dc73c6433768code_foldedäcodelet
i = iteration_vfem + 1
# Get the current iteration's parameters
current_clusters = vfem_result.clusters_vb[:, i]
current_responsibilities = vfem_result.R[:, :, i]
# Compute mixture weights (π) from responsibilities
π = vec(mean(current_responsibilities, dims=2))
π = π ./ sum(π) # Normalize to ensure they sum to 1
# Generate synthetic data of the same size as original
N_synthetic = size(X, 2) # Same number of points as original data
X_synthetic = Matrix{Float64}(undef, 2, N_synthetic)
# Sample from the mixture model
Random.seed!(42) # For reproducibility
components = rand(Categorical(π), N_synthetic)
for n in 1:N_synthetic
# Sample which component to use
# Sample from that component
X_synthetic[:, n] = rand(current_clusters[components[n]])
end
# Create the plot
scatter(
X_synthetic[1,:], X_synthetic[2,:];
title="Synthetic data from this model",
size=(600,300),
legend=nothing,
plot_lims...,
data_plot_kwargs...,
color=show_vfem_synthetic_labels ? cluster_colors[components] : :purple,
)
end$b4d965d1-91d8-43f8-84a5-b37ad5c6cafametadatadisabled©show_logsîskip_as_script§cell_id$b4d965d1-91d8-43f8-84a5-b37ad5c6cafacode_foldedäcodemd"""
Note that the inference bound is a [Kullback-Leibler (KL) divergence](https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence) between an (approximate) posterior ``q(z)`` and the (perfect) Bayesian posterior ``p(z|x)``. To learn more:
"""$32f0bbb4-dfc7-431e-9a3d-80162439edacmetadatadisabled©show_logsîskip_as_script§cell_id$32f0bbb4-dfc7-431e-9a3d-80162439edaccode_foldedäcode md"""
Apparently (see previous slide), there are three contributions to the free energy for ``x_j``:
* one entropy term for the edge ``x_j``
* two energy terms: one for each node that attaches to ``x_j`` (in the figure: nodes ``p_a`` and ``p_b``)
The local free energy for ``x_j`` can be written as (exercise)
```math
F[q_j] \propto \sum_{x_j} q(x_j) \log \frac{q_j(x_j)}{\nu_a(x_j)\cdot \nu_b(x_j)}
```
where
```math
\begin{align*}
\nu_a(x_j) &\propto \exp\left( \mathbb{E}_{q_{k}}\left[ \log p_a(x_a)\right]\right) \\
\nu_b(x_j) &\propto \exp\left( \mathbb{E}_{q_{l}}\left[ \log p_b(x_b)\right]\right)
\end{align*}
```
and ``\mathbb{E}_{q_{k}}\left[\cdot\right]`` is an expectation w.r.t. all ``q(x_k)`` with ``k \in N(a)\setminus {j}``.
``\nu_a(x_j)`` and ``\nu_b(x_j)`` can be locally computed in nodes ``a`` and ``b`` respectively and can be interpreted as colliding messages over edge ``x_j``.
Local free energy minimization is achieved by setting
```math
q_j(x_j) \propto \nu_a(x_j) \cdot \nu_b(x_j)
```
Note that message ``\nu_a(x_j)`` depends on posterior beliefs over incoming edges (``k``) for node ``a``, and in turn, the message from node ``a`` towards edge ``x_k`` depends on the belief ``q_j(x_j)``. I.o.w., direct mutual dependencies exist between posterior beliefs over edges that attach to the same node.
These considerations lead to the [Variational Message Passing](https://en.wikipedia.org/wiki/Variational_message_passing) procedure, which is an iterative free energy minimization procedure that can be executed completely through locally computable messages.
Procedure VMP, see [Dauwels (2007), section 3](https://github.com/bmlip/course/blob/main/assets/files/Dauwels-2007-on-variational-message-passing-on-factor-graphs.pdf)
> 1. Initialize all messages ``q`` and ``ν``, e.g., ``q(\cdot) \propto 1`` and ``\nu(\cdot) \propto 1``. FE Minimization by the Expectation-Maximization (EM) Algorithm "))
The EM algorithm is a special case of VFE minimization that focuses on Maximum-Likelihood estimation for models with latent variables.
Consider a model
```math
p(x,z,\theta)
```
with observations ``x = \{x_n\}``, latent variables ``z=\{z_n\}`` and parameters ``\theta``.
We can write the following VFE functional for this model:
```math
\begin{align*}
F[q] = \sum_z \sum_\theta q(z,\theta) \log \frac{q(z,\theta)}{p(x,z,\theta)}
\end{align*}
```
The EM algorithm makes the following simplifying assumptions:
1. The prior for the parameters is uninformative (uniform). This implies that
```math
p(x,z,\theta) = p(x,z|\theta) p(\theta) \propto p(x,z|\theta)
```
2. A factorization constraint
```math
q(z,\theta) = q(z) q(\theta)
```
3. The posterior for the parameters is a delta function:
```math
q(\theta) = \delta(\theta - \hat{\theta})
```
Basically, these three assumptions turn VFE minimization into maximum likelihood estimation for the parameters ``\theta`` and the VFE simplifies to
```math
\begin{align*}
F[q,\theta] = \sum_z q(z) \log \frac{q(z)}{p(x,z|\theta)}
\end{align*}
```
The EM algorithm minimizes this FE by iterating (iteration counter: ``i``) over
```math
\begin{align} \mathcal{L}^{(i)}(\theta) &= \sum_z \overbrace{p(z|x,\theta^{(i-1)})}^{q^{(i)}(z)} \log p(x,z|\theta) \tag{the E-step} \\
\theta^{(i)} &= \arg\max_\theta \mathcal{L}^{(i)}(\theta) \tag{the M-step} \end{align}
```
These choices are optimal for the given FE functional. In order to see this, consider the two decompositions
```math
\begin{align*}
F[q,\theta] &= \underbrace{-\sum_z q(z) \log p(x,z|\theta)}_{\text{energy}} - \underbrace{\sum_z q(z) \log \frac{1}{q(z)}}_{\text{entropy}} \qquad &&\text{(EE)}\\
&= \underbrace{\sum_z q(z) \log \frac{q(z)}{p(z|x,\theta)}}_{\text{divergence}} - \underbrace{\log p(x|\theta)}_{\text{log-likelihood}} \qquad &&\text{(DE)}
\end{align*}
```
The DE decomposition shows that the FE is minimized for the choice ``q(z) := p(z|x,\theta)``. Also, for this choice, the FE equals the (negative) log-evidence (, which is this case simplifies to the log-likelihood).
The EE decomposition shows that the FE is minimized wrt ``\theta`` by minimizing the energy term. The energy term is computed in the E-step and optimized in the M-step.
* Note that in the EM literature, the energy term is often called the *expected complete-data log-likelihood*.)
In order to execute the EM algorithm, it is assumed that we can analytically execute the E- and M-steps. For a large set of models (including models whose distributions belong to the exponential family of distributions), this is indeed the case and hence the large popularity of the EM algorithm.
The EM algorihm imposes rather severe assumptions on the FE (basically approximating Bayesian inference by maximum likelihood estimation). Over the past few years, the rise of Probabilistic Programming languages has dramatically increased the range of models for which the parameters can by estimated autmatically by (approximate) Bayesian inference, so the popularity of EM is slowly waning. (More on this in the Probabilistic Programming lessons).
Bishop (2006) works out EM for the GMM in section 9.2.
"""$d208f60e-56ea-4767-bd18-f29a853b7536metadatadisabled©show_logsîskip_as_script§cell_id$d208f60e-56ea-4767-bd18-f29a853b7536code_folded¤code3const cluster_colors = color_list(:seaborn_bright6)$26c7696e-d294-11ef-25f2-dbc0946c0858metadatadisabled©show_logsîskip_as_script§cell_id$26c7696e-d294-11ef-25f2-dbc0946c0858code_foldedäcode١md"""
Below we exemplify training of a Gaussian Mixture Model on the Old Faithful data set by VFE minimization, with the constraints as specified above.
"""$58bd0d43-743c-4745-b353-4a89b35e85bametadatadisabled©show_logsîskip_as_script§cell_id$58bd0d43-743c-4745-b353-4a89b35e85bacode_folded¤code&using Distributions, Plots, StatsPlots$b91bc3b6-b815-4942-b297-c0e2b4b99654metadatadisabled©show_logsîskip_as_script§cell_id$b91bc3b6-b815-4942-b297-c0e2b4b99654code_foldedäcodeQhide_solution(
md"""
Let's start with a prior estimate ``\theta^{(i-1)}`` and we want to minimize the free energy functional wrt ``q``. This leads to
```math
\begin{align*}
q^{(i)}(z) &= \arg\min_q F[q](\theta^{(i-1)}) \\
&= \arg\min_q \sum_z q(z) \log \frac{q(z)}{p(D,z|\theta^{(i-1)})} \\
&= \arg\min_q \sum_z q(z) \log \frac{q(z)}{p(z|D,\theta^{(i-1)}) \cdot p(D|\theta^{(i-1)})} \\
&= p(z|D,\theta^{(i-1)})
\end{align*}
```
Next, we use ``q^{(i)}(z)=p(z|D,\theta^{(i-1)})`` and minimize the free energy w.r.t. ``\theta``, leading to
```math
\begin{align*}
\theta^{(i)} &= \arg\min_\theta F[q^{(i)}](\theta) \\
&= \arg\min_\theta \sum_z p(z|D,\theta^{(i-1)}) \log \frac{p(z|D,\theta^{(i-1)})}{p(D,z|\theta)} \\
&= \arg\max_\theta \sum_z \underbrace{p(z|D,\theta^{(i-1)})}_{q^{(i)}(z)} \log p(D,z|\theta)
\end{align*}
```
""")$87d94630-c90b-4379-91bd-88641ee7b508metadatadisabled©show_logsîskip_as_script§cell_id$87d94630-c90b-4379-91bd-88641ee7b508code_foldedäcodeٗmd"""
### Implementation
_Reading this lecture online? Click **"View code"** in the top right to read the implementation of this visualisation._
"""$56bea391-b812-4fc4-8f27-fcb4cb984cf4metadatadisabled©show_logsîskip_as_script§cell_id$56bea391-b812-4fc4-8f27-fcb4cb984cf4code_foldedäcodeexercises(header_level=1)$de0c41a3-6319-4ae3-8a1a-ae6935910fa3metadatadisabled©show_logsîskip_as_script§cell_id$de0c41a3-6319-4ae3-8a1a-ae6935910fa3code_foldedäcodea@bindname iteration_vfem Slider([0:10..., 10:10:max_iterations-1...]; default=2, show_value=true)$a663d6d1-2d73-40c6-8a4c-9e41df74bc84metadatadisabled©show_logsîskip_as_script§cell_id$a663d6d1-2d73-40c6-8a4c-9e41df74bc84code_foldedäcodemd"""
##### Problem
You are asked to build a density model for a data set ([Old Faithful](https://en.wikipedia.org/wiki/Old_Faithful), Bishop pg. 681) that clearly is not distributed as a single Gaussian:
"""$26c5a1f6-d294-11ef-3565-39d027843fbbmetadatadisabled©show_logsîskip_as_script§cell_id$26c5a1f6-d294-11ef-3565-39d027843fbbcode_foldedäcodemd"""
Classification problems with unobserved classes are called **Clustering** problems. In clustering problems, the learning algorithm needs to *discover the underlying classes from the observed data*.
"""$2f490e1f-e495-4f55-a3f8-60d6fd716d4emetadatadisabled©show_logsîskip_as_script§cell_id$2f490e1f-e495-4f55-a3f8-60d6fd716d4ecode_foldedäcodemd"""
#### The Expectation-Maximization (EM) algorithm (**)
Consider a model ``p(x,z|\theta)``, where ``D=\{x_1,x_2,\ldots,x_N\}`` is observed, ``z`` are unobserved variables, and ``\theta`` are parameters. The **Expectation-Maximization** (EM) algorithm estimates the parameters by iterating over the following two equations (``i`` is the iteration index):
```math
\begin{align*}
q^{(i)}(z) &= p(z|D,\theta^{(i-1)}) \\
\theta^{(i)} &= \arg\max_\theta \sum_z q^{(i)}(z) \cdot \log p(D,z|\theta)
\end{align*}
```
Prove that this algorithm minimizes the Free Energy functional
```math
\begin{align*}
F[q](\theta) = \sum_z q(z) \log \frac{q(z)}{p(D,z|\theta)}
\end{align*}
```
"""$26c796c8-d294-11ef-25be-17dcd4a9d315metadatadisabled©show_logsîskip_as_script§cell_id$26c796c8-d294-11ef-25be-17dcd4a9d315code_foldedäcodeTmd"""
The generated figure resembles Figure 10.6 in Bishop. The plots show VFEM results for a GMM of ``K = 6`` Gaussians applied to the Old Faithful data set. The ellipses denote the two standard-deviation density contours for each of the components, and the color coding of the data points reflects the "soft" class label assignments.
"""$98a0ed70-a627-48d6-a1f8-3dec7aba2bb2metadatadisabled©show_logsîskip_as_script§cell_id$98a0ed70-a627-48d6-a1f8-3dec7aba2bb2code_foldedäcode@function updateMeanPrecisionPi(m_0,β_0,W_0,ν_0,α_0,r) #variational maximisation function
m = Array{Float64}(undef,2,K) #mean of the clusters
β = Array{Float64}(undef,K) #precision scaling for Gausian distribution
W = Array{Float64}(undef,2,2,K) #precision prior for Wishart distributions
ν = Array{Float64}(undef,K) #degrees of freedom parameter for Wishart distribution
α = Array{Float64}(undef,K) #Dirichlet distribution parameter
for k=1:K
sst = sufficientStatistics(X,r,k)
α[k] = α_0[k] + sst[1]; β[k] = β_0[k] + sst[1]; ν[k] = ν_0[k] .+ sst[1]
m[:,k] = (1/β[k])*(β_0[k].*m_0[:,k] + sst[1].*sst[2])
W[:,:,k] = inv(inv(W_0[:,:,k])+sst[3]*sst[1] + ((β_0[k]*sst[1])/(β_0[k]+sst[1])).*(sst[2]-m_0[:,k])*(sst[2]-m_0[:,k])')
end
return m,β,W,ν,α
end$26c867d8-d294-11ef-2372-d75ed0bcc02dmetadatadisabled©show_logsîskip_as_script§cell_id$26c867d8-d294-11ef-2372-d75ed0bcc02dcode_foldedäcode١md"""
We'll perform clustering on the data set from the [challenge](#challenge_hello) by fitting a GMM consisting of two Gaussians using the EM algorithm.
"""$26c8160c-d294-11ef-2a74-6f7009a7c51emetadatadisabled©show_logsîskip_as_script§cell_id$26c8160c-d294-11ef-2a74-6f7009a7c51ecode_foldedäcodemd"""
# Optional Slides
"""$26c64174-d294-11ef-2bbc-ab1a84532311metadatadisabled©show_logsîskip_as_script§cell_id$26c64174-d294-11ef-2bbc-ab1a84532311code_foldedäcodemd"""
There is no known conjugate prior for the latent variables in the GMM likelihood. Therefore, Bayes rule does not yield a closed-form expression for the posterior over the latent variables:
```math
\underbrace{p(\{z_{nk}\},\{\mu_k,\Sigma_k,\pi_k\} | D)}_{\text{posterior (no analytical solution)}} \propto \underbrace{p(D\,|\,\{z_{nk}\},\{\mu_k,\Sigma_k,\pi_k\})}_{\text{likelihood}} \cdot \underbrace{p( \{z_{nk}\},\{\mu_k,\Sigma_k,\pi_k\} )}_{\text{prior (no known conjugate)}}
```
"""$208ba1bb-a4bf-4b8c-93d2-0d6c6c8d16d4metadatadisabled©show_logsîskip_as_script§cell_id$208ba1bb-a4bf-4b8c-93d2-0d6c6c8d16d4code_foldedäcodehide_solution(
md"""
We see here an example of "precision-weighted means add" when two sources of information are fused, just like precision-weighted means add when two Gaussians are multiplied, eg a prior and likelihood. In this case, the prior is ``m_0`` and the likelihood estimate is ``\bar{x}``. ``\beta_0`` can be interpreted as the number of pseudo-observations in the prior.
""")$26c6f63c-d294-11ef-1090-e9238dd6ad3fmetadatadisabled©show_logsîskip_as_script§cell_id$26c6f63c-d294-11ef-1090-e9238dd6ad3fcode_foldedäcodemd"""
## Constrained VFE Minimization
It is common to add simplifying constraints to an optimization problem to make a difficult optimization task tractible. This is also standard practice when approximating Bayesian inference by FE minimization.
There are three important cases of adding constraints to the VFE functional that often alleviate the VFE minimization task:
- form constraints on ``q(z)``
- factorization constraints on ``q(z)``
- other (ad hoc) constraints
We will shortly discuss these simplifications below.
"""$aea77d69-9ecd-4be0-b6fd-c944d27d68dfmetadatadisabled©show_logsîskip_as_script§cell_id$aea77d69-9ecd-4be0-b6fd-c944d27d68dfcode_foldedäcodemd"""
##### 1. Form constraints
For almost every practical setting, we constrain the posterior ``q(z)`` to belong to a **specific parameterized family** of probability distributions, e.g.,
```math
q(z) = \mathcal{N}\left( z | \mu, \Sigma \right)\,.
```
In this case, the *functional* minimization problem for ``F[q]`` simplifies to the minimization of an ordinary *function*
```math
F(\mu,\Sigma) = \int \mathcal{N}\left( z | \mu, \Sigma \right) \log \frac{\mathcal{N}\left( z | \mu, \Sigma \right)}{p(x,z)}\mathrm{d}z
```
with respect to the parameters ``\mu`` and ``\Sigma``.
We can often use standard gradient-based optimization methods to minimize the function ``F(\mu,\Sigma)\,.``
"""$26c82f16-d294-11ef-0fe1-07326b56282fmetadatadisabled©show_logsîskip_as_script§cell_id$26c82f16-d294-11ef-0fe1-07326b56282fcode_foldedäcode nmd"""
## VFE Minimization with Mean-field Factorization Constraints: $(HTML("the CAVI Approach "))
Let's work out VFE minimization with additional mean-field constraints (=full factorization) constraints:
```math
q(z) = \prod_{j=1}^m q_j(z_j)\,.
```
In other words, the posteriors for ``z_j`` are all considered independent. This is a strong constraint but often leads to good solutions.
Given the mean-field constraints, it is possible to derive the following expression for the $(HTML("optimal solutions ")) ``q_j^*(z_j)``, for ``j=1,\ldots,m``:
```math
\begin{align}
\log q_j^*(z_j) &\propto \mathrm{E}_{q_{-j}^*}\left[ \log p(x,z) \right] \\
&= \underbrace{\sum_{z_{-j}} q_{-j}^*(z_{-j}) \underbrace{\log p(x,z)}_{\text{"field"}}}_{\text{"mean field"}}
\end{align}
```
where we defined ``q_{-j}^*(z_{-j}) \triangleq q_1^*(z_1)q_2^*(z_2)\cdots q_{j-1}^*(z_{j-1})q_{j+1}^*(z_{j+1})\cdots q_m^*(z_m)``.
**Proof** (from [Blei, 2017](https://doi.org/10.1080/01621459.2017.1285773)): We first rewrite the FE as a function of ``q_j(z_j)`` only:
```math
F[q_j] = \mathbb{E}_{q_{j}}\left[ \mathbb{E}_{q_{-j}}\left[ \log p(x,z_j,z_{-j})\right]\right] - \mathbb{E}_{q_j}\left[ \log q_j(z_j)\right] + \mathtt{const.}\,,
```
where the constant holds all terms that do not depend on ``z_j``. This expression can be written as
```math
F[q_j] = \sum_{z_j} q_j(z_j) \log \frac{q_j(z_j)}{\exp\left( \mathbb{E}_{q_{-j}}\left[ \log p(x,z_j,z_{-j})\right]\right)}
```
which is a KL-divergence that is minimized by Eq. B-10.9. (end proof)
This is not yet a full solution to the FE minimization task since the solution ``q_j^*(z_j)`` depends on expectations that involve other solutions ``q_{i\neq j}^*(z_{i \neq j})``, and each of these other solutions ``q_{i\neq j}^*(z_{i \neq j})`` depends on an expection that involves ``q_j^*(z_j)``.
In practice, we solve this chicken-and-egg problem by an iterative approach: we first initialize all ``q_j(z_j)`` (for ``j=1,\ldots,m``) to an appropriate initial distribution and then cycle through the factors in turn by solving eq.B-10.9 and update ``q_{-j}^*(z_{-j})`` with the latest estimates. (See [Blei, 2017](https://doi.org/10.1080/01621459.2017.1285773), Algorithm 1, p864).
This algorithm for approximating Bayesian inference is known **Coordinate Ascent Variational Inference** (CAVI).
"""$87715fcd-c92e-4bfa-8b0a-15e469548f3dmetadatadisabled©show_logsîskip_as_script§cell_id$87715fcd-c92e-4bfa-8b0a-15e469548f3dcode_foldedäcode$keyconcept("",
md"""
For a model ``p(x,z) = p(x|z) p(x)``, the **variational free energy** functional
```math
F[q] = \underbrace{ \int q(z) \log \frac{q(z)}{p(z)} \mathrm{d}z }_{\text{complexity}} - \underbrace{\int q(z) \log p(x|z) \mathrm{d}z}_{\text{accuracy}}
```
is considered a functional of a variational posterior ``q(z)`` over the latent variables ``z``. Minimizing ``F[q]`` accomplishes **approximate Bayesian inference** by
```math
\begin{align*}
\hat{q}(z) &\approx p(z|x) \\
F[\hat{q}] &\approx -\log p(x)
\end{align*}
```
""")last_hot_reload_time shortpath Latent Variable Models and VB.jlnbpkgbusy_packages,waiting_for_permission_but_probably_disabled§enabledðterminal_outputsLinearAlgebrap
[0m[1mWaiting for other notebooks to finish Pkg operations...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[36m[1m Project[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Project.toml`
[36m[1m Manifest[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Manifest.toml`
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mPrecompiling...[22m
[90m===[39m
Waiting for notebook process to start... Done. Starting precompilation...PDMatsp
[0m[1mWaiting for other notebooks to finish Pkg operations...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[36m[1m Project[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Project.toml`
[36m[1m Manifest[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Manifest.toml`
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mPrecompiling...[22m
[90m===[39m
Waiting for notebook process to start... Done. Starting precompilation...Randomp
[0m[1mWaiting for other notebooks to finish Pkg operations...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[36m[1m Project[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Project.toml`
[36m[1m Manifest[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Manifest.toml`
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mPrecompiling...[22m
[90m===[39m
Waiting for notebook process to start... Done. Starting precompilation...SpecialFunctionsp
[0m[1mWaiting for other notebooks to finish Pkg operations...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[36m[1m Project[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Project.toml`
[36m[1m Manifest[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Manifest.toml`
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mPrecompiling...[22m
[90m===[39m
Waiting for notebook process to start... Done. Starting precompilation...nbpkg_syncp
[0m[1mWaiting for other notebooks to finish Pkg operations...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[36m[1m Project[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Project.toml`
[36m[1m Manifest[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Manifest.toml`
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mPrecompiling...[22m
[90m===[39m
Waiting for notebook process to start... Done. Starting precompilation...BmlipTeachingToolsp
[0m[1mWaiting for other notebooks to finish Pkg operations...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[36m[1m Project[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Project.toml`
[36m[1m Manifest[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Manifest.toml`
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mPrecompiling...[22m
[90m===[39m
Waiting for notebook process to start... Done. Starting precompilation...Distributionsp
[0m[1mWaiting for other notebooks to finish Pkg operations...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[36m[1m Project[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Project.toml`
[36m[1m Manifest[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Manifest.toml`
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mPrecompiling...[22m
[90m===[39m
Waiting for notebook process to start... Done. Starting precompilation...Plotsp
[0m[1mWaiting for other notebooks to finish Pkg operations...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[36m[1m Project[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Project.toml`
[36m[1m Manifest[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Manifest.toml`
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mPrecompiling...[22m
[90m===[39m
Waiting for notebook process to start... Done. Starting precompilation...StatsPlotsp
[0m[1mWaiting for other notebooks to finish Pkg operations...[22m
[90m===[39m
[0m[1mResolving...[22m
[90m===[39m
[36m[1m Project[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Project.toml`
[36m[1m Manifest[22m[39m No packages added to or removed from `~/.julia/scratchspaces/c3e4b0f8-55cb-11ea-2926-15256bba5781/pkg_envs/env_zgtbzjsicr/Manifest.toml`
[0m[1mInstantiating...[22m
[90m===[39m
[0m[1mPrecompiling...[22m
[90m===[39m
Waiting for notebook process to start... Done. Starting precompilation...waiting_for_permission·restart_recommended_msgrestart_required_msginstalled_versionsLinearAlgebrastdlibPDMats0.11.41__internal_julia_version1.13.3RandomstdlibSpecialFunctions2.9.1BmlipTeachingTools1.4.1Distributions0.25.132Plots1.41.7StatsPlots0.15.8!__internal_julia_manifest_version1.13.3install_time_nsXu:instantiatedípluto_versionv1.1.0 
persist_js_state$f8c8013a-3e87-4d01-a3ae-86b39cf1f002running§runtime Âdepends_on_skipped_cellsµpublished_object_keyserrored§cell_id$f8c8013a-3e87-4d01-a3ae-86b39cf1f002depends_on_disabled_cells¦queued¤logsoutputmimetext/htmlrootassigneelast_run_timestampAڹ&has_pluto_hook_features¤body`
