KL divergence for a hierarchical prior structure e.g. Linear Regression

Question

For a Linear Regression $\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \epsilon$ with $\epsilon \sim \mathcal{N}(0, \sigma^2\mathbb{I})$, suppose the prior set on $\beta_k$ is $\sim \mathcal{N}(0, l_k)$ and prior on $l_k$ is $\mathbb{G}(\alpha_k, \beta_k)$ which is a Gamma Distribution.

Suppose the approximate posterior on $\beta_k$ is $q(\beta_k)$ and that of $l_k$ is $q(l_k)$ then if the KL divergence from $q$ to $p$ i.e. $KL(q||p)$ is to be computed then what is the recommended steps? This is required in Variational Bayesian Inference.

I have attempted like this...

$KL(q||p) = \displaystyle \prod_{k = 1}^P \int_{\beta_k, l_k} q(\beta_k) q(l_k) \log \frac{p(\beta_k|l_k) p(l_k)}{q(\beta_k) q(l_k)} d\beta_k dl_k $ $ \displaystyle \int_{\beta_k, l_k} q(\beta_k) q(l_k) \log \frac{p(\beta_k|l_k) p(l_k)}{q(\beta_k) q(l_k)} d\beta_k dl_k\\ = \displaystyle \int_{ l_k} \left[ q(l_k)\left\{\int_{\beta_k} q(\beta_k)\log q(\beta_k|l_k)p(l_k) d\beta_k\right\} - \log q(l_k)\right] dl_k - \int_{\beta_k} q(\beta_k) \log q(\beta_k) d\beta_k$

After this how to proceed? Can anybody help me? Thanks in advance.

Answer 1

Getting a closed-form solution to this problem may be quite difficult, but a Monte Carlo approach can allow you to solve a much simpler problem and simulate in order to estimate the impact of variation in l_k with regard to the KL divergence.

Since your residuals are normally-distributed and your parameter priors are likewise normally-distributed, congratulations! You're in conjugate Gaussian prior territory which leads to very straightforward estimation formulation (and corresponding KL-divergence calcs). The estimation itself from the posterior basically equates to penalized least squares (when the model is linear) with an L2-penalty on deviation from the prior.

Start by fixing your parameter prior distribution with respect to l_k (pretend that l_k is precisely known at the outset using the mean of the gamma distribution). Taking the log-likelihood of the posterior distribution leads to a very friendly estimation form. You can use the Fisher information (from the second derivative of the log posterior probability) to estimate posterior parameter uncertainty and compare to your prior distribution - this question better describes the link between the Fisher matrix and the KL-divergence:

Connection between Fisher metric and the relative entropy

If you then want to see how variation in l_k impacts your solution, Monte Carlo simulation from the gamma distribution would be a good approach (allow l_k to vary in accordance with the gamma distribution, fix it for each of N trials, then calculate the KL divergence for each trial, accumulating the results to yield the estimated distribution). This will allow you to keep the friendly formulation and calculation of KL-div but to generate through simulation the distribution of KL divergence subject to the variational assumptions underlying your priors.

Because your priors are Normal and your data is Normal, the posterior is also normal. When the prior covariance varies in accordance with l_k, you will get a family of posterior distributions that are likewise normal.

Update: it's worth mentioning that MCMC packages such as Stan, PyMC and Edward are very well-suited to estimating posterior distributions on complex problems with non-friendly, non-analytical, non-conjugate probability forms. These won't get you a closed form solution either (for really difficult problems these can be difficult or intractable to attain), but you can often get to a sampled posterior distribution that can be approximated with a multivariate Gaussian for which closed-form KL-divergence is straightforward to estimate.