KSG互信息估计器的原理详细推导(Kraskov, 2004)

Remote Sensing

于 2023-09-26 12:12:50 发布

阅读量574

点赞数 1

分类专栏：机器学习文章标签：概率论机器学习人工智能线性代数互信息信息论深度学习

本文链接：https://blog.csdn.net/rsstudent/article/details/133303208

版权

机器学习专栏收录该内容

30 篇文章 7 订阅

订阅专栏

KSG Estimators

KSG估计器是Kraskov在2004年提出的互信息估计器，其原文过于简略，我参考大量文献对细节进行了补充。原文位置：[https://arxiv.org/pdf/cond-mat/0305641.pdf]

KSG估计的基本方法是首先对互信息进行如下分解 (1)
$I (X, Y) = H (X) + H (Y) - H (X, Y)$
因此，只需要对 $H (x)$ 等三项进行逐个估计即可。

KLE Esitmator

KLE估计器用于估计微分熵。

利用蒙特卡洛估计有(2)
$H[x]=-\int\mu(x)log\mu(x)dx \approx -\frac{1}{N}log\hat{\mu}(x)$
Then the key problem is estimating $\hat{\mu}$ (x).

考虑一个 $\epsilon-ball$ 以 $x_i$ 为中心, 则 $(0,\frac{\epsilon}{2})$ 的范围内有 $k - 1$ 个点， $(\frac{\epsilon}{2},\frac{\epsilon}{2}+d\epsilon)$ 中间有一个点（the kth neighbor）， $(\frac{\epsilon}{2}+d\epsilon)$ 之外有剩余的 $N - k - 1$ 个点，因此首先在除了 $x_i$ 的N-1个点中挑选k个点，其中一个点位于边缘上，这k个点中剩余选取k-1个点位于最中心范围内,则 $p(\epsilon)$ 是 $x_i$ 周围 $\epsilon-ball$ 的概率质量(mass) (3)：
$\begin{aligned} p_i(\epsilon)=&\tbinom{N-1}{k}\tbinom{k}{k-1}(1-p_i(\epsilon))^{N-k-1}p_i(\epsilon)^{k-1}\frac{dp_i(\epsilon)^1}{d\epsilon}\\ & = \frac{(N-1)!}{(N-k-1)!1!(k-1)!}(1-p_i(\epsilon))^{N-k-1}p_i(\epsilon)^{k-1}\frac{dp_i(\epsilon)^1}{d\epsilon} \end{aligned}$
令 $u=p_i(\epsilon)$ ，则 $du=dp_i(\epsilon)d\epsilon$
Beta函数定义为 $beta(p,q)=\int_0^{1}x^{p-1}(1-x)^{q-1}dx=\frac{(p-1)!(q-1)!}{(p+q-1)!=\frac{\frac{p+q}{pq}}{\tbinom{p+q}{p}}}$
因此，验证积分为1： (4)
$\begin{aligned} \int_0^{\infin}p_i(\epsilon)d\epsilon&=\int_0^1k\tbinom{N-1}{k}\mu^{k-1}(1-\mu)^{N-k-1}d\mu=k\tbinom{N-1}{k}beta(k,N-k)\\ &=k\tbinom{N-1}{k}\frac{N}{k(N-k)}\frac{1}{\tbinom{N}{k}}=1 \end{aligned}$

Lemma 1.

$X\sim Beta(\alpha, \beta),\mathbb{E}[log(X)]=\psi(\alpha)-\psi(\alpha+\beta),\psi()=\frac{\Gamma'()}{\Gamma()}$

证明： (5)
$\begin{aligned} E[log(X)]&=\int_0^1B(\alpha,\beta)^{-1}log(x)x^{\alpha-1}(1-x)^{\beta-1}dx\\ &=B(\alpha,\beta)^{-1}\int_0^1\frac{\partial}{\partial\alpha}x^{\alpha-1}(1-x)^{\beta-1}dx \end{aligned}$
交换微分和积分的顺序（类似于莱布尼兹法则）: (6)
$E[log(X)]=Beta(\alpha,\beta)^{-1}\frac{\partial}{\partial\alpha}Beta(\alpha,\beta)=\frac{\partial}{\partial\alpha}logB(\alpha,\beta)$
展开Beta函数： (7)
$=\frac{\partial}{\partial\alpha}[log\Gamma(\alpha)+log\Gamma(\beta)-log\Gamma(\alpha+\beta)]=\psi(\alpha)-\psi(\alpha+\beta)$

Proposition 2.

$E[logp_i(\epsilon)]=\psi(k)-\psi(N)$

使用Lemma 1可得。

熵的估计

假设在 $x_i$ 的邻域范围内，概率密度 $\mu(x)$ 是均匀的，则有 (8)
$p_i(\epsilon)\approx c_d\epsilon^d\mu(x_i)$
其中， $c_d$ 是n维的单位球体体积，定义为 $\int_{\Vert s \Vert<1/2}ds$ ,对于某些范数。

如果是2范数，则n维球体积公式 (9)
$V_n=\frac{\pi^{\frac{n}{2}}R^d}{\Gamma(1+\frac{n}{2})}$
$R = 1$ 时即为 $c_d$ ，再乘以 $R^d$ 。

因此，该区域概率密度 $\mu(x_i)$ 乘以体积 $c_d\epsilon^d$ ，得到该区域概率质量 $p_i(\epsilon)$

根据（8）,有 (10)
$\begin{aligned} E[-log\mu(X_i)]&\approx-E[logp_i(\epsilon)]+logc_d+dE[log\epsilon]\\ &=-\phi(k)+\phi(N)+logc_d+dE[log\epsilon]\\ &\approx-\phi(k)+\phi(N)+logc_d+\frac{d}{N}\sum_{i=1}^{N}log\epsilon_i \end{aligned}$
其中，最后的 $\epsilon_i$ 表示第 $i$ 个点到他第kth最近邻的距离。

KSG1估计器

首先， $H (X)$ 和 $H (Y)$ 可以直接使用（10）进行估计。

其次，使用(10)，联合熵的一个估计是 (11)：
$\hat{H}(X,Y)=-\psi(k)+\psi(N)+log(c_{dz})+\frac{dz}{N}\sum log\epsilon_i$
其中， $d z = d x + d y$ ，即Z的维度是x和y维度之和。且对于球的体积，KSG采用的最大范数进行的度量，也就是最大范数下N维空间球的体积，因此区域看上去是n维的立方体。

其中可以注意到： (12)
$c_{dz}=c_{dx}c_{dy}$
我思考时类比于正方体的底面积乘以高，这一点在下面（23）式中使用。

然而，使用固定的k估计的 $H (X) H (Y) H (X, Y)$ 直接对 $I (X, Y)$ 进行估计，存在问题。原因是，（10）式的偏差主要来自于对（8）式假设的均匀概率密度的违反，因此这三项估计都因此会存在偏差。而对于一个固定的k来说，联合空间上的第k近邻的距离要远于边际空间，因此实际上同样是kth近邻，空间的大小不同，因此概率密度违反（8）式假设的大小也不同，因此会导致偏差非常不同从而无法抵消（且观看下图图示，以1近邻为例）。
在这里插入图片描述

为了避免这一点，由于（10）对任意k都成立，因此不必使用固定的k。KSG1的方法是首先找到Joint上kth近邻，然后在X marginnal space上确定其距离 $\epsilon(i)$ ，在该距离内的点的个数作为 $n_x(i)$ ，则到joint上第kth点的距离正是 $\epsilon(i)$ ，而这个点在Xmarginal上则必然是第 $n_x(i)+1$ 个点。对Y marginal space用同样大小的距离 $\epsilon(i)$ 做同等操作。具体的情况如下图所示，此时Joint上是1th近邻，而此时x方向上产生5近邻，而y方向上是3近邻。
在这里插入图片描述

此时对于 $H (X)$ 来说，则有（考虑（2）式和（10）式中第一个约等号的蒙特卡洛近似，我未必写的非常非常详细，如果这里还需过于详细的介绍，则无需再看此文，而是重新补习统计相关知识）： (13)
$\hat{H}(X)=-\frac{1}{N}\sum_{i=1}^{N}\psi(n_x(i)+1)+\psi(N)+log(c_{dz})+\frac{dx}{N}\sum log\epsilon_i$
对于Y，则利用同样的距离找到 $n_y(i)+1$ 进行估计，导致Y的估计不准，因为y到其kth点的距离不是 $\frac{\epsilon}{2}$ 。

KSG2估计器

KSG2估计器更为复杂，其考虑在两个marginal的方向上采用不同的距离，此时图示joint的kth区域不再是超立方体区域，而是超长方体区域。首先根据joint上kth近邻的位置分别确定x marginal上的距离和y marginal上的距离（已经与KSG1不同，该方法确定x marginal上的距离后，直接将该距离应用于 y marginal上，因此是超正方体），此时点落在 $x_i\pm\epsilon_x(i)/2$ 和 $y_i\pm\epsilon_y(i)/2$ 之内。

此时存在两种情况，即 $\epsilon_x(i)$ 和 $\epsilon_y(i)$ 由同一点决定或由两点决定，具体见下图。
在这里插入图片描述

故考虑在将原先的超正方体概率替换为超长方体的概率，该超长方体的概率是两种情况的混合，则 (14)
$P_k(\epsilon_x,\epsilon_y)=P_k^b(\epsilon_x,\epsilon_y)+P_k^c(\epsilon_x,\epsilon_y)$
其中第一种时，有k-1个点在最内，中间的微分区域（逻辑与（3）相同）有1个点，外部有N-1-k个点在最外，若不失一般性的假设 $\epsilon_x<\epsilon_y$ ，该种状况发生的概率表达为 (15)
$\begin{aligned} P_k^b(\epsilon_x,\epsilon_y)&=\tbinom{N-1}{k}\tbinom{k}{k-1}q_i^{k-1}(\frac{\partial ^2 q_i}{\partial\epsilon_x\partial\epsilon_y})(1-p_i(\epsilon_y))^{N-k-1}\\ &=k\tbinom{N-1}{k}q_i^{k-1}(\frac{\partial ^2 q_i}{\partial\epsilon_x\partial\epsilon_y})(1-p_i(\epsilon_y))^{N-k-1} \end{aligned}$
原文的表述将 $q_i$ 吸进了偏导项中，其中 $q_i$ 表示的是上述超长方体的区域，而 $p_i$ 则仍然是由较大的 $\epsilon_y$ 确定的超正方体区域。使用 $p_i$ 是必须的，因为使用了最大范数，即最远的点的距离，保证了如果该正方形内的点都在矩形中，反之则不行。

对于第二种情况，此时有k-2个点在最内，而2个点在微分区域，N-k-1个点在外围区域，因此此种状况发生的概率为 (16)
$\begin{aligned} P_k^c(\epsilon_x,\epsilon_y)&=\tbinom{N-1}{k}\tbinom{k}{k-2}q_i^{k-2}(\frac{\partial ^2 q_i^2}{\partial\epsilon_x\partial\epsilon_y})(1-p_i(\epsilon_y))^{N-k-1}\\ &=\tbinom{N-1}{k}\tbinom{k}{k-2}2q_iq_i^{k-2}(\frac{\partial ^2 q_i}{\partial\epsilon_x\partial\epsilon_y})(1-p_i(\epsilon_y))^{N-k-1}\\ &=k(k-1)\tbinom{N-1}{k}q_i^{k-1}(\frac{\partial ^2 q_i}{\partial\epsilon_x\partial\epsilon_y})(1-p_i(\epsilon_y))^{N-k-1} \end{aligned}$
同时注意到，由于X和Y上是独立的，根据 (17)
$u(\frac{\partial q}{\partial x\partial y})=(\frac{\partial q}{\partial x})(\frac{\partial q}{\partial y})$
故（16）进一步拆解为 (18)
$k(k-1)\tbinom{N-1}{k}q_i^{k-2}(\frac{\partial q}{\partial x})(\frac{\partial q}{\partial y})(1-p_i(\epsilon_y))^{N-k-1}$
验证 $P_k(\epsilon_x,\epsilon_y)=P_k^b(\epsilon_x,\epsilon_y)+P_k^c(\epsilon_x,\epsilon_y)$ 的积分为1. (19)
$\begin{aligned} &\int_0^{\infin}\int_0^{\epsilon_y}P_k(\epsilon_x,\epsilon_y)d\epsilon_xd\epsilon_y\\ &=\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}\int_0^{\epsilon_y}\frac{\partial}{\partial x}({kq_i^{k-1}(\frac{\partial q_i}{\partial\epsilon_y}}))d\epsilon_xd\epsilon_y\\ &=k\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}p_i(\epsilon_y)^{k-1}d\epsilon_y\\ &=k\tbinom{N-1}{k}Beta(k,N-k)=1 \end{aligned}$
下面证明 $E[logq_i(\epsilon_x,\epsilon_y)]=\psi(k)-\psi(N)-1/k$ ，首先在第二个积分号进行分部积分： (20)
$\begin{aligned} E[logq_i(\epsilon_x,\epsilon_y)]&=\int_0^{\infin}\int_0^{\epsilon_y}log(q_i)P_k(\epsilon_x,\epsilon_y)d\epsilon_xd\epsilon_y\\ &=\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}\int_0^{\epsilon_y}\frac{\partial}{\partial x}({kq_i^{k-1}(\frac{\partial q_i}{\partial\epsilon_y}}))d\epsilon_xd\epsilon_y\\ &=\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}[log(q_i)kq_i^{k-1}(\frac{\partial q_i}{\partial \epsilon_y})|_{0}^{\epsilon_y}\\&-\int_0^{\epsilon_y}\{kq_i^{k-2}\frac{\partial q_i}{\partial \epsilon_x}(\frac{\partial q_i}{\partial \epsilon_y})\}d\epsilon_x]d\epsilon_y\\ \end{aligned}$
由最开始的Lemma1可知 $X\sim Beta(\alpha, \beta),\mathbb{E}[log(X)]=\psi(\alpha)-\psi(\alpha+\beta),\psi()=\frac{\Gamma'()}{\Gamma()}$

可以观察到第一部分： (21)
$\begin{aligned} &=\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}[log(q_i)kq_i^{k-1}(\frac{\partial q_i}{\partial \epsilon_y})|_{0}^{\epsilon_y}d\epsilon_xd\epsilon_y\\ &=\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}log(p_i(\epsilon_y))kp_i(\epsilon_y)^{k-1}(\frac{\partial p_i(\epsilon_y)}{\partial \epsilon_y})d\epsilon_y\\ &=\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}log(p_i(\epsilon_y))kp_i(\epsilon_y)^{k-1}(\frac{\partial p_i(\epsilon_y)}{\partial \epsilon_y})d\epsilon_y\\ &=\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}log(p_i(\epsilon_y))kp_i(\epsilon_y)^{k-1}dp_i(\epsilon_y)\\ &=\tbinom{N-1}{k}\int_0^{1}(1-u)^{N-k-1}log(u)ku^{k-1}du=\psi(k)-\psi(N) \end{aligned}$
其中，注意到 $q_i(0)=0$ 和当 $\epsilon_x=\epsilon_y$ 时有 $q(\epsilon_y,\epsilon_y)$ 为一个超立方体，恰是 $p_i(\epsilon_y)$ ，最后的部分由前面Prososition 2.的结论得出。

