# Variational Bayes Made Easy

Mohammad Emtiyaz Khan<EMAIL_ADDRESS>
RIKEN Center for Advanced Intelligence Project
1-4-1 Nihonbashi, Chuo-ku, Tokyo 103-0027, Japan

###### Abstract

Variational Bayes is a popular method for approximate inference but its
derivation can be cumbersome. To simplify the process, we give a 3-step recipe
to identify the posterior form by explicitly looking for linearity with
respect to expectations of well-known distributions. We can then directly
write the update by simply “reading-off” the terms in front of those
expectations. The recipe makes the derivation easier, faster, shorter, and
more general.

##  1 Introduction

Since its introduction in the early 90s [Hinton and Van Camp, 1993, Saul et
al., 1996, Jaakkola and Jordan, 1996], variational Bayes (VB) has become a
prominent method for approximate inference but, despite significant progress
in probabilistic programming, deriving VB updates remains cumbersome for many.
The main source of difficulty is the need to derive closed-form expressions
for the integrals. For example, consider observed data $\mathbf{y}$ and latent
vector
$\mbox{$\mbox{$\mathbf{z}$}$}=(\mbox{$\mbox{$\mathbf{z}$}$}_{1},\mbox{$\mbox{$\mathbf{z}$}$}_{2},\ldots,\mbox{$\mbox{$\mathbf{z}$}$}_{K})$
for a model $p(\mbox{$\mbox{$\mathbf{y}$}$},\mbox{$\mbox{$\mathbf{z}$}$})$,
and assume that we seek a mean-field posterior approximation:
$p(\mbox{$\mbox{$\mathbf{z}$}$}|\mbox{$\mbox{$\mathbf{y}$}$})\approx
q(\mbox{$\mbox{$\mathbf{z}$}$})=\prod_{i}q_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})$.
Then, a stationary point $q^{*}(\mbox{$\mbox{$\mathbf{z}$}$})$ of the evidence
lower-bound (ELBO) satisfies the following [Bishop, 2006, Eq. 10.9],

$\displaystyle\log
q_{i}^{*}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})=\mathbb{E}_{q^{*}_{\backslash
i}}\left[\log
p(\mbox{$\mbox{$\mathbf{y}$}$},\mbox{$\mbox{$\mathbf{z}$}$})\right]+\text{const.},$
(1)

where the integral in the right hand side marginalizes out the rest of the
variables (denoted by $\mbox{$\mbox{$\mathbf{z}$}$}_{\backslash i}$) using the
marginal $q^{*}_{\backslash i}(\mbox{$\mbox{$\mathbf{z}$}$}_{\backslash i})$;
see a proof in App. A.

For conjugate-exponential (CE) family models, the integrals have closed-form
expressions, but deriving them can still be cumbersome. It is not uncommon to
find papers where these derivations take many pages of work. With such long
derivations, it is easy to lose the main intuition and elegance of the
updates. Difficulty increases when non-conjugate factors are also present. For
such cases, despite using automatic differentiation, we still need to make
several difficult choices: Which optimizer to choose and are they compatible
with the conjugate updating? How to set the learning rate so that overall
updates converge? And, when using coordinate-wise updates, in what sequence
should we update the variables? Here, we present a recipe to derive the
updates that simplify these difficulties.

Our recipe consists of 3-steps based on an alternate way to express Eq. 1. We
assume an exponential-family form for
$q_{\text{\mbox{$\mbox{$\boldsymbol{\lambda}$}$}}_{i}}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})\propto\exp\left(\langle\text{\mbox{$\mbox{$\mathbf{T}$}$}}_{i}(\text{\mbox{$\mbox{$\mathbf{z}$}$}}_{i}),\text{\mbox{$\mbox{$\boldsymbol{\lambda}$}$}}_{i}\rangle\right)$
with natural parameter $\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{i}$,
sufficient statistics $\mbox{$\mbox{$\mathbf{T}$}$}_{i}(\cdot)$, and a
constant base-measure. We can then write Eq. 1 in a form shown below that uses
the pair
$(\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{i},\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{i})$
where
$\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{i}(\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{i})=\mathbb{E}_{q_{\text{\mbox{$\mbox{$\boldsymbol{\lambda}$}$}}_{i}}}[\mbox{$\mbox{$\mathbf{T}$}$}_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})]$
is the expectation parameter:

$\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{i}^{*}=\nabla_{\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{i}}\left.\mathbb{E}_{q_{\text{\mbox{$\mbox{$\boldsymbol{\lambda}$}$}}}}\left[\log
p(\mbox{$\mbox{$\mathbf{y}$}$},\mbox{$\mbox{$\mathbf{z}$}$})\right]\right|_{\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{i}=\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{i}^{*}}.$
(2)

Here, we marginalize with respect to
$q_{\text{\mbox{$\mbox{$\boldsymbol{\lambda}$}$}}}(\mbox{$\mbox{$\mathbf{z}$}$})$
where
$\mbox{$\mbox{$\boldsymbol{\lambda}$}$}=(\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{1},\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{2},\ldots,\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{K})$
and denote by
$(\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{i}^{*},\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{i}^{*})$
the natural and expectation parameter pair for the optimal
$q_{i}^{*}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})$. A proof is given in App. A by
using results from Khan and Lin [2017], along with an extension for non-
constant base-measures. The advantage of this expression is that we do not
have to think about the distributional forms of $q_{i}$ or any integrals.
Rather, we can simply look up the definition of
$(\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{i},\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{i})$;
see Table 1 for a list.

We give an example below. For notational ease, we will write
$\mathbb{E}_{q_{\text{\mbox{$\mbox{$\boldsymbol{\lambda}$}$}}}}[\cdot]$ as
$\mathbb{E}_{q}[\cdot]$.

###### Ex 1

(A simple mixture model) We consider two components
$p_{a}(\mbox{$\mbox{$\mathbf{y}$}$})$ and
$p_{b}(\mbox{$\mbox{$\mathbf{y}$}$})$ for an observation $\mathbf{y}$ where
mixture indicator $z\in\\{0,1\\}$ is sampled from a Bernoulli prior,

$\displaystyle
p(\mbox{$\mbox{$\mathbf{y}$}$}|z)=p_{a}(\mbox{$\mbox{$\mathbf{y}$}$})^{z}p_{b}(\mbox{$\mbox{$\mathbf{y}$}$})^{1-z},\text{
where }p(z)=\pi_{0}^{z}(1-\pi_{0})^{1-z}.$

The posterior over $z$ can be obtained by using Bayes’ rule,
$p(z=1|\mbox{$\mbox{$\mathbf{y}$}$})=\frac{\pi_{0}p_{a}(\text{\mbox{$\mbox{$\mathbf{y}$}$}})}{\pi_{0}p_{a}(\text{\mbox{$\mbox{$\mathbf{y}$}$}})+(1-\pi_{0})p_{b}(\text{\mbox{$\mbox{$\mathbf{y}$}$}})}$,
but let us suppose that we do not know the posterior and we want to recover it
from Eq. 2.

Denoting the posterior by $q(z)$, we can figure out its form by expanding

$\begin{split}\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{y}$}$},z)]&=\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{y}$}$}|z)+\log p(z)]\\\ &=\mathbb{E}_{q}\left[z\log
p_{a}(\mbox{$\mbox{$\mathbf{y}$}$})+(1-z)\log
p_{b}(\mbox{$\mbox{$\mathbf{y}$}$})+z\log\pi_{0}+(1-z)\log(1-\pi_{0})\right]\\\
&=\underbrace{\mathbb{E}_{q}(z)}_{=\mu}\underbrace{\log\frac{\pi_{0}p_{a}(\mbox{$\mbox{$\mathbf{y}$}$})}{(1-\pi_{0})p_{b}(\mbox{$\mbox{$\mathbf{y}$}$})}}_{\text{Coeff.
in front of
$\mu$}}+\underbrace{\log[(1-\pi_{0})p_{b}(\mbox{$\mbox{$\mathbf{y}$}$})]}_{\text{constant}}.\end{split}$
(3)

The last line suggests to choose $q(z)$ with expectation parameter
$\mu=\mathbb{E}_{q}(z)$. From Table 1, we find Bernoulli distribution to have
this expectation parameter, therefore we set
$q(z)\propto\pi_{1}^{z}(1-\pi_{1}^{1-z})$ where $\pi_{1}>0$ is the probability
of $z=1$ under $q$.

Having figured out the form of $q$, we can now use Eq. 2 to find its
parameter. For this we first look up from Table 1 the natural parameter:
$\lambda=\log\frac{\pi_{1}}{1-\pi_{1}}$ and set it according to Eq. 2. It is
clear from Eq. 3 that $\mathbb{E}_{q}\left[\log
p(\mbox{$\mbox{$\mathbf{y}$}$},z)\right]$ is linear in $\mu$, therefore the
gradient $\nabla_{\mu}\mathbb{E}_{q}\left[\log
p(\mbox{$\mbox{$\mathbf{y}$}$},z)\right]$ is simply the coefficient in front
of $\mu$. This gives us

$\lambda^{*}=\log\frac{\pi_{0}p_{a}(\mbox{$\mbox{$\mathbf{y}$}$})}{(1-\pi_{0})p_{b}(\mbox{$\mbox{$\mathbf{y}$}$})}\quad\implies\pi_{1}^{*}=\frac{\pi_{0}p_{a}(\mbox{$\mbox{$\mathbf{y}$}$})}{\pi_{0}p_{a}(\mbox{$\mbox{$\mathbf{y}$}$})+(1-\pi_{0})p_{b}(\mbox{$\mbox{$\mathbf{y}$}$})}.$
(4)

where in the second equality we rewrite the update in terms of $\pi_{1}$. This
recovers the posterior $p(z=1|\mbox{$\mbox{$\mathbf{y}$}$})$ obtained with the
Bayes’ rule.

The procedure is an instance of our 3-step recipe shown below where we
identify the posterior form by explicitly looking for linearity with respect
to $\boldsymbol{\mu}$ and use it to update $\boldsymbol{\lambda}$:

  1. 1.

Identify all $q_{i}$ at once, by looking for linearity of $\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{y}$}$},\mbox{$\mbox{$\mathbf{z}$}$})]$ w.r.t. some
$\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{i}$.

  2. 2.

Compute
$\nabla_{\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{i}}\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{y}$}$},\mbox{$\mbox{$\mathbf{z}$}$})]$ (hint: for CE,
just read-off the coefficient in front of
$\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{i}$).

  3. 3.

Update
$\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{i}\leftarrow(1-\rho_{i})\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{i}+\rho_{i}\nabla_{\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{i}}\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{y}$}$},\mbox{$\mbox{$\mathbf{z}$}$})]$ with a learning
rate $\rho_{i}$ (often set to 1).

For most cases, there is no need to derive any integrals: we just expand
$\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{y}$}$},\mbox{$\mbox{$\mathbf{z}$}$})]$ as a multi-
linear function of all $\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{i}$ at once and
look-up the rest from a table (similar to Table 1). We can then directly write
the update by simply “reading-off” the coefficient in front of
$\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{i}$, which is now a function of the
$\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{\backslash i}$ and observation
$\mathbf{y}$.

The steps differ fundamentally from the standard way of deriving VB where each
$q_{i}$ is identified _separately_. Such derivations can be very long and
tedious because we need to investigate each node separately by looking for a
known distributional form for each of them [Blei et al., 2017, App. A][Bishop,
2006, Sec. 10.2]. The distributional form is not used anyways because the
updates are ultimately written and implemented using the parameters of the
distributions. Our recipe makes the direct use of the
$(\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{i},\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{i})$
pair, and leads to an easier, faster, and shorter derivation where all nodes
are handled altogether.

The recipe is also more general and covers a wide-variety of cases, including
conjugate, non-conjugate, and even deterministic factors. This is due to Step
3 which uses the Bayesian learning rule (BLR) [Khan and Rue, 2021] that
contains many algorithms as special cases. Below, we give additional sub-steps
for Step 3 to derive many variants:

  1. (3a)

Derive Bayes-rule by using $\rho_{i}=1$; rewrite it in your parameter of
choice (see 1).

  2. (3b)

Derive Coordinate Ascent VI (CAVI) and variational message passing (VMP) by
setting $\rho_{i}=1$ and doing a coordinate-wise update (see 2 and 3).

  3. (3c)

Derive Stochastic VI (SVI) by setting $\rho_{i}=1$ for the local variables and
updating the global variable after a local update on any node $i$ (see 4).

  4. (3d)

Derive updates of deterministic nodes (and Laplace’s approximation) by
invoking the delta method to approximate the expectation (see 5).

  5. (3e)

In presence of the non-conjugate terms, just use the gradient and simplify it
using reparameterization trick whenever feasible [Khan and Rue, 2021] (see 6).

Another benefit is that, Steps 3a-3d do not only not require any integrals but
we do not even need any derivatives either; due to linearity we just need the
coefficients in front. With Step 3d, _even MAP estimates can be derived
without any derivatives_. Moreover, there is no need to use different
optimizers for conjugate and non-conjugate parts, and the learning rate and
update schedule need not be set carefully because the update converges under
fairly general conditions [Khan et al., 2016]. In contrast, coordinate-wise
update may not always converge and require the solution to be unique for each
coordinate [Paquet, 2014].

Table 1: Exponential-family used in this paper. For Gaussian-Wishart, we use two sets of variables: $\mbox{$\mbox{$\mathbf{z}$}$}_{1}$ is a real vector and $\mbox{$\mbox{$\mathbf{Z}$}$}_{2}$ is a positive-definite matrix. We do not give the expressions for the expectation parameters explicitly but these can be found in Wikipedia, with a more exhaustive list covering many other distributions. Name |  Distribution $q(\mbox{$\mbox{$\mathbf{z}$}$})$ | Expectation Param $\boldsymbol{\mu}$ | Natural param $\boldsymbol{\lambda}$
---|---|---|---
Bernoulli |  $\propto\pi^{z}(1-\pi^{1-z})$ | $\mathbb{E}_{q}(z)$ | $\log\frac{\pi}{1-\pi}$
Beta |  $\propto z^{\alpha-1}e^{-\beta z}$ | $\mathbb{E}_{q}\left[\begin{array}[]{l}\log z\\\ \log(1-z)\end{array}\right]$ | $\left[\begin{array}[]{l}\alpha-1\\\ \beta-1\end{array}\right]$
Gaussian |  $\propto e^{-\frac{1}{2}(\text{\mbox{$\mbox{$\mathbf{z}$}$}}-\text{\mbox{$\mbox{$\mathbf{m}$}$}})^{\top}\text{\mbox{$\mbox{$\mathbf{S}$}$}}(\text{\mbox{$\mbox{$\mathbf{z}$}$}}-\text{\mbox{$\mbox{$\mathbf{m}$}$}})}$ | $\mathbb{E}_{q}\left[\begin{array}[]{l}\mbox{$\mbox{$\mathbf{z}$}$}\\\ \mbox{$\mbox{$\mathbf{z}$}$}\mbox{$\mbox{$\mathbf{z}$}$}^{\top}\end{array}\right]$ | $\left[\begin{array}[]{l}\mbox{$\mbox{$\mathbf{S}$}$}\mbox{$\mbox{$\mathbf{m}$}$}\\\ -\mbox{$\frac{1}{2}$}\mbox{$\mbox{$\mathbf{S}$}$}\end{array}\right]$
Gaussian-Wishart |  $\propto|\mbox{$\mbox{$\mathbf{Z}$}$}_{2}|^{\frac{1}{2}}e^{-\frac{1}{2}(\text{\mbox{$\mbox{$\mathbf{z}$}$}}_{1}-\text{\mbox{$\mbox{$\mathbf{m}$}$}})^{\top}\gamma\text{\mbox{$\mbox{$\mathbf{Z}$}$}}_{2}(\text{\mbox{$\mbox{$\mathbf{z}$}$}}_{1}-\text{\mbox{$\mbox{$\mathbf{m}$}$}})}$ $|\mbox{$\mbox{$\mathbf{Z}$}$}_{2}|^{\frac{\nu-D-1}{2}}e^{-\frac{1}{2}\text{\mbox{Tr}}{(\text{\mbox{$\mbox{$\mathbf{W}$}$}}^{-1}\text{\mbox{$\mbox{$\mathbf{Z}$}$}}_{2})}}$ | $\mathbb{E}_{q}\left[\begin{array}[]{l}\log|\mbox{$\mbox{$\mathbf{Z}$}$}_{2}|\\\ \mbox{$\mbox{$\mathbf{Z}$}$}_{2}\\\ \mbox{$\mbox{$\mathbf{Z}$}$}_{2}\mbox{$\mbox{$\mathbf{z}$}$}_{1}\\\ \mbox{$\mbox{$\mathbf{z}$}$}_{1}^{\top}\mbox{$\mbox{$\mathbf{Z}$}$}_{2}\mbox{$\mbox{$\mathbf{z}$}$}_{1}\end{array}\right]$ | $\left[\begin{array}[]{l}\mbox{$\frac{1}{2}$}(\nu-D)\\\ -\mbox{$\frac{1}{2}$}(\mbox{$\mbox{$\mathbf{W}$}$}^{-1}+\gamma\mbox{$\mbox{$\mathbf{m}$}$}\mbox{$\mbox{$\mathbf{m}$}$}^{\top})\\\ \gamma\mbox{$\mbox{$\mathbf{m}$}$}\\\ -\mbox{$\frac{1}{2}$}\gamma\end{array}\right]$

##  2 Examples

Let us see another example with a slightly more complex model where a mean-
field approximation is required. We will see that, unlike standard derivation
where each factor requires a separate expansion of the log-joint, we can
identify both distributions simultaneously. Another advantage is that we do
not have to start with a mean-field structure, rather we can figure this out
by simply looking for linearity.

###### Ex 2

(A mixture-model with two levels) We will add an additional latent variable:
$\pi_{0}\rightarrow z_{i}\rightarrow\mbox{$\mbox{$\mathbf{y}$}$}_{i}$ where we
have multiple _vector_ observations $\mbox{$\mbox{$\mathbf{y}$}$}_{i}$, and
two sets of latent variables: $z_{i}$ and $\pi_{0}$. We assume the following
model with a beta prior over $\pi_{0}$,

$\begin{split}p(\mbox{$\mbox{$\mathbf{Y}$}$}|\mbox{$\mbox{$\mathbf{z}$}$})&=\prod_{i=1}^{N}p_{a}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})^{z_{i}}p_{b}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})^{1-z_{i}},\,\,p(z_{i}|\pi_{0})=\pi_{0}^{z_{i}}(1-\pi_{0})^{1-z_{i}},\,\,p(\pi_{0})\propto\pi_{0}^{\alpha_{0}-1}(1-\pi_{0})^{\beta_{0}-1}\end{split}$

Our goal is to estimate
$p(\mbox{$\mbox{$\mathbf{z}$}$},\pi_{0}|\mbox{$\mbox{$\mathbf{Y}$}$})$ but the
model is not a conjugate one, therefore we have to use a mean-field
approximation for tractability.

In the $1^{\text{st}}$ step we can figure out $q$ for both $z_{i}$ and
$\pi_{0}$ at once by using linearity in the expected log-joint; there is no
need to decide the factorization yet. We show this below where the second-line
is obtained by using Eq. 3,

$\displaystyle\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{Y}$}$}|\mbox{$\mbox{$\mathbf{z}$}$})p(\mbox{$\mbox{$\mathbf{z}$}$})p(\pi_{0})]=\mathbb{E}_{q}\left[\sum_{i=1}^{N}\left(\log
p(\mbox{$\mbox{$\mathbf{y}$}$}_{i}|z_{i})+\log p(z_{i})\right)+\log
p(\pi_{0})\right]$
$\displaystyle=\mathbb{E}_{q}\left[\sum_{i=1}^{N}\left(z_{i}\log\frac{\pi_{0}p_{a}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})}{(1-\pi_{0})p_{b}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})}+\log[(1-\pi_{0})p_{b}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})]\right)+\log\left(\pi_{0}^{(\alpha_{0}-1)}(1-\pi_{0})^{(\beta_{0}-1)}\right)\right]$
$\displaystyle=\sum_{i=1}^{N}\underbrace{\mathbb{E}_{q}(z_{i})}_{=\mu_{i}}\Bigg{(}\underbrace{\mathbb{E}_{q}{\left[\begin{array}[]{c}\log\pi_{0}\\\
\log(1-\pi_{0})\end{array}\right]}}_{=\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{0}}+\left[\begin{array}[]{c}\log
p_{a}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})\\\ \log
p_{b}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})\end{array}\right]\Bigg{)}^{\top}\left[\begin{array}[]{c}+1\\\
-1\end{array}\right]+\sum_{i}\log p_{b}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})$
(11)
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad+{\underbrace{\mathbb{E}_{q}{\left[\begin{array}[]{c}\log\pi_{0}\\\
\log(1-\pi_{0})\end{array}\right]}}_{=\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{0}}}^{\top}\left[\begin{array}[]{c}\alpha_{0}-1\\\
N+\beta_{0}-1\end{array}\right],$ (16)

where in the last step we assumed that
$q(\mbox{$\mbox{$\mathbf{z}$}$},\pi_{0})=\prod_{i}q(z_{i})q(\pi_{0})$ to get
linearity with respect to both $\mu_{1}$ and
$\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{0}$. Using Table 1, we have $q(z_{i})$ as
Bernoulli and $q(\pi_{0})$ as Beta.

After this, the $2^{\text{nd}}$ step is straightforward to simply read the
coefficient in front of $\mu_{i}$ and $\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{0}$
respectively, and use them in the $3^{\text{rd}}$ step with $\rho_{i}=1$. This
is shown below in both natural parameters and posterior parameters of
Bernoulli (denoted by $\pi_{i,1}$) and Beta distributions (denoted by
$(\alpha_{1},\beta_{1})$) obtained by using Table 1:

$\displaystyle\lambda_{i}$
$\displaystyle\leftarrow\Bigg{(}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{0}+\left[\begin{array}[]{c}\log
p_{a}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})\\\ \log
p_{b}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})\end{array}\right]\Bigg{)}^{\top}\left[\begin{array}[]{c}1\\\
-1\end{array}\right]\implies\pi_{i,1}\leftarrow\frac{e^{\mathbb{E}_{q}[\log\pi_{0}]}p_{a}(\mbox{$\mbox{$\mathbf{y}$}$})}{e^{\mathbb{E}_{q}[\log\pi_{0}]}p_{a}(\mbox{$\mbox{$\mathbf{y}$}$})+e^{\mathbb{E}_{q}[\log(1-\pi_{0})]}p_{b}(\mbox{$\mbox{$\mathbf{y}$}$})}$
(21) $\displaystyle\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{0}$
$\displaystyle\leftarrow\sum_{i=1}^{N}\mu_{i}\left[\begin{array}[]{c}+1\\\
-1\end{array}\right]+\left[\begin{array}[]{c}\alpha_{0}-1\\\
N+\beta_{0}-1\end{array}\right]\implies\left[\begin{array}[]{c}\alpha_{1}\\\
\beta_{1}\end{array}\right]\leftarrow\left[\begin{array}[]{c}\alpha_{0}+\sum_{i}\mathbb{E}_{q}(z_{i})\\\
\beta_{0}+\sum_{i}(1-\mathbb{E}_{q}(z_{i}))\end{array}\right]$ (30)

The recipe identifies the distributions for all nodes at once by looking for
linearity; see the last line in Eq. 16. In addition, it reuses the already-
known integrals by looking-up the required expectations. For instance, in the
example above, we need to find $\mathbb{E}_{q}(z_{i})$,
$\mathbb{E}_{q}(\log\pi_{0})$, and $\mathbb{E}_{q}(\log(1-\pi_{0}))$, whose
expressions can be found in Wikipedia pages for Bernoulli and Beta
distributions.

We are now ready to discuss the most-commonly used example for VB, which is
the Gaussian mixture model. Unlike the standard derivation covered in [Bishop,
2006], we can derive all the updates with just a single expansion of the
expected log-joint. The example clearly demonstrates the usefulness of the
recipe in simplifying the derivation.

###### Ex 3

(Gaussian mixture model) The model uses a Gaussian likelihood for the two
components and a Gaussian-Wishart prior on their means and covariances,

$\displaystyle p_{a}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})$
$\displaystyle=\mbox{${\cal
N}$}(\mbox{$\mbox{$\mathbf{y}$}$}_{i}|\mbox{$\mbox{$\mathbf{m}$}$}_{a},\mbox{$\mbox{$\mathbf{S}$}$}_{a}^{-1}),\quad\quad$
$\displaystyle
p(\mbox{$\mbox{$\mathbf{m}$}$}_{a},\mbox{$\mbox{$\mathbf{S}$}$}_{a})=\mbox{${\cal
N}$}\left(\mbox{$\mbox{$\mathbf{m}$}$}_{a}|0,(\gamma_{0}\mbox{$\mbox{$\mathbf{S}$}$}_{a})^{-1}\right)\mathcal{W}(\mbox{$\mbox{$\mathbf{S}$}$}_{a}|\mbox{$\mbox{$\mathbf{W}$}$}_{0},\nu_{0})$
$\displaystyle p_{b}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})$
$\displaystyle=\mbox{${\cal
N}$}(\mbox{$\mbox{$\mathbf{y}$}$}_{i}|\mbox{$\mbox{$\mathbf{m}$}$}_{b},\mbox{$\mbox{$\mathbf{S}$}$}_{b}^{-1}),\quad\quad$
$\displaystyle
p(\mbox{$\mbox{$\mathbf{m}$}$}_{b},\mbox{$\mbox{$\mathbf{S}$}$}_{b})=\mbox{${\cal
N}$}\left(\mbox{$\mbox{$\mathbf{m}$}$}_{b}|0,(\gamma_{0}\mbox{$\mbox{$\mathbf{S}$}$}_{b})^{-1}\right)\mathcal{W}(\mbox{$\mbox{$\mathbf{S}$}$}_{b}|\mbox{$\mbox{$\mathbf{W}$}$}_{0},\nu_{0})$

for scalars $\gamma_{0}>0$ and $\nu_{0}>D-1$, and positive-definite matrix
$\mbox{$\mbox{$\mathbf{W}$}$}_{0}$. For this, we simply need to make two
changes in Eq. 16. First, we replace $\log
p_{a}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})$ and $\log
p_{b}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})$ by their expected values. We show one
of them below in terms of the required expectations,

$\displaystyle\mathbb{E}_{q}[\log
p_{a}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})]=\mbox{$\frac{1}{2}$}\mathbb{E}_{q}[\log|\mbox{$\mbox{$\mathbf{S}$}$}_{a}|]-\mbox{$\frac{1}{2}$}\mathbb{E}_{q}\left[(\mbox{$\mbox{$\mathbf{y}$}$}_{i}-\mbox{$\mbox{$\mathbf{m}$}$}_{a})^{\top}\mbox{$\mbox{$\mathbf{S}$}$}_{a}(\mbox{$\mbox{$\mathbf{y}$}$}_{i}-\mbox{$\mbox{$\mathbf{m}$}$}_{a})\right]+c$
$\displaystyle\quad\quad=\mbox{$\frac{1}{2}$}\underbrace{\mathbb{E}_{q}[\log|\mbox{$\mbox{$\mathbf{S}$}$}_{a}|]}_{\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{p_{a},1}}-\mbox{$\frac{1}{2}$}\mbox{Tr}\big{(}\mbox{$\mbox{$\mathbf{y}$}$}_{i}\mbox{$\mbox{$\mathbf{y}$}$}_{i}^{\top}\underbrace{\mathbb{E}_{q}[\mbox{$\mbox{$\mathbf{S}$}$}_{a}]}_{\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{p_{a},2}}\big{)}+\mbox{$\mbox{$\mathbf{y}$}$}_{i}^{\top}\underbrace{\mathbb{E}_{q}[\mbox{$\mbox{$\mathbf{S}$}$}_{a}\mbox{$\mbox{$\mathbf{m}$}$}_{a}]}_{\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{p_{a},3}}-\mbox{$\frac{1}{2}$}\underbrace{\mathbb{E}_{q}[\mbox{$\mbox{$\mathbf{m}$}$}_{a}^{\top}\mbox{$\mbox{$\mathbf{S}$}$}_{a}\mbox{$\mbox{$\mathbf{m}$}$}_{a}]}_{\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{p_{a},4}}+c.$

The term is linear in terms of the expectation-parameter
$\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{a},1:4}$ which corresponds to the
Gaussian-Wishart distribution (Table 1); see the Wikipedia pages for their
expressions. Second, we need to add the contribution of the prior, which is
also linear in $\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{a},1:4}$,

$\displaystyle\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{m}$}$}_{a},\mbox{$\mbox{$\mathbf{S}$}$}_{a})]=\mbox{$\frac{1}{2}$}(\nu_{0}-D)\mathbb{E}_{q}[\log|\mbox{$\mbox{$\mathbf{S}$}$}_{a}|]-\mbox{$\frac{1}{2}$}\mbox{Tr}\big{(}\mbox{$\mbox{$\mathbf{W}$}$}_{0}^{-1}\mathbb{E}_{q}[\mbox{$\mbox{$\mathbf{S}$}$}_{a}]\big{)}-\mbox{$\frac{1}{2}$}\mathbb{E}_{q}[\mbox{$\mbox{$\mathbf{m}$}$}_{a}^{\top}(\gamma_{0}\mbox{$\mbox{$\mathbf{S}$}$}_{a})\mbox{$\mbox{$\mathbf{m}$}$}_{a}]+c$
$\displaystyle\quad\quad=\mbox{$\frac{1}{2}$}(\nu_{0}-D)\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{a},1}-\mbox{$\frac{1}{2}$}\mbox{Tr}\big{(}\mbox{$\mbox{$\mathbf{W}$}$}_{0}^{-1}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{a},3}\big{)}-\mbox{$\frac{1}{2}$}\gamma_{0}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{a},4}+c.$

Then, we expand the expected log-joint by using model definition in the first
line, then by plugging Eq. 16 to get the second line, and do rearrangement
afterward to get linearity,

$\displaystyle\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{Y}$}$}|\mbox{$\mbox{$\mathbf{z}$}$},\mbox{$\mbox{$\mathbf{m}$}$}_{a:b},\mbox{$\mbox{$\mathbf{S}$}$}_{a:b})p(\mbox{$\mbox{$\mathbf{z}$}$})p(\pi_{0})p(\mbox{$\mbox{$\mathbf{m}$}$}_{a:b},\mbox{$\mbox{$\mathbf{S}$}$}_{a:b})]$
$\displaystyle=\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{Y}$}$}|\mbox{$\mbox{$\mathbf{z}$}$},\mbox{$\mbox{$\mathbf{m}$}$}_{a:b},\mbox{$\mbox{$\mathbf{S}$}$}_{a:b})p(\mbox{$\mbox{$\mathbf{z}$}$})p(\pi_{0})]+\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{m}$}$}_{a:b},\mbox{$\mbox{$\mathbf{S}$}$}_{a:b})]$
$\displaystyle=\sum_{i=1}^{N}\mu_{i}\Bigg{(}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{0}+\left[\begin{array}[]{c}\mathbb{E}_{q}(\log
p_{a}(\mbox{$\mbox{$\mathbf{y}$}$}_{i}))\\\ \mathbb{E}_{q}(\log
p_{b}(\mbox{$\mbox{$\mathbf{y}$}$}_{i}))\end{array}\right]\Bigg{)}^{\top}\left[\begin{array}[]{c}+1\\\
-1\end{array}\right]+\sum_{i}\mathbb{E}_{q}[\log
p_{b}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})]+\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{0}^{\top}\left[\begin{array}[]{c}\alpha_{0}-1\\\
N+\beta_{0}-1\end{array}\right]$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{m}$}$}_{a},\mbox{$\mbox{$\mathbf{S}$}$}_{a})]+\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{m}$}$}_{b},\mbox{$\mbox{$\mathbf{S}$}$}_{b})]$
$\displaystyle=\sum_{i=1}^{N}\mu_{i}\Bigg{(}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{0}+\left[\begin{array}[]{c}\mbox{$\frac{1}{2}$}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{a},1}-\mbox{$\frac{1}{2}$}\mbox{Tr}\left(\mbox{$\mbox{$\mathbf{y}$}$}_{i}\mbox{$\mbox{$\mathbf{y}$}$}_{i}^{\top}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{a},2}\right)+\mbox{$\mbox{$\mathbf{y}$}$}_{i}^{\top}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{a},3}-\mbox{$\frac{1}{2}$}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{a},4}\\\
\mbox{$\frac{1}{2}$}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{b},1}-\mbox{$\frac{1}{2}$}\mbox{Tr}\left(\mbox{$\mbox{$\mathbf{y}$}$}_{i}\mbox{$\mbox{$\mathbf{y}$}$}_{i}^{\top}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{b},2}\right)+\mbox{$\mbox{$\mathbf{y}$}$}_{i}^{\top}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{b},3}-\mbox{$\frac{1}{2}$}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{a},4}\end{array}\right]\Bigg{)}^{\top}\left[\begin{array}[]{c}+1\\\
-1\end{array}\right]$
$\displaystyle\quad\quad\quad\quad+\sum_{i}\left(\mbox{$\frac{1}{2}$}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{b},1}-\mbox{$\frac{1}{2}$}\mbox{Tr}\big{[}\mbox{$\mbox{$\mathbf{y}$}$}_{i}\mbox{$\mbox{$\mathbf{y}$}$}_{i}^{\top}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{b},2}\big{]}+\mbox{$\mbox{$\mathbf{y}$}$}_{i}^{\top}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{b},3}-\mbox{$\frac{1}{2}$}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{b},4}\right)+\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{0}^{\top}\left[\begin{array}[]{c}\alpha_{0}-1\\\
N+\beta_{0}-1\end{array}\right]$
$\displaystyle\quad\quad\quad\quad+\mbox{$\frac{1}{2}$}(\nu_{0}-D)\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{a},1}-\mbox{$\frac{1}{2}$}\mbox{Tr}\left(\mbox{$\mbox{$\mathbf{W}$}$}_{0}^{-1}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{a},2}\right)-\mbox{$\frac{1}{2}$}\gamma_{0}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{a},4}$
$\displaystyle\quad\quad\quad\quad+\mbox{$\frac{1}{2}$}(\nu_{0}-D)\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{b},1}-\mbox{$\frac{1}{2}$}\mbox{Tr}\left(\mbox{$\mbox{$\mathbf{W}$}$}_{0}^{-1}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{b},2}\right)-\mbox{$\frac{1}{2}$}\gamma_{0}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{b},4}+\text{const.}$

The last equation looks complicated but it is linear with respect to each
$\mu_{i},\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{0},\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{a}},\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{b}}$,
therefore we can simply look up the coefficient in front to write the update.
The updates of $q(\pi_{0})$ remain same as before because the coefficient is
unchanged. For $q(z_{i})$, the coefficient takes an expectation over $\log
p_{a}$ and $\log p_{b}$, giving us the following,

