Applying variational inference to this model I am basically trying to do a weighted linear regression in a bayesian way. This is to ensure that the I can take care of the heretoscedastic noise.
So, my model is like:
$$
y_i \sim \mathcal{N}(\beta^Tx_i, \sigma^2/w_i)
$$
Here, $y_i$ are each of the output for the input $x_i$. Now, I have the following distributional assumptions on $\beta$ and $w_i$.
$$
\beta \sim \mathcal{N}(\beta_0, \Sigma_0)
$$
$$
w_i \sim \textrm{Gamma}(a, b)
$$  
Now, I want to approximate the posterior $p(\beta, w|Y, X)$ using VB. So the full log joint model can be written as:
$$
\ln p(Y, \beta, w|X) = \sum_{i=1}^N \ln p(y_i|w_i, \beta, x_i) + \sum_{i=1}^N \ln p(w_i) + \ln p(\beta)
$$
Now, looking at variational inference, I need to minimise the KL-divergence between some appropriately chosen distribution $q(\beta, w)$ and the joint model i.e.
$$
E_q\bigg[\ln \frac{p(Y, \beta, w|X)}{q(\beta, w)}\bigg]
$$
We can expand this as:
$$
E_q\bigg[\sum_{i=1}^N \ln P(y_i|w_i, \beta, x_i)\bigg] + E_q\bigg[\sum_{i=1}^N \ln P(w_i)\bigg] + E_q\bigg[\ln p(\beta)\bigg] - E_q\bigg[\ln q(\beta, w)\bigg]
$$
I can now apply the mean field approximation i.e. $q(\beta, w) \approx q(\beta) q(w)$.
$$
E_{q_{\beta,w}}\bigg[\sum_{i=1}^N \ln P(y_i|w_i, \beta, x_i)\bigg] + E_{q_w}\bigg[\sum_{i=1}^N \ln P(w_i)\bigg] + E_{q_{\beta}}\bigg[\ln p(\beta)\bigg] - E_{q_{\beta}}\bigg[\ln q(\beta)\bigg] - E_{q_{w}}\bigg[\ln q(w)\bigg]
$$
I think the reasoning so far is correct but now I am quite lost as to how to proceed. I am not looking for a full solution here but an idea into how to proceed and what tools I would need to generate the update equations for estimating the parameters of the $\beta$ and $w$ variables.
 A: Check this paper Variational Bayesian inference for linear and logistic regression from Jan Drugowitsch.
A: You derived the variational lower bound in your equation. 
Have a look at the wiki page for variational approximation. https://en.wikipedia.org/wiki/Variational_Bayesian_methods
A: tried to calculate it without any warranty.
$P(y|x\beta,\frac{\sigma}{w_i})P(\beta|\beta_0,\Sigma_0)$
at one point I have
$-\beta(zy\sum{x_i }+ \Sigma_0\beta_0)+0.5(z\sum{x_i}+\Sigma_0)\beta^2$ where I made $z=\sum{\frac{\sigma}{w_i}}$
can you figure out the update equations ?
I guess its like(hope R is fine):
$var = diag(N)*w X'X + \Sigma_0 \\
\beta = solve(var)*(diag(N)*w XY + \Sigma_0 *\beta_0 )$
for the gamma normal update:
$P(Y|x\beta)P(w_i|a,b)$ separated for $ln w_i$ and $w_i$ i got stuff like
$((a-1)-0.5\frac{\sigma}{w^2_i}) ln w_i + (\frac{\sigma}{w^2_i}( y^2-yx\beta + (x\beta)^2)) \frac{w_i}{2} $
ok this might very well be wrong but here the update equations I got.
$a_i =a_0 +0.5*\frac{\sigma}{w^2_i}$ 
$b_i =b_0 +0.5*(\frac{\sigma}{w^2_i} ((y_i- x_i\beta)^2 + var(\beta) ))$
let me know where I made mistakes. and if you understand this.
this $ln\frac{\sigma}{w_i}$ term confuses me. on second thought there should not be a $w_i$ value in the update equation. on second thought $\frac{\sigma}{w^2_i}$ could be very well be wrong and only be $\sigma$