第二部分： (22)
$\begin{aligned} &-\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}(\int_0^{\epsilon_y}\{kq_i^{k-2}\frac{\partial q_i}{\partial \epsilon_x}(\frac{\partial q_i}{\partial \epsilon_y})\}d\epsilon_x)d\epsilon_y\\ &=-\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}(\int_0^{\epsilon_y}\{kq_i^{k-2}\frac{\partial q_i}{\partial \epsilon_x}(\frac{\partial q_i}{\partial \epsilon_y})\}d\epsilon_x)d\epsilon_y \end{aligned}$
由局部均匀假设： (23)
$\begin{aligned} &q_i(\epsilon_x,\epsilon_y)\approx\mu(x_i,y_i)c_{dx}c_{dy}\epsilon_x^{dx}\epsilon_y^{dy}\\ \rightarrow&\frac{\partial}{\partial \epsilon_y}q_i(\epsilon_x,\epsilon_y)=\mu(x_i,y_i)c_{dx}c_{dy}\epsilon_x^{dx}dy\epsilon_y^{dy-1}=q_i(\epsilon_x,\epsilon_y)\frac{dy}{\epsilon_y} \end{aligned}$
并且由刚才分部积分时提及的当 $\epsilon_x=\epsilon_y$ 时有 $q(\epsilon_y,\epsilon_y)$ 为一个超立方体，因此 (24)
$\frac{\partial}{\partial \epsilon_y}q_i(\epsilon_x,\epsilon_y)|_{\epsilon_x=\epsilon_y}=p_i(\epsilon_y)\frac{dy}{\epsilon_y}$