$\displaystyle\pi_{i,1}$
$\displaystyle\leftarrow\frac{e^{\mathbb{E}_{q}[\log\pi_{0}]+\mathbb{E}_{q}[\log
p_{a}(\text{\mbox{$\mbox{$\mathbf{y}$}$}})]}}{e^{\mathbb{E}_{q}[\log\pi_{0}]+\mathbb{E}_{q}[\log
p_{a}(\text{\mbox{$\mbox{$\mathbf{y}$}$}})]}+e^{\mathbb{E}_{q}[\log(1-\pi_{0})]+\mathbb{E}_{q}[\log
p_{b}(\text{\mbox{$\mbox{$\mathbf{y}$}$}})]}}$ (31)

The update for the mean and covariance of the component is also obtained in a
straightforward manner. For
$q(\mbox{$\mbox{$\mathbf{m}$}$}_{a},\mbox{$\mbox{$\mathbf{S}$}$}_{a})$, for
instance, we set its 4 natural parameters (given in Table 1) to the
coefficients in front of $\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{a},1}$ to
$\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{p_{a},4}$ respectively, and rearrange to
get

$\displaystyle\left[\begin{array}[]{l}\mbox{$\frac{1}{2}$}(\nu_{a}-D)\\\
-\mbox{$\frac{1}{2}$}(\mbox{$\mbox{$\mathbf{W}$}$}_{a}^{-1}+\gamma_{a}\mbox{$\mbox{$\mathbf{m}$}$}_{a}\mbox{$\mbox{$\mathbf{m}$}$}_{a}^{\top})\\\
\gamma_{a}\mbox{$\mbox{$\mathbf{m}$}$}_{a}\\\
-\mbox{$\frac{1}{2}$}\gamma_{a}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{l}\mbox{$\frac{1}{2}$}\sum_{i}\mu_{i}+\mbox{$\frac{1}{2}$}(\nu_{0}-D)\\\
-\mbox{$\frac{1}{2}$}\sum_{i}\mu_{i}\mbox{$\mbox{$\mathbf{y}$}$}_{i}\mbox{$\mbox{$\mathbf{y}$}$}_{i}^{\top}-\mbox{$\frac{1}{2}$}\mbox{$\mbox{$\mathbf{W}$}$}_{o}^{-1}\\\
\sum_{i}\mu_{i}\mbox{$\mbox{$\mathbf{y}$}$}_{i}\\\
-\mbox{$\frac{1}{2}$}\sum_{i}\mu_{i}-\mbox{$\frac{1}{2}$}\gamma_{0}\end{array}\right]$
(40) $\displaystyle\implies\left[\begin{array}[]{l}\nu_{a}\\\
\mbox{$\mbox{$\mathbf{W}$}$}_{a}^{-1}\\\ \mbox{$\mbox{$\mathbf{m}$}$}_{a}\\\
\gamma_{a}\end{array}\right]$
$\displaystyle=\left[\begin{array}[]{l}\sum_{i}\mu_{i}+\nu_{0}\\\
\sum_{i}\mu_{i}\mbox{$\mbox{$\mathbf{y}$}$}_{i}\mbox{$\mbox{$\mathbf{y}$}$}_{i}^{\top}-\frac{(\sum_{i}\mu_{i}\text{\mbox{$\mbox{$\mathbf{y}$}$}}_{i})(\sum_{i}\mu_{i}\text{\mbox{$\mbox{$\mathbf{y}$}$}}_{i})^{\top}}{\sum_{i}\mu_{i}+\gamma_{0}}+\mbox{$\mbox{$\mathbf{W}$}$}_{0}^{-1}\\\
\frac{1}{\gamma_{a}}\sum_{i}\mu_{i}\mbox{$\mbox{$\mathbf{y}$}$}_{i}\\\
\sum_{i}\mu_{i}+\gamma_{0}\end{array}\right]$ (49)

We encourage the reader to verify that these updates are same as ones derived
in [Bishop, 2006, Eqs 10.49, 10.58, 10.60-10.63], but unlike that derivation,
all updates here are derived by using just one expansion of expected log-
joint, which leads to a shorter and much faster derivation.

We now discuss Sub-steps (3b)-(3f) and show that, by using the BLR in Step 3,
we can derive many other types of algorithms from the VB updates. It makes
sense to start from a more Bayesian approach and make approximation when
required, rather than doing it the other way. Using the BLR in Step 3 does
exactly this. For example, as we saw in the previous examples, we can write a
CAVI or VMP style update as special cases by choosing $\rho_{i}=1$, and doing
coordinate-wise updates. But the same approach can be used to derive SVI
[Hoffman et al., 2013], described below for completeness.

###### Ex 4

(SVI updates for the mixture model with two levels) Instead of a coordinate-
wise update that goes through all $z_{i}$, we can pick randomly just one
$z_{i}$, do the update with $\rho_{i}=1$ using Eq. 4, and after this we update
the global variable $\pi_{0}$ with $\rho_{0}<1$,

$\displaystyle\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{0}$
$\displaystyle\leftarrow(1-\rho_{0})\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{0}+\rho_{0}\sum_{i=1}^{N}\mu_{i}\left[\begin{array}[]{c}+1\\\
-1\end{array}\right]+\left[\begin{array}[]{c}\alpha_{0}-1\\\
N+\beta_{0}-1\end{array}\right]$ (54)
$\displaystyle\implies\left[\begin{array}[]{c}\alpha_{1}\\\
\beta_{1}\end{array}\right]\leftarrow(1-\rho_{0})\left[\begin{array}[]{c}\alpha_{1}\\\
\beta_{1}\end{array}\right]+\rho_{0}\left[\begin{array}[]{c}\alpha_{0}+\sum_{i}\mathbb{E}_{q}(z_{i})\\\
\beta_{0}+N-\sum_{i}\mathbb{E}_{q}(z_{i})\end{array}\right]$ (61)

Updating $\pi_{0}$ after every $z_{i}$ update can speed up convergence
[Hoffman et al., 2013].

Now, we give another example on matrix factorization where we show how to
derive updates for deterministic nodes and also for expectation maximization.
We derive the well-known alternating least-squares and probabilistic PCA
[Tipping and Bishop, 1999]. Both of these are obtained by using the delta
method where the expectation parameter is approximated by using its mean:
$\mathbb{E}_{q}(\mbox{$\mbox{$\mathbf{z}$}$}\mbox{$\mbox{$\mathbf{z}$}$}^{\top})\approx\mbox{$\mbox{$\mathbf{m}$}$}\mbox{$\mbox{$\mathbf{m}$}$}^{\top}$,
where $\mathbf{m}$ is the mean of $q$. More details on such delta method is
given in Khan and Rue [2021].

###### Ex 5

Matrix factorization and latent factor models aim to fit data matrix
$\mathbf{Y}$ of size $N\times D$ by using two sets of factors $\mathbf{U}$ and
$\mathbf{V}$ respectively of sizes $N\times K$ and $D\times K$ where $K$ is
the number of factors. The likelihood and the prior are given by,

$p(\mbox{$\mbox{$\mathbf{Y}$}$}|\mbox{$\mbox{$\mathbf{U}$}$},\mbox{$\mbox{$\mathbf{V}$}$})=\prod_{i=1}^{N}\prod_{j=1}^{D}\mbox{${\cal
N}$}(y_{ij}|\mbox{$\mbox{$\mathbf{u}$}$}_{i}^{\top}\mbox{$\mbox{$\mathbf{v}$}$}_{j},\mbox{$\mbox{$\mathbf{I}$}$}),\quad
p(\mbox{$\mbox{$\mathbf{u}$}$}_{i})=\mbox{${\cal
N}$}(\mbox{$\mbox{$\mathbf{u}$}$}_{i}|0,\mbox{$\mbox{$\mathbf{I}$}$}/\delta_{u}),\quad
p(\mbox{$\mbox{$\mathbf{v}$}$}_{j})=\mbox{${\cal
N}$}(\mbox{$\mbox{$\mathbf{v}$}$}_{j}|0,\mbox{$\mbox{$\mathbf{I}$}$}/\delta_{v})$

where $y_{ij}$ is the $ij$’th entry, $\mbox{$\mbox{$\mathbf{u}$}$}_{i}$ is the
$i$’th row of $\mathbf{U}$, and $\mbox{$\mbox{$\mathbf{v}$}$}_{j}$ is the
$j$’th row of $\mathbf{V}$. The most popular procedure is to use the
alternating least-squares (ALS) procedure, but there is also expectation
maximization [Tipping and Bishop, 1999] and VMP Paquet [2014]. These can all
be derived by using our recipe for VB, as we show now.

Following the $1^{\text{st}}$ step, we expand the expected log-joint and look
for linearity,

