Help with Variational Bayes on a weighted linear regression model

Question

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.

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)