KL divergence for multivariate samples

kl.py

# https://mail.python.org/pipermail/scipy-user/2011-May/029521.html

import numpy as np

def KLdivergence(x, y):
  """Compute the Kullback-Leibler divergence between two multivariate samples.

  Parameters
  ----------
  x : 2D array (n,d)
    Samples from distribution P, which typically represents the true
    distribution.
  y : 2D array (m,d)
    Samples from distribution Q, which typically represents the approximate
    distribution.

  Returns
  -------
  out : float
    The estimated Kullback-Leibler divergence D(P||Q).

  References
  ----------
  Pérez-Cruz, F. Kullback-Leibler divergence estimation of
continuous distributions IEEE International Symposium on Information
Theory, 2008.
  """
  from scipy.spatial import cKDTree as KDTree

  # Check the dimensions are consistent
  x = np.atleast_2d(x)
  y = np.atleast_2d(y)

  n,d = x.shape
  m,dy = y.shape

  assert(d == dy)


  # Build a KD tree representation of the samples and find the nearest neighbour
  # of each point in x.
  xtree = KDTree(x)
  ytree = KDTree(y)

  # Get the first two nearest neighbours for x, since the closest one is the
  # sample itself.
  r = xtree.query(x, k=2, eps=.01, p=2)[0][:,1]
  s = ytree.query(x, k=1, eps=.01, p=2)[0]

  # There is a mistake in the paper. In Eq. 14, the right side misses a negative sign
  # on the first term of the right hand side.
  return -np.log(r/s).sum() * d / n + np.log(m / (n - 1.))

@patrickshas @atabakd, I and my collegues have found where the sign issue is generated.
The eq. 14 in the reference paper (DOI: 10.1109/ISIT.2008.4595271) has a missing minus sign because the fraction p_k/q_k (cfr. eq. 15,16) has s(x) and r(x) respectively in the denominators and hence, in order to have s(x)/r(x), they should have been powered to -1, which is traduced in a minus sign before the log.
So the formula in line 52 is correct.
Ph.D. M. Cesarini

@mtcesarini I'm still getting negative values for my dataset Have you been able to use it successfully

@mtcesarini same here, using this function some KL divergences are negative for very simple cases like:

np.random.seed(23) 

dist1 = np.random.multivariate_normal(np.array([1, 1]), np.identity(2), 10000)
dist2 = np.random.multivariate_normal(np.array([1, 1]), np.identity(2), 10000)

KLdivergence(dist1,dist2)

returns -0.010910543767634614

Remember that this is an estimation of the KL divergence @GalRaz @palinode-maker. The estimation can be negative. I don't know what the estimation error is, though.

There is indeed an error in the paper and the result found by this code is correct.

In the paper by Perez-Crus, the definition for P(x) is given in terms of c / r(x) in equation 15. Q(x) is given in terms of c / s(x) in equation 16. c is a constant and is the same for both terms. Therefore the ratio P(x) / Q(x) from the definition of Kullback-Leibler Divergence should become c/r(x) * s(x)/c = s(x) / r(x) in equation 14. However, the author made a typo and in the paper this term is r(x) / s(x) in equation 14 which is incorrect.

Please see this paper by Wang et al. (2006) which provides the correct form in their equation 11.

The equation in line 52 of the code is correct, the negative sign corrects for the fact that r(x) and s(x) are flipped. Alternatively, the negative sign could be removed and r(x) and s(x) could be flipped which would be in the same form as equation 11 from the Wang et al. paper referenced above. This is shown below:

return np.log(s/r).sum() * d / n + np.log(m / (n - 1.))

@GalRaz and @TessaHayman-Aimsun Do keep in mind that this is an estimation but that it will converge near the correct value. When comparing two identical distributions such as in your example @GalRaz , the estimation is indeed negative but it is close enough to 0 to be a good enough estimate. Hope this clears up any confusions.