$\displaystyle\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{Y}$}$},\mbox{$\mbox{$\mathbf{U}$}$},\mbox{$\mbox{$\mathbf{V}$}$})]=\mathbb{E}_{q}\left[\sum_{i=1}^{N}\left(\sum_{j=1}^{D}-\frac{1}{2}(y_{ij}-\mbox{$\mbox{$\mathbf{u}$}$}_{i}^{\top}\mbox{$\mbox{$\mathbf{v}$}$}_{j})^{2}-\frac{\delta_{v}}{2}\mbox{$\mbox{$\mathbf{v}$}$}_{j}^{\top}\mbox{$\mbox{$\mathbf{v}$}$}_{j}\right)-\frac{\delta_{u}}{2}\mbox{$\mbox{$\mathbf{u}$}$}_{i}^{\top}\mbox{$\mbox{$\mathbf{u}$}$}_{i}\right]+\textrm{const}$
$\displaystyle=\sum_{i=1}^{N}\Bigg{(}\sum_{j=1}^{D}-\frac{1}{2}\Bigg{(}\mbox{Tr}\Bigg{(}\underbrace{\mathbb{E}_{q}(\mbox{$\mbox{$\mathbf{u}$}$}_{i}\mbox{$\mbox{$\mathbf{u}$}$}_{i}^{\top})}_{=\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{u_{i}}^{(2)}}\underbrace{\mathbb{E}_{q}(\mbox{$\mbox{$\mathbf{v}$}$}_{j}\mbox{$\mbox{$\mathbf{v}$}$}_{j}^{\top})}_{=\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{v_{j}}^{(2)}}\Bigg{)}-2y_{ij}{\underbrace{\mathbb{E}_{q}(\mbox{$\mbox{$\mathbf{u}$}$}_{i})}_{=\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{u_{i}}^{(1)}}}^{\top}\underbrace{\mathbb{E}_{q}(\mbox{$\mbox{$\mathbf{v}$}$}_{j})}_{\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{v_{j}}^{(1)}}\Bigg{)}$
$\displaystyle\qquad\qquad\qquad\qquad-\frac{\delta_{v}}{2}\mbox{Tr}\underbrace{\mathbb{E}_{q}\left(\mbox{$\mbox{$\mathbf{v}$}$}_{j}\mbox{$\mbox{$\mathbf{v}$}$}_{j}^{\top}\right)}_{=\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{v_{j}}^{(2)}}\Bigg{)}-\frac{\delta_{u}}{2}\mbox{Tr}\underbrace{\mathbb{E}_{q}\left(\mbox{$\mbox{$\mathbf{u}$}$}_{i}\mbox{$\mbox{$\mathbf{u}$}$}_{i}^{\top}\right)}_{=\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{u_{i}}^{(2)}}+\textrm{const}$

Here, we use mean-field approximation for all
$\mbox{$\mbox{$\mathbf{u}$}$}_{i}$ and $\mbox{$\mbox{$\mathbf{v}$}$}_{j}$
because we want linearity with respect to all. The sufficient statistics
correspond to Gaussian (third row in Table 1).

The $2^{\text{nd}}$ step is to simply read off as the coefficients in the
front, and, by using the natural parameters given in Table 1, we can directly
write the update from the $3^{\text{rd}}$ step. The first two updates below
use the natural parameters denoted by
$(\mbox{$\mbox{$\mathbf{S}$}$}_{u_{i}}\mbox{$\mbox{$\mathbf{m}$}$}_{u_{i}},-\mbox{$\mbox{$\mathbf{S}$}$}_{u_{i}}/2)$
corresponding to
$(\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{u_{i}}^{(1)},\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{u_{i}}^{(2)})$,
while the next two updates do the same for the variable
$\mbox{$\mbox{$\mathbf{v}$}$}_{j}$,

$\displaystyle\mbox{$\mbox{$\mathbf{S}$}$}_{u_{i}}\mbox{$\mbox{$\mathbf{m}$}$}_{u_{i}}$
$\displaystyle\leftarrow(1-\rho_{i})\mbox{$\mbox{$\mathbf{S}$}$}_{u_{i}}\mbox{$\mbox{$\mathbf{m}$}$}_{u_{i}}+\rho_{i}{\sum_{j=1}^{D}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{v_{j}}^{(1)}y_{ij}}$
$\displaystyle\mbox{$\mbox{$\mathbf{S}$}$}_{u_{i}}$
$\displaystyle\leftarrow(1-\rho_{i})\mbox{$\mbox{$\mathbf{S}$}$}_{u_{i}}+\rho_{i}{\Bigg{(}\sum_{j=1}^{D}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{v_{j}}^{(2)}+\delta_{u}\mbox{$\mbox{$\mathbf{I}$}$}_{K}\Bigg{)}}$
$\displaystyle\mbox{$\mbox{$\mathbf{S}$}$}_{v_{j}}\mbox{$\mbox{$\mathbf{m}$}$}_{v_{j}}$
$\displaystyle\leftarrow(1-\rho_{j})\mbox{$\mbox{$\mathbf{S}$}$}_{v_{j}}\mbox{$\mbox{$\mathbf{m}$}$}_{v_{j}}+\rho_{j}{\sum_{i=1}^{N}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{u_{i}}^{(1)}y_{ij}}$
$\displaystyle\mbox{$\mbox{$\mathbf{S}$}$}_{v_{j}}$
$\displaystyle\leftarrow(1-\rho_{j})\mbox{$\mbox{$\mathbf{S}$}$}_{v_{j}}+\rho_{j}\Bigg{(}{\sum_{i=1}^{N}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{u_{i}}^{(2)}+\delta_{v}\mbox{$\mbox{$\mathbf{I}$}$}_{K}}\Bigg{)}.$

We can then specialize these updates to derive ALS, EM, VB etc.

  1. 1.

Setting $\rho_{i}=\rho_{j}=1$ and updating posteriors of $\mathbf{U}$ and
$\mathbf{V}$ in alternate iterations (that is coordinate wise updates), we get
the VMP update.

$\begin{split}\mbox{$\mbox{$\mathbf{S}$}$}_{u_{i}}\mbox{$\mbox{$\mathbf{m}$}$}_{u_{i}}&\leftarrow{\sum_{j=1}^{D}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{v_{j}}^{(1)}y_{ij}},\qquad\qquad\mbox{$\mbox{$\mathbf{S}$}$}_{u_{i}}\leftarrow\sum_{j=1}^{D}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{v_{j}}^{(2)}+\delta_{u}\mbox{$\mbox{$\mathbf{I}$}$}_{K}\\\
\mbox{$\mbox{$\mathbf{S}$}$}_{v_{j}}\mbox{$\mbox{$\mathbf{m}$}$}_{v_{j}}&\leftarrow{\sum_{i=1}^{N}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{u_{i}}^{(1)}y_{ij}},\qquad\qquad\mbox{$\mbox{$\mathbf{S}$}$}_{v_{j}}\leftarrow{\sum_{i=1}^{N}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{u_{i}}^{(2)}+\delta_{v}\mbox{$\mbox{$\mathbf{I}$}$}_{K}}.\end{split}$
(62)

  2. 2.

We can get the probabilistic PCA updates [Tipping and Bishop, 1999], where we
only compute the posterior wrt $\mathbf{U}$ and assume $\mathbf{V}$ to be
deterministic. We can do this by making two more additional approximations,

    1. (a)

denote $\mbox{$\mbox{$\mathbf{m}$}$}_{v_{j}}$ by
$\hat{\mbox{$\mbox{$\mathbf{v}$}$}}_{j}$, and

    2. (b)

use the delta approximation for
$\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{v_{j}}^{(2)}=\mathbb{E}_{q}[\mbox{$\mbox{$\mathbf{v}$}$}_{j}\mbox{$\mbox{$\mathbf{v}$}$}_{j}^{\top}]\,\,\approx\,\,\mathbb{E}_{q}[\mbox{$\mbox{$\mathbf{v}$}$}_{j}]\mathbb{E}_{q}[\mbox{$\mbox{$\mathbf{v}$}$}_{j}]^{\top}=\hat{\mbox{$\mbox{$\mathbf{v}$}$}}_{j}\hat{\mbox{$\mbox{$\mathbf{v}$}$}}_{j}^{\top}$.

With these, the update in Eq. 62 reduces to

$\begin{split}\mbox{$\mbox{$\mathbf{S}$}$}_{u_{i}}\mbox{$\mbox{$\mathbf{m}$}$}_{u_{i}}&\leftarrow\sum_{j=1}^{D}\hat{\mbox{$\mbox{$\mathbf{v}$}$}}_{j}y_{ij},\quad\mbox{$\mbox{$\mathbf{S}$}$}_{u_{i}}\leftarrow\sum_{j=1}^{D}\hat{\mbox{$\mbox{$\mathbf{v}$}$}}_{j}\hat{\mbox{$\mbox{$\mathbf{v}$}$}}_{j}^{\top}+\delta_{u}\mbox{$\mbox{$\mathbf{I}$}$}_{K}\\\
\hat{\mbox{$\mbox{$\mathbf{v}$}$}}_{j}&\leftarrow\left(\sum_{i=1}^{N}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{u_{i}}^{(2)}+\delta_{v}\mbox{$\mbox{$\mathbf{I}$}$}_{K}\right)^{-1}{\sum_{i=1}^{N}\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{u_{i}}^{(1)}\,y_{ij}}.\end{split}$
(63)

