Skip to content

Instantly share code, notes, and snippets.

@atabakd
Last active July 5, 2024 19:47
Show Gist options
  • Save atabakd/ed0f7581f8510c8587bc2f41a094b518 to your computer and use it in GitHub Desktop.
Save atabakd/ed0f7581f8510c8587bc2f41a094b518 to your computer and use it in GitHub Desktop.
KL divergence for multivariate samples
# 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.))
@mtcesarini
Copy link

@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

@TessaHayman-Aimsun
Copy link

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

@GalRaz
Copy link

GalRaz commented Jun 7, 2022

@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

@VolodyaCO
Copy link

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.

@dm2232000
Copy link

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment