I am trying to setup VB to do a weighted linear regression for vector observations. My setup is that I have $N$ numbers of $d$-dimensional vector observations. I would like to model the noise as being independent for an output observation but I would like each point to have its own variance. In view of that, my model is as follows:

$$ y_i \sim \mathcal{N}(T(x_i; \beta), \frac{\sigma^2}{w_i} {\textbf{I}}) \\ \beta \sim \mathcal{N} (\beta_{0}, \Sigma_{0}) \\ w_i \sim \mathcal{G} (a, b) \\ $$

Here $y_i$ is the one $d$-dimensional output observations and $x_i$ is the corresponding independent variable. They are related by some transformation $T$ which is parameterized by $\beta$. The output variance for a point $i$ is scaled by $w_i$ and $\sigma^2$ is some given global variance.

The full joint model can be written as:

$$ p(\beta, w, y |x) = \prod_{i=1}^{N} \bigg[p(y_i|x_i, w_i, H) \ p(w_i)\bigg] \ p(\beta) $$

Applying the mean field approximation $Q(w, \beta) = Q(w)Q(\beta)$ and getting the optimal distribution for $w_i$'s we need to take the expectation with respect to the $\beta$ parameters. We can work with the logarithm to make things a bit simpler

$$ \ln q(w^*) = \bigg\langle \ln \prod_{i=1}^{N} \bigg[p(y_i|x_i, w_i, \beta) \ p(w_i)\bigg] \ p(\beta)\bigg\rangle_{\beta} \\ =\bigg\langle \ln \prod_{i=1}^{N} \bigg[p(y_i|x_i, w_i, H) \ p(w_i)\bigg]\bigg\rangle_{\beta} + C $$

Now, we look at the term $p(y_i|x_i, w_i, H) \ p(w_i)$ in more detail. Expanding we have:

$$ \prod_{i=1}^{N} \frac{1}{\sqrt{2 \pi}} \big(\frac{\sigma^2}{w_i} {\textbf{I}}\big)^{-\frac{1}{2}} \exp \bigg(-0.5 \times [y_i - T(x_i)]^T \big(\frac{\sigma^2}{w_i} {\textbf{I}}\big)^{-1} \ [y_i - T(x_i)]\bigg) \ \bigg(\frac{b^a}{\Gamma(a)} w_i^{a-1} \exp{(-b w_i)}\bigg) $$

Dropping the constants we get:

$$ \propto \prod_{i=1}^{N} \big(\frac{\sigma^2}{w_i} {\textbf{I}}\big)^{-\frac{1}{2}} \exp \bigg(-0.5 \times [y_i - T(x_i)]^T \big(\frac{\sigma^2}{w_i} {\textbf{I}}\big)^{-1} \ [y_i - T(x_i)]\bigg) \bigg(w_i^{a-1} \exp{(-b w_i)}\bigg) \\ $$

Gathering some like terms we get:

$$ \propto \prod_{i=1}^{N} \big(\frac{\sigma^2}{w_i} {\textbf{I}}\big)^{-\frac{1}{2}} \ w_i^{a-1} \exp \bigg(-0.5 \times [y_i - T(x_i)]^T \big(\frac{\sigma^2}{w_i} {\textbf{I}}\big)^{-1} \ [y_i - T(x_i)] - b w_i\bigg) $$

Not sure if I should simplify this further and try and absorb the scalar $w_i$s into the matrix terms but it still does not help me somehow.

Now I have two issues. I was hoping at this point, this will start to look like some familiar distribution. It sort of looks like some Gamma distribution perhaps but I am not sure how I can simplify it to get it in that form.

My second issue is that after I have done that, I need to still take the expectation with respect to $Q(\beta)$ and I am not sure how to do that.

Any pointers on how to proceed next would be very appreciated!

[EDIT] So, I have been thinking about this. I am assuming I can write the joint model for each $i$ term as:

