在 Python 中计算 Kullback–Leibler 散度的有效方法答案

【问题标题】：Efficient way of computing Kullback–Leibler divergence in Python在 Python 中计算 Kullback–Leibler 散度的有效方法
【发布时间】：2016-03-04 14:04:11
【问题描述】：

我必须计算数千个离散概率向量之间的Kullback-Leibler Divergence (KLD)。目前我正在使用以下代码，但它对我的目的来说太慢了。我想知道是否有更快的方法来计算 KL Divergence？

import numpy as np
import scipy.stats as sc

    #n is the number of data points
    kld = np.zeros(n, n)
        for i in range(0, n):
            for j in range(0, n):
                if(i != j):
                    kld[i, j] = sc.entropy(distributions[i, :], distributions[j, :])

【问题讨论】：

如何使用 scipy 在 python 中生成具有最小 KL 散度的概率分布生成器？

标签： python performance numpy scipy statistics

【解决方案1】：

Scipy 的stats.entropy 在其默认意义上邀请输入为一维数组，为我们提供一个标量，这是在列出的问题中完成的。在内部，此函数还允许broadcasting，我们可以在此处滥用以获得矢量化解决方案。

来自docs-

scipy.stats.entropy(pk, qk=None, base=None)

如果只有概率 pk 给定，熵计算为 S = -sum(pk * log(pk), 轴=0)。

如果 qk 不是 None，则计算 Kullback-Leibler 散度 S = 总和（pk * log（pk / qk），轴=0）。

在我们的例子中，我们针对所有行对每一行执行这些熵计算，执行求和以在这两个嵌套循环的每次迭代中获得一个标量。因此，输出数组的形状为(M,M)，其中M 是输入数组中的行数。

现在，这里的问题是 stats.entropy() 将与 axis=0 相加，因此我们将为其提供两个版本的 distributions，这两个版本都会将 rowth-dimension 带到 axis=0 以减少它其他两个轴交错 - (M,1) 和 (1,M) 使用 broadcasting 为我们提供 (M,M) 形状的输出阵列。

因此，解决我们的案例的矢量化且更有效的方法是 -

from scipy import stats
kld = stats.entropy(distributions.T[:,:,None], distributions.T[:,None,:])

运行时测试和验证 -

In [15]: def entropy_loopy(distrib):
    ...:     n = distrib.shape[0] #n is the number of data points
    ...:     kld = np.zeros((n, n))
    ...:     for i in range(0, n):
    ...:         for j in range(0, n):
    ...:             if(i != j):
    ...:                 kld[i, j] = stats.entropy(distrib[i, :], distrib[j, :])
    ...:     return kld
    ...: 

In [16]: distrib = np.random.randint(0,9,(100,100)) # Setup input

In [17]: out = stats.entropy(distrib.T[:,:,None], distrib.T[:,None,:])

In [18]: np.allclose(entropy_loopy(distrib),out) # Verify
Out[18]: True

In [19]: %timeit entropy_loopy(distrib)
1 loops, best of 3: 800 ms per loop

In [20]: %timeit stats.entropy(distrib.T[:,:,None], distrib.T[:,None,:])
10 loops, best of 3: 104 ms per loop

【讨论】：

美丽。谢谢。
@Amir 很高兴为您提供帮助。甚至我都不知道 scipy 的熵会支持广播，它确实做到了，这是 NumPy 拥有的最好的东西！
这正是我的想法，这让我什至不去尝试！您能否给我一个参考，说明您为什么编写像“distrib.T[:,:,None]”和“distrib.T[:,None,:]”这样的矩阵变量？ None 有什么作用，为什么要转置这两个矩阵？
@Amir That None 或np.newaxis 是引入单一维度或将维度从现有的 2D 数组形状扩展到 3D 数组的方法。需要该转置，因为在内部该函数在轴 = 0 上求和，这是为了复制我们在每次迭代的问题的循环代码中获得的所有元素相乘元素的总和。希望这是有道理的！
@KillianTattan 我猜是stats.entropy(distributions[0,None].T, distributions[1:].T)。