These are equivalent to those derived in Bishop [2006, Eqs. 25-17].

  3. 3.

Similarly, we get the ALS updates by further

    1. (a)

denoting $\mbox{$\mbox{$\mathbf{m}$}$}_{u_{i}}$ by
$\hat{\mbox{$\mbox{$\mathbf{u}$}$}}_{i}$, and

    2. (b)

adding the delta approximation:
$\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{u_{i}}^{(2)}=\mathbb{E}_{q}[\mbox{$\mbox{$\mathbf{u}$}$}_{i}\mbox{$\mbox{$\mathbf{u}$}$}_{i}^{\top}]\,\,\approx\,\,\mathbb{E}_{q}[\mbox{$\mbox{$\mathbf{u}$}$}_{i}]\mathbb{E}_{q}[\mbox{$\mbox{$\mathbf{u}$}$}_{i}]^{\top}=\hat{\mbox{$\mbox{$\mathbf{u}$}$}}_{i}\hat{\mbox{$\mbox{$\mathbf{u}$}$}}_{i}^{\top}$.

With these, Eq. 63 reduces to the following ALS scheme,

$\begin{split}\hat{\mbox{$\mbox{$\mathbf{u}$}$}}_{i}&\leftarrow\Bigg{(}\sum_{j=1}^{D}\hat{\mbox{$\mbox{$\mathbf{v}$}$}}_{j}\hat{\mbox{$\mbox{$\mathbf{v}$}$}}_{j}^{\top}+\delta_{u}\mbox{$\mbox{$\mathbf{I}$}$}_{K}\Bigg{)}^{-1}\sum_{j=1}^{D}\hat{\mbox{$\mbox{$\mathbf{v}$}$}}_{j}y_{ij},\\\
\hat{\mbox{$\mbox{$\mathbf{v}$}$}}_{j}&\leftarrow\Bigg{(}\sum_{i=1}^{N}\hat{\mbox{$\mbox{$\mathbf{u}$}$}}_{i}\hat{\mbox{$\mbox{$\mathbf{u}$}$}}_{i}^{\top}+\delta_{v}\mbox{$\mbox{$\mathbf{I}$}$}_{K}\Bigg{)}^{-1}\sum_{i=1}^{N}\hat{\mbox{$\mbox{$\mathbf{u}$}$}}_{i}^{\top}y_{ij},\end{split}$
(64)

No derivatives are used to derive ALS (a MAP estimation procedure). It is also
possible to update all quantities in parallel by using all learning rates
$<1$, and convergence is less of an issue unlike VMP [Paquet, 2014].

Finally, we briefly discuss the inclusion of non-conjugate terms.

###### Ex 6

Suppose that, in 1, we decided to use a non-conjugate prior, say, a logit-
normal distribution with mean $m$,
$p(\pi_{0})\propto\frac{1}{\pi_{0}(1-\pi_{0})}\exp\left[-\frac{(\text{logit}(\pi_{0})-m)^{2}}{2}\right].$
Still, the derivation proceeds in the same way by expanding
$\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{Y}$}$}|\mbox{$\mbox{$\mathbf{z}$}$})p(\mbox{$\mbox{$\mathbf{z}$}$})p(\pi_{0})]=$

$\displaystyle\Bigg{(}\sum_{i=1}^{N}\underbrace{\mathbb{E}_{q}(z_{i})}_{=\mu_{i}}\underbrace{\mathbb{E}_{q}{\left[\begin{array}[]{c}\log\pi_{0}\\\
\log(1-\pi_{0})\end{array}\right]}^{\top}}_{=\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{0}^{\top}}\underbrace{\left[\begin{array}[]{c}p_{a}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})\\\
p_{b}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})\end{array}\right]}_{\text{Coeff. in
front}}\Bigg{)}+N\underbrace{\mathbb{E}_{q}[\log(1-\pi_{0})]}_{\mu_{0}^{(2)}}+\underbrace{\mathbb{E}_{q}[\log
p(\pi_{0})]}_{\textrm{non-conjugate term}}.$

and combine the last two terms by collecting the linear terms together,

$\underbrace{\mathbb{E}_{q}{\left[\begin{array}[]{c}\log\pi_{0}\\\
\log(1-\pi_{0})\end{array}\right]}^{\top}}_{=\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{0}^{\top}}\underbrace{\left[\begin{array}[]{c}-1\\\
N-1\end{array}\right]}_{\text{Coeff. in
front}}+\underbrace{\mathbb{E}_{q}\left[-\frac{(\text{logit}(\pi_{0})-m)^{2}}{2}\right]}_{\textrm{non-
conjugate term}}.$

We will now just have to add the last term to the derivative of
$\nabla_{\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{0}}$, and rewrite the
update. The form of the update does not change much, rather the non-conjugate
prior can simply be written as a “pseudo” conjugate Beta prior,

$\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{0}\leftarrow(1-\rho_{0})\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{0}+\rho_{0}\left(\left[\begin{array}[]{c}\hat{\alpha}_{0}-1\\\
N+\hat{\beta}_{0}-1\end{array}\right]+\sum_{i=1}^{N}\mu_{i}\left[\begin{array}[]{c}p_{a}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})\\\
p_{b}(\mbox{$\mbox{$\mathbf{y}$}$}_{i})\end{array}\right]\right).$

where
$\hat{\alpha}_{0}=\nabla_{\mu_{0}^{(1)}}\mathbb{E}_{q}\left[-\mbox{$\frac{1}{2}$}(\text{logit}(\pi_{0})-m)^{2}\right]$
and
$\hat{\beta}_{0}=\nabla_{\mu_{0}^{(2)}}\mathbb{E}_{q}\left[-\mbox{$\frac{1}{2}$}(\text{logit}(\pi_{0})-m)^{2}\right]$
are the parameters of the pseudo beta prior. The gradients automatically
convert the non-conjugate prior form into a conjugate one, which is due to
Khan and Lin [2017].

The elegance of the update above is not a coincidence, rather it is by design.
Khan and Lin [2017] linearize the whole VB objective with respect to all
$\boldsymbol{\mu}$. This linearizes the non-conjugate terms while preserving
the linearity of the conjugate terms; see [Khan and Lin, 2017, Lemma 1]. This
is a remarkable property of such linearizations. The BLR builds on this but it
is the linearization which enables generalization to all sorts of algorithms.
This insight is under-appreciated so far, but we do hope that the examples
shown in this paper motivates many readers to exploit such (Bayesian)
linearization procedure.

##  3 Future Work

One drawback of the recipe is that it depends heavily on the
$(\mbox{$\mbox{$\boldsymbol{\lambda}$}$},\mbox{$\mbox{$\boldsymbol{\mu}$}$})$
parameterization. There are several works that extend beyond this, for
example, to mixture distributions [Lin et al., 2019], structured Gaussian
distributions [Lin et al., 2021], and also to transformation families [Kıral
et al., 2023]. Such approaches have shown to be useful in deriving existing
algorithms, as well as designing new ones [Khan and Rue, 2021]. We believe
that such extensions exist for generic distributions and that it is possible
to derive all sorts of updates in the same fashion as described in this paper.

Acknowledgments

I would like to thank Dr. Negar Safinianaini and other members of our group
for discussions which not only compelled me to write this paper but also
allowed me to see the value of techniques I have been using to derive VB for
many years.