则（22）继续简化有： (25)
$\begin{aligned} &=-\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}(\int_0^{\epsilon_y}\{kq_i^{k-2}\frac{\partial q_i}{\partial \epsilon_x}(\frac{\partial q_i}{\partial \epsilon_y})\}d\epsilon_x)d\epsilon_y\\ &=-\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}(\int_0^{\epsilon_y}\{kq_i^{k-2}\frac{\partial q_i}{\partial \epsilon_x}(q_i\frac{dy}{\epsilon_y})\}d\epsilon_x)d\epsilon_y\\ &=-\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}(\frac{dy}{\epsilon_y})(\int_0^{\epsilon_y}\{kq_i^{k-1}\frac{\partial q_i}{\partial \epsilon_x}\}d\epsilon_x)d\epsilon_y\\ &=-\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}(\frac{dy}{\epsilon_y})p_i(\epsilon_y)^k)d\epsilon_y\\ &=-\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}(\frac{dy}{\epsilon_y})p_i(\epsilon_y)^k)d\epsilon_y \end{aligned}$
其中第三个等号注意到了 $\frac{dy}{\epsilon_y}$ 是常数，注意 $d y$ 是 $Y$ 的维度不是微分。第四个等号是将 $kq_i^{k-1}$ 吸收到微分，再由积分所得。

