Estimate the Kullback Leibler (KL) divergence with monte carlo

I want to estimate the KL divergence between two continuous distributions f and g. However, I can't write down the density for either f or g. I can sample from both f and g via some method (for example, markov chain monte carlo).

The KL divergence from f to g is defined like this

$$D_{KL}(f || g) = \int_{-\infty}^{\infty} f(x) \log\left(\frac{f(x)}{g(x)}\right) dx$$

This is the expectation of $\log\left(\frac{f(x)}{g(x)}\right)$ with respect to f so you could imagine some monte carlo estimate

$$\frac{1}{N}\sum_i^N \log\left(\frac{f(x_i)}{g(x_i)}\right)$$

Where i indexes N samples that are drawn from f (i.e. $x_i \sim f()$ for i = 1, ..., N)

However, since I don't know f() and g(), I can't even use this monte carlo estimate. What is the standard way of estimating the KL in this situation?

EDIT: I do NOT know the unnormalized density for either f() or g()

• Have you considered using the ecdfs? – Toby May 21 '17 at 20:53
• this will work but it can be arbitrarily slow for hard choice of f and g (close, or close tails). If you decide to ignore samples away from tails then you might have more luck with upper bounding the roc. – enthdegree May 24 '17 at 5:39
• Essentially a duplicate: stats.stackexchange.com/questions/211175/… – kjetil b halvorsen May 24 '17 at 8:21

Here I assume that you can only sample from the models; an unnormalized density function is not available.

You write that

$$D_{KL}(f || g) = \int_{-\infty}^{\infty} f(x) \log\left(\underbrace{\frac{f(x)}{g(x)}}_{=: r}\right) dx,$$

where I have defined the ratio of probabilities to be $r$. Alex Smola writes, although in a different context that you can estimate these ratios "easily" by just training a classifier. Let us assume you have obtained a classifier $p(f|x)$, which can tell you the probability that an observation $x$ has been generated by $f$. Note that $p(g|x) = 1 - p(f|x)$. Then:

$$r = \frac{p(x|f)}{p(x|g)} \\ = \frac{p(f|x) {p(x) p(g)}}{p(g|x)p(x) p(f)} \\ = \frac{p(f|x)}{p(g|x)},$$

where the first step is due to Bayes and the last follows from the assumption that $p(g) = p(f)$.

Getting such a classifier can be quite easy for two reasons.

First, you can do stochastic updates. That means that if you are using a gradient-based optimizer, as is typical for logistic regression or neural networks, you can just draw a samples from each $f$ and $g$ and make an update.

Second, as you have virtually unlimited data–you can just sample $f$ and $g$ to death–you don't have to worry about overfitting or the like.

I assume you can evaluate $f$ and $g$ up to a normalizing constant. Denote $f(x) = f_u(x)/c_f$ and $g(x) = g_u(x)/c_g$.

A consistent estimator that may be used is $$\widehat{D_{KL}}(f || g) = \left[n^{-1} \sum_j f_u(x_j)/\pi_f(x_j)\right]^{-1}\frac{1}{N}\sum_i^N \left[\log\left(\frac{f_u(z_i)}{g_u(z_i)}\right)\frac{f_u(z_i)}{\pi_r(z_i)}\right] - \log (\hat{r})$$ where $$\hat{r} = \frac{1/n}{1/n}\frac{\sum_j f_u(x_j)/\pi_f(x_j)}{\sum_j g_u(y_j)/\pi_g(y_j)} \tag{1}.$$ is an importance sampling estimator for the ratio $c_f/c_g$. Here you use $\pi_f$ and $\pi_g$ as instrumental densities for $f_u$ and $g_u$ respectively, and $\pi_r$ to target the log ratio of unnormalized densities.

So let $\{x_i\} \sim \pi_f$, $\{y_i\} \sim \pi_g$, and $\{z_i\} \sim \pi_r$. The numerator of (1) converges to $c_f$. The denominator converges to $c_g$. The ratio is consistent by the continuous mapping theorem. The log of the ratio is consistent by continuous mapping again.

Regarding the other part of the estimator, $$\frac{1}{N}\sum_i^N \left[\log\left(\frac{f_u(z_i)}{g_u(z_i)}\right)\frac{f_u(z_i)}{\pi_r(z_i)}\right] \overset{\text{as}}{\to} c_f E\left[ \log\left(\frac{f_u(z_i)}{g_u(z_i)}\right) \right]$$ by the law of large numbers.

My motivation is the following:

\begin{align*} D_{KL}(f || g) &= \int_{-\infty}^{\infty} f(x) \log\left(\frac{f(x)}{g(x)}\right) dx \\ &= \int_{-\infty}^{\infty} f(x)\left\{ \log \left[\frac{f_u(x)}{g_u(x)} \right] + \log \left[\frac{c_g}{c_f} \right]\right\} dx \\ &= E_f\left[\log \frac{f_u(x)}{g_u(x)} \right] + \log \left[\frac{c_g}{c_f} \right] \\ &= c_f^{-1} E_{\pi_r}\left[\log \frac{f_u(x)}{g_u(x)}\frac{f_u(x)}{\pi_r(x)} \right] + \log \left[\frac{c_g}{c_f} \right]. \end{align*} So I just break it up into tractable pieces.

For more ideas on how to simulate the likelhood ratio, I found a paper that has a few: https://projecteuclid.org/download/pdf_1/euclid.aos/1031594732

• (+1) It's worth noting here that importance sampling can have extremely high variance (even infinite variance) if the target distribution has fatter tails than the distribution you're sampling from and/or the number of dimensions is at all large. – David J. Harris May 24 '17 at 18:22
• @DavidJ.Harris very very true – Taylor May 24 '17 at 19:46