# Kullback-Leibler Divergence for two samples

I tried to implement a numerical estimate of the Kullback-Leibler Divergence for two samples. To debug the implementation draw the samples from two normal distributions $\mathcal N (0,1)$ and $\mathcal N (1,2)$.

For a simple estimate I generated two histograms and tried to numerically approximate the integral. I got stuck with handling those parts of the histogram where the bins of one of the histogram are zero such that I either end up with dividing by zero or the logarithm of zero. How do I handle this issue?

A related question came in my mind: How to exactly compute the KL-Divergence between two different uniform distributions? Do I have to restrict the integral to the union of the support of both distributions?

• Well, the support of the normal distribution is the set of real numbers. There's no problem in pure math, but yeah, for your numerical approximation, you need to make sure your sample size is large enough relative to the region you want to integrate over. You won't be able to integrate over (-inf, +inf) like you can in pure math... Go for something reasonable? If you're more than 3 standard deviations away from the mean, it's going to be pretty thin... – Matthew Gunn May 6 '16 at 8:29
• With respect to your second question, the KL-divergence between two different uniform distributions is undefined ($\log(0)$ is undefined). Similarly, the KL-divergence for two empirical distributions is undefined unless each sample has at least one observation with the same value as every observation in the other sample. – jbowman Nov 29 '16 at 17:28
• @jbowman Small note. Although you are right that $\log(0)$ is undefined (or $-\infty$), it is customary in information theory to treat $\log(0) \cdot 0$ as $0$. – Luca Citi Apr 12 '17 at 20:45
• A similar question: mathoverflow.net/questions/119752/… – kjetil b halvorsen Jul 27 '17 at 14:43

The Kullback-Leibler divergence is defined as $$\DeclareMathOperator{\KL}{KL} \KL(P || Q) = \int_{-\infty}^\infty p(x) \log \frac{p(x)}{q(x)} \; dx$$ so to calculate (estimate) this from empirical data we would need, maybe, some estimates of the density functions $p(x), q(x)$. So a natural starting point could be via density estimation (and after that, just numerical integration). How good or stable such a method would be I don't know.

But first your second question, then I will return to the first one. Lets say $p$ and $q$ are uniform densities on $[0,1]$ and $[0,10]$ respectively. Then $\KL(p || q)=\log 10$ while $\KL(q || p)$ is more difficult to define, but the only reasonable value to give it is $\infty$, as far as I can see, since it involves integrating $\log(1/0)$ which we can choose to interprete as $\log \infty$. This results are reasonable from the interpretation I give in Intuition on the Kullback-Leibler (KL) Divergence

Returning to the main question. It is asked in a very nonparametric way, and no assumptions are stated on the densities. Probably some assumptions are needed. But assuming the two densities are proposed as competing models for the same phenomenon, we can probably assume they have the same dominating measure: KL divergence between a continuous and a discrete probability distribution would always be infinity, for example. One paper addressing this question is the following: https://pdfs.semanticscholar.org/1fbd/31b690e078ce938f73f14462fceadc2748bf.pdf They propose a method which do not need preliminary density estimation, and analyses its properties.

(There are many other papers). I will come back and post some details from that paper, the ideas.

 EDIT


Some ideas from that paper, which is about estimation of KL divergence with iid samples from absolutely continuous distributions. I show their proposal for one-dimensional distributions, but they give a solution for vectors also (using nearest neighbor density estimation). For proofs read the paper!

They propose to use a version of the empirical distribution function, but interpolated linearly between sample points to get a continuous version. They define $$P_e(x) = \frac1{n}\sum_{i=1}^n U(x-x_i)$$ where $U$ is the Heavyside step function, but defined so that $U(0)=0.5$. Then that function interpolated linearly (and extended horizontally beyond the range) is $P_c$ ($c$ for continuous). Then they propose to estimate the Kullback-Leibler divergence by $$\hat{D}(P \| Q) = \frac1{n}\sum_{i=1}^n \log\left(\frac{\delta P_c(x_i)}{\delta Q_c(x_i)}\right)$$ where $\delta P_c = P_c(x_i)-P_c(x_i-\epsilon)$ and $\epsilon$ is a number smaller than the smallest spacing of the samples.

R code for the version of the empirical distribution function that we need is

my.ecdf  <-  function(x)   {
x   <-   sort(x)
x.u <-   unique(x)
n  <-  length(x)
x.rle  <-  rle(x)$lengths y <- (cumsum(x.rle)-0.5) / n FUN <- approxfun(x.u, y, method="linear", yleft=0, yright=1, rule=2) FUN }  note that rle is used to take care of the case with duplicates in x. Then the estimation of the KL divergence is given by KL_est <- function(x, y) { dx <- diff(sort(unique(x))) dy <- diff(sort(unique(y))) ex <- min(dx) ; ey <- min(dy) e <- min(ex, ey)/2 n <- length(x) P <- my.ecdf(x) ; Q <- my.ecdf(y) KL <- sum( log( (P(x)-P(x-e))/(Q(x)-Q(x-e)))) / n KL }  Then I show a small simulation: KL <- replicate(1000, {x <- rnorm(100) y <- rt(100, df=5) KL_est(x, y)}) hist(KL, prob=TRUE)  which gives the following histogram, showing (an estimation) of the sampling distribution of this estimator: For comparison, we calculate the KL divergence in this example by numerical integration: LR <- function(x) dnorm(x,log=TRUE)-dt(x,5,log=TRUE) 100*integrate(function(x) dnorm(x)*LR(x),lower=-Inf,upper=Inf)$value
 3.337668


hmm ... the difference being large enough that there is much here to investigate!

Expanding a little bit on kjetil-b-halvorsen's answer, and sorry for not commenting, I don't have the reputation:

1. I have the feeling that the analytical computation should be (without multiplication by 100):

LR <- function(x) dnorm(x,log=TRUE)-dt(x,5,log=TRUE) integrate(function(x) dnorm(x)*LR(x),lower=-Inf,upper=Inf)$value  1. If I'm right, the estimator$\hat D(P||Q)$does not converge to the KL divergence, but the convergence is stated as:$\hat D(P||Q)-1 \to D(P||Q)\$. The arrow represents a.s convergence.

Once those two corrections are made, the results seem more realistic.

• Thanks, I will look into this and update my answer. – kjetil b halvorsen Apr 16 '18 at 11:47