进一步根据（24）有：(26)
$\frac{d}{d\epsilon_y}p_i(\epsilon_y)=\frac{d}{\partial \epsilon_y}q_i(\epsilon_x,\epsilon_y)|_{\epsilon_x=\epsilon_y}=q_i(\epsilon_x,\epsilon_y)|_{\epsilon_x=\epsilon_y}\frac{dy}{d\epsilon_y}=p_i(\epsilon_y)\frac{dy}{\epsilon_y}$
因此，（25）简化为(26)
$\begin{aligned} &=-\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}(\frac{dy}{\epsilon_y})p_i(\epsilon_y)^k)d\epsilon_y\\ &=-\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}p_i(\epsilon_y)^{k-1}\frac{\partial}{\partial y}p_i(\epsilon_y)d\epsilon_y\\ &=-\tbinom{N-1}{k}\int_0^{\infin}(1-p_i(\epsilon_y))^{N-k-1}p_i(\epsilon_y)^{k-1}dp_i=-\tbinom{N-1}{k}Beta(k,N-k)\\ &=-\frac{1}{k} \end{aligned}$
联立（21）和（26），得到(27)
$E[logq_i(\epsilon_x,\epsilon_y)]=\psi(k)-\psi(N)-\frac{1}{k}$
进一步的，根据公式（23），自然可得(28)
$\begin{aligned} log\mu(x_i,y_i)\approx logq_i(\epsilon_x,\epsilon_y)-log(c_{dx}c_{dy})-dxlog\epsilon_x-dylog\epsilon_y \end{aligned}$
求期望有：(29)
$E[log\mu(x_i,x_y)]=\\\psi(k)-\psi(N)-\frac{1}{k}-log(c_{dx}c_{dy})-dx\frac{1}{N}\sum log\epsilon_x-dy\frac{1}{N}\sum log\epsilon_y$
再进一步利用（1）式可得KSG2估计器：(30)
$I^{(2)}(X,Y)=\psi(k)+\psi(N)-\frac{1}{k}-\frac{1}{N}\sum\psi(n_x)-\frac{1}{N}\sum\psi(n_y)$