$$ p(y_i, w_i, \beta|x_i) = \prod_{i=1}^{N}p(y_i|x_i, w_i, \beta) p(w_i) \frac{p(\beta)}{N} $$

Now, for example, when computing the distribution for $w_i$ with the mean field approximation, I need to take the expectation wrt to $\beta$. So, expanding this, I have:

$$ \propto \sqrt{w_i} \exp{\bigg(-[y_i - T(x_i, \beta)]^T w_i [y_i - T(x_i, \beta)]\bigg)} \\ \bigg( w_i^{a-1} \exp{(-b w_i)}\bigg) \exp{\bigg(-(\beta - \beta_0)^T \Sigma_{ 0}(\beta - \beta_0) \bigg)} $$

Now, the first two exponentials come together nicely to make a Gamma distribution and I am not sure though how to manipulate the last exponential so that the posterior over $w_i$ gets some nice form.


1 Answer 1


Variational Bayes is a tricky beast, because it is extremely often very poorly explained. Hopefully I can do a satisfactory job.

The variational bayes method start from the variational lower bound. It's true that, for any probability distribution $q$

$$ \forall q, \log ( \int p ) \geq \int q \log ( \frac{p}{q} ) $$

This is obtained quite simply by applying Jensen's inequality on the concave function $\log$. The optimal q is of course the normalized version of $p$:

$$ q^\star = \frac{p}{\int p} $$

A basic idea that is shared by a lot of slightly different methods (which are all refered to as "variational bayes" in a complete abuse of language) is then to maximize the quantity on the right for $q \in \mathcal Q$ restricted to a sub-ensemble of all probability distributions $\mathcal Q$. Let me list some popular choices:

  • what you are doing here: $\mathcal Q$ is the space of all probability distributions which factorize in a certain way. I don't know what to call this method though I think that it is most often called VB-EM (?). This only works in some very specialized cases: only when the expected values which you struggle to compute give you simple probability distributions

  • using a simple approximating family: for example $\mathcal Q$ is all Gaussian distributions, or all gaussian distributions which factorize in a certain way. This variational method is much more general than the preceding one, but you need to use gradient ascent to find its local maxima. Since those gradients can be hard to compute, it might be useful to use stochastic gradient methods.

If you are in a case where VB-EM can't be applied (which is most cases: it isn't a very general method; but I haven't looked at your math very carefully), but you want to use the variational lower-bound, you should look at the second class of VB methods.

You could also look at other approximate inference methods (Newton methods and expectation propagation)

  • $\begingroup$ Thanks. This bit I understand. I think my issue is that as my model stands, I might not have used the right conjugate priors. As far as my calculations go, I do not see the posterior $q(w)$ come out as a good/managable distribution that I can write down and that is what I am trying to figure out. $\endgroup$
    – Luca
    Oct 22, 2015 at 14:03
  • $\begingroup$ I looked a little bit at your math. Your prior for w_i is a Gamma right ? That's the conjugate prior so it should work. Your hard step is the relationship between beta and the (conditional) mean of y_i. Is T simple enough so that you can compute the mean value under simple distributions over beta ? $\endgroup$ Oct 26, 2015 at 10:31
  • $\begingroup$ It is a gamma and that is what I thought. For this purpose, we can assume $T$ is simple enough but I could not even get it to show the posterior over $w_i$ as a gamma distribution, so far... $\endgroup$
    – Luca
    Oct 26, 2015 at 10:33
  • $\begingroup$ $\log p(y_i | T_i, w_i ) = - 0.5 w_i * \alpha(y_i, T_i) - 0.5 \log (w_i)$. Then it's just a matter of computing the average of $\alpha(y_i, T_i)$ under the probability distribution over $\beta$ $\endgroup$ Oct 26, 2015 at 12:18
  • $\begingroup$ hmmmmm... I will need to think about this. Any chance you could elaborate a bit? I am quite new to this. I do not understand what the $\alpha$ parameter is. $\endgroup$
    – Luca
    Oct 26, 2015 at 12:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.