## References

  * Bishop [2006] C. M. Bishop.  _Pattern Recognition and Machine Learning (Information Science and Statistics)_.  Springer-Verlag, Berlin, Heidelberg, 2006.  ISBN 0387310738.
  * Blei et al. [2017] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe.  Variational inference: A review for statisticians.  _Journal of the American statistical Association_ , 112(518):859–877, 2017.
  * Hinton and Van Camp [1993] G. E. Hinton and D. Van Camp.  Keeping the neural networks simple by minimizing the description length of the weights.  In _Annual Conference on Computational Learning Theory_ , pages 5–13, 1993.
  * Hoffman et al. [2013] M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley.  Stochastic variational inference.  _The Journal of Machine Learning Research_ , 14(1):1303–1347, 2013.
  * Jaakkola and Jordan [1996] T. Jaakkola and M. Jordan.  A variational approach to Bayesian logistic regression problems and their extensions.  In _International conference on Artificial Intelligence and Statistics_ , 1996.
  * Khan and Lin [2017] M. E. Khan and W. Lin.  Conjugate-computation variational inference: converting variational inference in non-conjugate models to inferences in conjugate models.  In _International Conference on Artificial Intelligence and Statistics_ , pages 878–887, 2017.
  * Khan and Rue [2021] M. E. Khan and H. Rue.  The Bayesian learning rule.  _arXiv preprint arXiv:2107.04562_ , 2021.
  * Khan et al. [2016] M. E. Khan, R. Babanezhad, W. Lin, M. Schmidt, and M. Sugiyama.  Faster stochastic variational inference using proximal-gradient methods with general divergence functions.  In _Proceedings of the Conference on Uncertainty in Artificial Intelligence_ , 2016.
  * Kıral et al. [2023] E. M. Kıral, T. Möllenhoff, and M. E. Khan.  The Lie-group Bayesian learning rule.  In _International conference on Artificial Intelligence and Statistics_ , 2023.
  * Lin et al. [2019] W. Lin, M. E. Khan, and M. Schmidt.  Fast and simple natural-gradient variational inference with mixture of exponential-family approximations.  2019\.
  * Lin et al. [2021] W. Lin, F. Nielsen, M. E. Khan, and M. Schmidt.  Tractable structured natural-gradient descent using local parameterizations.  In _International Conference on Machine Learning_ , pages 6680–6691. PMLR, 2021.
  * Paquet [2014] U. Paquet.  On the convergence of stochastic variational inference in Bayesian networks.  _NIPS Workshop on variational inference_ , 2014.
  * Saul et al. [1996] L. K. Saul, T. Jaakkola, and M. I. Jordan.  Mean field theory for sigmoid belief networks.  _Journal of Artificial Intelligence Research_ , 4:61–76, 1996.
  * Tipping and Bishop [1999] M. E. Tipping and C. M. Bishop.  Probabilistic principal component analysis.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 61(3):611–622, 1999.
  * Winn and Bishop [2005] J. Winn and C. M. Bishop.  Variational message passing.  _Journal of Machine Learning Research_ , 6(Apr):661–694, 2005.

##  Appendix A VB Stationary Point, The Linearity Property and Equivalence of
Eqs. 1 and 2

VB stationary point: We first give a short proof of the stationarity condition
Eq. 1, which essentially follows by rewriting the ELBO as a function of
$q_{i}$ alone and expressing it as a Kullback-Leibler divergence (KLD) term,

$\begin{split}\mathcal{L}(q)&=\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{y}$}$},\mbox{$\mbox{$\mathbf{z}$}$})]+\sum_{i}\mathbb{E}_{q}[-\log
q_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})]\\\
&=\mathbb{E}_{q_{i}}\left[\mathbb{E}_{q_{/i}}[\log
p(\mbox{$\mbox{$\mathbf{y}$}$},\mbox{$\mbox{$\mathbf{z}$}$})]-\log
q_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})\right]+\text{const.}\\\
&=\mathbb{D}_{\text{KL}}[q_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})\,\|\,\tilde{p}_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})]+\text{const.}\end{split}$
(65)

where
$\tilde{p}_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})\propto\exp[\mathbb{E}_{q_{\backslash
i}}[\log p(\mbox{$\mbox{$\mathbf{y}$}$},\mbox{$\mbox{$\mathbf{z}$}$})]]$ is a
distribution where the rest $\mbox{$\mbox{$\mathbf{z}$}$}_{\backslash i}$ is
marginalized out. The minimum occurs when the two arguments in the divergence
are equal, which gives us the condition in Eq. 1.

The linearity property: This is due to the fact that the CE terms in log-joint
are multi-linear in sufficient statistics [Winn and Bishop, 2005, p. 668].
Below, we give a formal statement which is the basis of Step 1.

###### Theorem 1

If the conditional
$p(\mbox{$\mbox{$\mathbf{z}$}$}_{i}|\mbox{$\mbox{$\mathbf{z}$}$}_{\backslash
i},\mbox{$\mbox{$\mathbf{y}$}$})$ is conjugate to
$q_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})$, that is, if there exists
$\mbox{$\mbox{$\boldsymbol{\eta}$}$}_{i}(\cdot)$ such that the log-conditional
can be expressed in terms of
$\mbox{$\mbox{$\mathbf{T}$}$}_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})$ as
follows,

$\log
p(\mbox{$\mbox{$\mathbf{z}$}$}_{i}|\mbox{$\mbox{$\mathbf{z}$}$}_{\backslash
i},\mbox{$\mbox{$\mathbf{y}$}$})=\langle\mbox{$\mbox{$\mathbf{T}$}$}_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})\,,\,\mbox{$\mbox{$\boldsymbol{\eta}$}$}_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{\backslash
i},\mbox{$\mbox{$\mathbf{y}$}$}_{i})\rangle+\text{const.},$ (66)

then $\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{y}$}$},\mbox{$\mbox{$\mathbf{z}$}$})]=\langle\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{i}\,,\,\mathbb{E}_{q_{\backslash
i}}\left[\mbox{$\mbox{$\boldsymbol{\eta}$}$}_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{\backslash
i},\mbox{$\mbox{$\mathbf{y}$}$}_{i})\right]\rangle+\text{const.}$, which is
linear with respect to $\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{i}$.

This also justifies Step 2 because the gradient with respect to
$\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{i}$ becomes

$\nabla_{\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{i}}\mathbb{E}_{q}[\log
p(\mbox{$\mbox{$\mathbf{y}$}$},\mbox{$\mbox{$\mathbf{z}$}$})]=\mathbb{E}_{q_{\backslash
i}}\left[\mbox{$\mbox{$\boldsymbol{\eta}$}$}_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{\backslash
i},\mbox{$\mbox{$\mathbf{y}$}$}_{i})\right]$ (67)

which is simply the term in front of $\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{i}$.

Equivalence: Due to Eq. 2, we can show that, at a fixed point, the optimal
natural parameter is equal to the coefficient in front of
$\mbox{$\mbox{$\boldsymbol{\mu}$}$}_{i}$,

$\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{i}^{*}=\nabla_{\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{i}}\left.\mathbb{E}_{q}\left[\log
p(\mbox{$\mbox{$\mathbf{y}$}$},\mbox{$\mbox{$\mathbf{z}$}$})\right]\right|_{\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{i}=\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{i}^{*}}=\mathbb{E}_{q_{\backslash
i}^{*}}\left[\mbox{$\mbox{$\boldsymbol{\eta}$}$}_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{\backslash
i},\mbox{$\mbox{$\mathbf{y}$}$}_{i})\right],$ (68)

Using this, we can show that Eq. 2 leads to Eq. 1 by using the definition of
$q_{i}^{*}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})$,

$\begin{split}\log
q_{i}^{*}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})&=\langle\mbox{$\mbox{$\mathbf{T}$}$}(\mbox{$\mbox{$\mathbf{z}$}$}_{i}),\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{i}^{*}\rangle+c,\\\
&=\langle\mbox{$\mbox{$\mathbf{T}$}$}(\mbox{$\mbox{$\mathbf{z}$}$}_{i}),\mathbb{E}_{q_{\backslash
i}^{*}}\left[\mbox{$\mbox{$\boldsymbol{\eta}$}$}_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{/i},\mbox{$\mbox{$\mathbf{y}$}$}_{i})\right]\rangle+c\\\
&=\mathbb{E}_{q_{\backslash i}^{*}}\left[\log
p(\mbox{$\mbox{$\mathbf{y}$}$},\mbox{$\mbox{$\mathbf{z}$}$})\right]+c^{\prime},\end{split}$
(69)

where the second and third line follow by using Eq. 68 and Eq. 66 respectively
($c$ and $c^{\prime}$ are constants that do not depend on
$\mbox{$\mbox{$\mathbf{z}$}$}_{i}$). The derivation is reversible and we can
also derive Eq. 2 from Eq. 1, making the two statements equivalent whenever
the posterior is a conjugate-exponential family.

Extension for non-constant base-measures: If the exponential-family includes a
base-measure $h_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})$ as shown below,

$q_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})\propto
h_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})\exp\left(\langle\text{\mbox{$\mbox{$\mathbf{T}$}$}}_{i}(\text{\mbox{$\mbox{$\mathbf{z}$}$}}_{i}),\text{\mbox{$\mbox{$\boldsymbol{\lambda}$}$}}_{i}\rangle\right),$

then the condition in Eq. 2 is modified to the following,

$\mbox{$\mbox{$\boldsymbol{\lambda}$}$}_{i}^{*}=\nabla_{\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{i}}\left.\mathbb{E}_{q}\left[\log
p(\mbox{$\mbox{$\mathbf{y}$}$},\mbox{$\mbox{$\mathbf{z}$}$})-\log
h_{i}(\mbox{$\mbox{$\mathbf{z}$}$}_{i})\right]\right|_{\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{i}=\text{\mbox{$\mbox{$\boldsymbol{\mu}$}$}}_{i}^{*}},$
(70)

The proof is exactly the same so we omit it.