Estimating Mutual Information for Discrete-Continuous Mixtures

Introduction

A fundamental quantity of interest in machine learning is mutual information (MI), which characterizes the shared information between a pair of random variables $(X,Y)$. MI obeys several appealing properties including the data-processing inequality, invariance under one-to-one transformations and the chain rule , which led to a wide use in canonical tasks such as classification , clustering and feature selection . Mutual information also emerges as the “correct” quantity in several graphical model inference problems (e.g., the Chow-Liu tree and conditional independence testing ). MI is also pervasively used in many data science application domains, such as sociology , computational biology , and computational neuroscience .

An important problem in any of these applications is to estimate mutual information effectively from samples. While mutual information has been the de facto measure of information in several applications for decades, the estimation of mutual information from samples remains an active research problem. Recently, there has been a resurgence of interest in entropy and mutual information estimators, on both the theoretical as well as practical fronts .

The previous estimators focus on either of two cases – the data is either purely discrete or purely continuous. In these special cases, the mutual information can be calculated based on the three (differential) entropies of $X$, $Y$ and $(X,Y)$. We term estimators based on this principle as $3H$-estimators (since they estimate three entropy terms), and a majority of previous estimators fall under this category .

In practical downstream applications, we often have to deal with a mixture of continuous and discrete random variables. Random variables can be mixed in several ways. First, one random variable can be discrete whereas the other is continuous. For example, we want to measure the strength of relationship between children’s age and height, here age $X$ is discrete and height $Y$ is continuous. Secondly, a single scalar random variable itself can be a mixture of discrete and continuous components. For example, consider $X$ taking a zero-inflated-Gaussian distribution, which takes value with probability $p$ and is a Gaussian distribution with mean $\mu$ with probability $1-p$. This distribution has both a discrete component as well as a component with density. Finally, $X$ and / or $Y$ can be high dimensional vector, each of whose components may be discrete, continuous or mixed.

Problem Formation

In this section, we define mutual information for general distributions as follows (e.g., ).

Let $P_{XY}$ be a probability measure on the space $\mathcal{X}\times\mathcal{Y}$, where $\mathcal{X}$ and $\mathcal{Y}$ are both Euclidean spaces. For any measurable set $A\subseteq\mathcal{X}$ and $B\subseteq\mathcal{Y}$, define $P_{X}(A)=P_{XY}(A\times\mathcal{Y})$ and $P_{Y}(B)=P_{XY}(\mathcal{X}\times B)$. Let $P_{X}P_{Y}$ be the product measure $P_{X}\times P_{Y}$. If $P_{XY}$ is absolutely continuous w.r.t. $P_{X}P_{Y}$, then the mutual information $I(X;Y)$ of $P_{XY}$ is defined as

where $\frac{dP_{XY}}{dP_{X}P_{Y}}$ is the Radon-Nikodym derivative.

Notice that this general definition includes the following cases of mixtures: (1) $X$ is discrete and $Y$ is continuous (or vice versa); (2) $X$ or $Y$ has many components each, where some components are discrete and some are continuous; (3) $X$ or $Y$ or their joint distribution is a mixture of continuous and discrete distributions.

Estimators of Mutual Information

The estimation problem is quite different depending on whether the underlying distribution is discrete, continuous or mixed. As pointed out earlier, most existing estimators for mutual information are based on the $3H$ principle: they estimate the three entropy terms first. This $3H$ principle can be applied only in the purely discrete or purely continuous case.

Discrete data: For entropy estimation of a discrete variable $X$, the straightforward approach to plug-in the estimated probabilities $\hat{p}_{X}(x)$ into the formula for entropy has been shown to be suboptimal . Novel entropy estimators with sub-linear sample complexity have been proposed . MI estimation can then be performed using the $3H$ principle, and such an approach is shown to be worst-case optimal for mutual-information estimation .

Continuous data: There are several estimators for differential entropy of continuous random variables, which have been exploited in a $3H$ principle to calculate the mutual information . One family of entropy estimators are based on kernel density estimators followed by re-substitution estimation. An alternate family of entropy estimators is based on $k$-Nearest Neighbor ($k$-NN) estimates, beginning with the pioneering work of Kozachenko and Leonenko (the so-called KL estimator). Recent progress involves an inspired mixture of an ensemble of kernel and $k$-NN estimators . Exponential concentration bounds under certain conditions are in .

Mixed Random Variables: Since the entropies themselves may not be well defined for mixed random variables, there is no direct way to apply the $3H$ principle. However, once the data is quantized, this principle can be applied in the discrete domain. That mutual information in arbitrary measure spaces can indeed be computed as a maximum over quantization is a classical result . However, the choice of quantization is complicated and while some quantization schemes are known to be consistent when there is a joint density , the mixed case is complex. Estimator of the average of Radon-Nikodym derivative $dP/dQ$ has been studied in . Very recent work generalizing the ensemble entropy estimator when some components are discrete and others continuous is in .

Beyond $3H$ estimation: In an inspired work proposed a direct method for estimating mutual information (KSG estimator) when the variables have a joint density. The estimator starts with the $3H$ estimator based on differential entropy estimates based on the $k$-NN estimates, and employs a heuristic to couple the estimates in order to improve the estimator. While the original paper did not contain any theoretical proof, even of consistency, its excellent practical performance has encouraged widespread adoption. Recent work has established the consistency of this estimator along with its convergence rate. Further, recent works involving a combination of kernel density estimators and $k$-NN methods have been proposed to further improve the KSG estimator. extends the KSG estimator to the case when one variable is discrete and another is scalar continuous.

None of these works consider a case even if one of the components has a mixture of continuous and discrete distribution, let alone for general probability distributions. There are two generic options: (1) one can add small independent noise on each sample to break the multiple samples and apply a continuous valued MI estimator (like KSG), or (2) quantize and apply discrete MI estimators but the performance for high-dimensional case is poor. These form baselines to compare against in our detailed simulations.

2 Mixed Regime

We first examine the behavior of other estimators in the mixed regime, before proceeding to develop our estimator. Let us consider the case when $X$ is discrete (but real valued) and $Y$ possesses a density. In this case, we will examine the consequence of using the $3H$ principle, with differential entropy estimated by the $k$-nearest neighbors. To do this, fix a parameter $k$, that determines the number of neighbors and let $\rho_{i,x}$, $\rho_{i,y}$ and $\rho_{i,xy}$ denote the distance of the $k$-nearest neighbor of $X_{i}$, $Y_{i}$ and $(X_{i},Y_{i})$, respectively. Then

where $\psi(\cdot)$ is the digamma function and $a(\cdot)=\log(\cdot)-\psi(\cdot)$. In the case that $X$ is discrete and $Y$ has a density, $I_{\rm 3H}(X;Y)=-\infty+a-b=-\infty$, which is clearly wrong.

The basic idea of the KSG estimator is to ensure that the $\rho$ is the same for both $x$, $y$ and $(x,y)$ and the difference is instead in the number of nearest neighbors. Let $n_{x,i}$ be the number of samples of $X_{i}$’s within distance $\rho_{i,xy}$ and $n_{y,i}$ be the number of samples of $Y_{i}$’s within distance $\rho_{i,xy}$. Then the KSG estimator is given by $\widehat{I}_{KSG}^{(N)}\equiv\frac{1}{N}\sum_{i=1}^{N}\left(\,\psi(k)+\log(N)-\log(n_{x,i}+1)-\log(n_{y,i}+1)\,\right)$ where $\psi(\cdot)$ is the digamma function.

In the case of $X$ being discrete and $Y$ being continuous, it turns out that the KSG estimator does not blow up (unlike the $3H$ estimator), since the distances do not go to zero. However, in the mixed case, the estimator has a non-trivial bias due to discrete points and is no longer consistent.

3 Proposed Estimator

We propose the following estimator for general probability distributions, inspired by the KSG estimator. The intuition is as follows. First notice that MI is the average of the logarithm of Radon-Nikodym derivative, so we compute the Radon-Nikodym derivative for each sample $i$ and take the empirical average. The re-substitution estimator for MI is then given as follows: $\widehat{I}(X;Y)\equiv\frac{1}{n}\sum_{i=1}^{n}\log\left(\frac{dP_{XY}}{dP_{X}P_{Y}}\right)_{(x_{i},y_{i})}.$ The basic idea behind our estimate of the Radon-Nikodym derivative at each sample point is as follows:

When the point is discrete (which can be detected by checking if the $k$-nearest neighbor distance of data $i$ is zero), then we can assert that data $i$ is in a discrete component, and we can use plug-in estimator for Radon-Nikodym derivative.

If the point is such that there is a joint density (locally), the KSG estimator suggests a natural idea: fix the radius and estimate the Radon-Nikodym derivative by $\left(\psi(k)+\log(N)-\log(n_{x,i}+1)-\log(n_{y,i}+1)\right)$.

If $k$-nearest neighbor distance is not zero, then it may be either purely continuous or mixed. But we show below that the method for purely continuous is also applicable for mixed.

Proof of Consistency

We show that under certain technical conditions on the joint probability measure, the proposed estimator is consistent. We begin with the following definitions. Let $f(x,y)=dP_{XY}/dP_{X}P_{Y}$ denote the Radon-Nikodym derivative and define

$k$ is chosen to be a function of $N$ such that $k_{N}\to\infty$ and $k_{N}\log N/N\to 0$ as $N\to\infty$.

The set of discrete points $\{(x,y):P_{XY}(x,y,0)>0\}$ is finite.

\int_{\mathcal{X}\times\mathcal{Y}}\big{|}\,\log\frac{dP_{XY}}{dP_{X}P_{Y}}\,\big{|}\,dP_{XY}<+\infty.

Notice that the assumptions are satisfied whenever (1) the distribution is (finitely) discrete; (2) the distribution is continuous; (3) some dimensions are (countably) discrete and some dimensions are continuous; (4) a (finite) mixture of the previous cases. Most real world data can be covered by these cases. A sketch of the proof is below with the full proof in the supplementary material.

(Sketch) We start with an explicit form of the Radon-Nikodym derivative $dP_{XY}/(dP_{X}P_{Y})$.

For almost every $(x,y)\in\mathcal{X}\times\mathcal{Y}$, we have

$\Omega_{2}=\{(x,y):f(x,y)>0,P_{XY}(x,y,0)>0\}\,;$

$\Omega_{3}=\{(x,y):f(x,y)>0,P_{XY}(x,y,0)=0\}\;.$

For $(x,y)\in\Omega_{3}$, it can be viewed as a continuous part. We use the similar proof technique as to prove that the mean of estimate $\xi_{1}$ is closed to $\log f(x,y)$.

The following theorem bounds the variance of the proposed estimator.

$(k_{N}\log N)^{2}/N\to 0$ as $N\to\infty$.

(Sketch) We use the Efron-Stein inequality to bound the variance of the estimator. For simplicity, let $\widehat{I}^{(N)}(Z)$ be the estimate based on original samples $\{Z_{1},Z_{2},\dots,Z_{N}\}$, where $Z_{i}=(X_{i},Y_{i})$, and $\widehat{I}^{(N)}(Z_{\setminus j})$ is the estimate from $\{Z_{1},\dots,Z_{j-1},Z_{j+1},\dots,Z_{N}\}$. Then a certain version of Efron-Stein inequality states that: {\textrm{ Var }}\left[\,\widehat{I}^{(N)}(Z)\,\right]\leq 2\sum_{j=1}^{N}\left(\,\sup_{Z_{1},\dots,Z_{N}}\Big{|}\,\widehat{I}^{(N)}(Z)-\widehat{I}^{(N)}(Z_{\setminus j})\,\Big{|}\,\right)^{2}\;. Now recall that

To upper bound the difference $|\,\xi_{i}(Z)-\xi_{i}(Z_{\setminus j})\,|$ created by eliminating sample $Z_{j}$ for different $i$ ’s we consider three different cases: (1) $i=j$; (2) $\rho_{k,i}=0$; (3) $\rho_{k,i}>0$, and conclude that $\sum_{i=1}^{N}|\,\xi_{i}(Z)-\xi_{i}(Z_{\setminus j})\,|\leq O(k\log N)$ for all $Z_{i}$’s. The detail of the case study is in Section. B in the supplementary material. Plug it into Efron-Stein inequality, we obtain:

By Assumption 6, we have $\lim_{N\to\infty}{\textrm{ Var }}\left[\,\widehat{I}^{(N)}(Z)\,\right]=0$. ∎

Simulations

We evaluate the performance of our estimator in a variety of (synthetic and real-world) experiments.

Experiment I. $(X,Y)$ is a mixture of one continuous distribution and one discrete distribution. The continuous distribution is jointly Gaussian with zero mean and covariance $\Sigma=\begin{pmatrix}1&0.9\\ 0.9&1\end{pmatrix}$, and the discrete distribution is $P(X=1,Y=1)=P(X=-1,Y=-1)=0.45$ and $P(X=1,Y=-1)=P(X=-1,1)=0.05$. These two distributions are mixed with equal probability. The scatter plot of a set of samples from this distribution is shown in the left panel of Figure. 1, where the red squares denote multiple samples from the discrete distribution. For all synthetic experiments, we compare our proposed estimator with a (fixed) partitioning estimator, an adaptive partitioning estimator implemented by , the KSG estimator and noisy KSG estimator (by adding Gaussian noise $N(0,\sigma^{2}I)$ on each sample to transform all mixed distributions into continuous one). We plot the mean squared error versus number of samples in Figure 2. The mean squared error is averaged over 250 independent trials.

The KSG estimator is entirely misled by the discrete samples as expected. The noisy KSG estimator performs better but the added noise causes the estimate to degrade. In this experiment, the estimate is less sensitive to the noise added and the line is indistinguishable with the line for KSG. The partitioning and adaptive partitioning method quantizes all samples, resulting in an extra quantization error. Note that only the proposed estimator has error decreasing with the sample size.

Experiment II. $X$ is a discrete random variable and $Y$ is a continuous random variable. $X$ is uniformly distributed over integers $\{0,1,\dots,m-1\}$ and $Y$ is uniformly distributed over the range $[X,X+2]$ for a given $X$. The ground truth $I(X;Y)=\log(m)-(m-1)\log(2)/m$. We choose $m=5$ and a scatter plot of a set of samples is in the right panel of Figure. 1. Notice that in this case (and the following experiments) our proposed estimator degenerates to KSG if the hyper parameter $k$ is chosen the same, hence KSG is not plotted. In this experiment our proposed estimator outperforms other methods.

Experiment III. Higher dimensional mixture. Let $(X_{1},Y_{1})$ and $(Y_{2},X_{2})$ have the same joint distribution as in experiment II and independent of each other. We evaluate the mutual information between $X=(X_{1},X_{2})$ and $Y=(Y_{1},Y_{2})$. Then ground truth $I(X;Y)=2(\log(m)-(m-1)\log(2)/m)$. We also consider $X=(X_{1},X_{2},X_{3})$ and $Y=(Y_{1},Y_{2},Y_{3})$ where $(X_{3},Y_{3})$ have the same joint distribution as in experiment II and independent of $(X_{1},Y_{1}),(X_{2},Y_{2})$. The ground truth $I(X;Y)=3(\log(m)-(m-1)\log(2)/m)$. The adaptive partitioning algorithm works only for one-dimensional $X$ and $Y$ and is not compared here.

We can see that the performance of partitioning estimator is very bad because the number of partitions grows exponentially with dimension. Proposed algorithm suffers less from the curse of dimensionality. For the right figure, noisy KSG method has smaller error, but we point out that it is unstable with respect to the noise level added: as the noise level is varied from $\sigma=0.5$ to $\sigma=0.7$ and the performance varies significantly (far from convergence).

Experiment IV. Zero-inflated Poissonization. Here $X\sim{\rm Exp}(1)$ is a standard exponential random variable, and $Y$ is zero-inflated Poissonization of $X$, i.e., $Y=0$ with probability $p$ and $Y\sim{\rm Poisson}(x)$ given $X=x$ with probability $1-p$. Here the ground truth is $I(X;Y)=(1-p)(2\log 2-\gamma-\sum_{k=1}^{\infty}\log k\cdot 2^{-k})\approx(1-p)0.3012$, where $\gamma$ is Euler-Mascheroni constant. We repeat the experiment for no zero-inflation ($p=0$) and for $p=15\%$. We find that the proposed estimator is comparable to adaptive partitioning for no zero-inflation and outperforms others for 15% zero-inflation.

We conclude that our proposed estimator is consistent for all these four experiments, and the mean squared error is always the best or comparable to the best. Other estimators are either not consistent or have large mean squared error for at least one experiment.

Gene regulatory network inference. Gene expressions form a rich source of data from which to infer gene regulatory networks; it is now possible to sequence gene expression data from single cells using a technology called single-cell RNA-sequencing . However, this technology has a problem called dropout, which implies that sometimes, even when the gene is present it is not sequenced . While we tested our algorithm on real single-cell RNA-seq dataset, it is hard to establish the ground truth on these datasets. Instead we resorted to a challenge dataset for reconstructing regulatory networks, called the DREAM5 challenge . The simulated (insilico) version of this dataset contains gene expression for 20 genes with 660 data point containing various perturbations. The goal is to reconstruct the true network between the various genes. We used mutual information as the test statistic in order to obtain AUROC for various methods. While the dataset did not have any dropouts, in order to simulate the effect of dropouts in real data, we simulated various levels of dropout and compared the AUROC (area under ROC) of different algorithms in the right of Figure 3 where we find the proposed algorithm to outperform the competing ones.

Acknowledgement

We thank Arman Rahimzamani and Himanshu Asnani for their constructive comments on the proofs of the lemmas, especially for the proof of Lemma A.2.

Appendix

Appendix A Proof of Theorem 1

To prove the asymptotic unbiasedness of the estimator, we need to write the Radon-Nikodym derivative in an explicit form. The following lemma gives the explicit form of $\frac{dP_{XY}}{dP_{X}P_{Y}}$.

For almost every $(x,y)\in\mathcal{X}\times\mathcal{Y}$, $\frac{dP_{XY}}{dP_{X}P_{Y}}=f(x,y)=\lim_{r\to 0}\frac{P_{XY}(x,y,r)}{P_{X}(x,r)P_{Y}(y,r)}$.

$\Omega_{2}=\{(x,y):f(x,y)>0,P_{XY}(x,y,0)>0\}\,;$

$\Omega_{3}=\{(x,y):f(x,y)>0,P_{XY}(x,y,0)=0\}\;.$

$(x,y)\in\Omega_{1}$: In this case, we will show that $\Omega_{1}$ has zero probability with respect to $P_{XY}$.

First, the probability of $\rho_{k,1}>0$ is upper bounded by:

If $X$ is distributed as $\text{Bino}(N,p)$ and $m\geq 0$ , then:

By Assumption 2, $k/N\rightarrow 0$ as $N\rightarrow\infty$, then $(Np-k)^{2}/N=N(p-k/N)^{2}\rightarrow\infty$, So for sufficiently large $N$, the RHS of Lemma A.2 is upper bounded by $\max\{\frac{C_{1}m}{Np},\frac{C_{2}}{Np}\}\leq\frac{C(m+1)}{Np}$, where $C=\max\{C_{1},C_{2}\}$ is some constant not depends on $N$. Therefore, by applying Lemma A.2 with $m=1$, the first term of (14) is bounded by:

Combine with the case that $\rho_{i,xy}>0$, we obtain that:

where the first term comes from triangle inequality and the fact that $|\xi_{1}|\leq 2\log N$. Integrating over $\Omega_{2}$, we have:

where $\mu$ denotes counting measure. By Assumption 1, $k$ goes to infinity as $N$ goes to infinity, so $1/k$ vanishes as $N$ increases. By Assumption 1 and 2, $k/N$ goes to 0 and $\Omega_{2}$ has finite counting measure, so the second term also vanishes. Since $\Omega_{2}$ has finite counting measure, so $\inf_{(x,y)\in\Omega_{2}}P_{XY}(x,y,0)=\epsilon>0$. By Assumption 3, $\int_{\Omega_{2}}|\,\log f(x,y)\,|dP_{XY}<+\infty$. Therefore, for sufficiently large $N$, the first term also vanishes. Therefore,

$(x,y)\in\Omega_{3}$: In this case, $P_{XY}(x,y,r)$ is a monotonic function of $r$ such that $P_{XY}(x,y,0)=0$ and $\lim_{r\to\infty}P_{XY}(x,y,r)=1$. Hence, we can view $\log\left(\,P_{XY}(x,y,r)/P_{X}(x,r)P_{Y}(y,r)\,\right)$ as a function of $P_{XY}(x,y,r)$, and it converges to $\log f(x,y)$ as $P_{XY}(x,y,r)\to 0$, for almost every $(x,y)$. Since $P_{XY}(\Omega_{3})\leq 1<+\infty$ and $\int_{\Omega_{3}}|\log f(x,y)|dP_{XY}<+\infty$. Then by Egoroff’s Theorem, for any $\epsilon>0$, there exists a subset $E\subseteq\Omega_{3}$ with $P_{XY}(E)<\epsilon$ and $\int_{E}|\log f(x,y)|dP_{XY}<\epsilon$, such that $\log\left(\,P_{XY}(x,y,r)/P_{X}(x,r)P_{Y}(y,r)\,\right)$ converges as $P_{XY}(x,y,r)\to 0$, uniformly on $\Omega_{3}\setminus E$. For $(x,y)\in E$, notice that $|\xi_{1}|\leq 2\log N$, so we have:

here $F_{\rho_{k,1}}(r)$ is the CDF of the $k$-nearest neighbor distance $\rho_{k,1}$, given $(X,Y)=(x,y)$. By results of order statistics, its derivative with respect to $P_{XY}(x,y,r)$ is given by:

Now we consider the four terms separately. For (25), since $\log\left(\,P_{XY}(x,y,r)/P_{X}(x,r)P_{Y}(y,r)\,\right)$ converges as $P_{XY}(x,y,r)\to 0$, uniformly on $\Omega_{3}\setminus E$. So for every $(x,y)\in\Omega_{3}\setminus E$, there exists an $r_{N}$ such that $P_{XY}(x,y,r_{N})=4k\log N/N$ and $|\log\left(\,P_{XY}(x,y,r)/P_{X}(x,r)P_{Y}(y,r)\,\right)-\log f(x,y)|<\delta_{N}$ for every $r\leq r_{N}$. Here $r_{N}$ may depend on $(x,y)$, but $\delta_{N}$ does not depend on $(x,y)$ and $\lim_{N\to\infty}\delta_{N}=0$. Therefore, (25) is upper bounded by:

for sufficiently large $N$ such that $N-k>N/2$. Therefore, (25) is upper bounded by

For (25), we simply plug in $F_{\rho_{k,1}}(r)$ and integrate over $P_{XY}(x,y,r)$ and obtain

where we use the fact that $\psi(k)-\psi(N)=\frac{(N-1)!}{(k-1)!(N-k-1)!}\int_{t=0}^{1}(\log t)t^{k-1}(1-t)^{N-k-1}dt$. Notice that $\psi(N)<\log N$ and $\lim_{N\to 0}(\psi(N)-\log N)=0$.

Now we deal with (25) and (25). The following lemmas establish the distribution of $n_{x,1}$ and $n_{y,1}$ given $(X,Y)=(x,y)$ and $\rho_{k,1}=r>0$.

Given $(X,Y)=(x,y)$ and $\rho_{k,1}=r>0$, then $n_{x,1}-k$ is distributed as ${\rm Bino}(N-k-1,\frac{P_{X}(x,r)-P_{XY}(x,y,r)}{1-P_{XY}(x,y,r)})$; $n_{y,1}-k$ is distributed as ${\rm Bino}(N-k-1,\frac{P_{Y}(y,r)-P_{XY}(x,y,r)}{1-P_{XY}(x,y,r)})$.

The following lemma is useful to establish the upper bound for (25) and (25).

Now we are ready to upper bound (25). First, we rewrite the term (25) as:

For (32), by the fact that $\log(x/y)\leq(x-y)/y$ for all $x,y>0$ and Cauchy-Schwarz inequality, we have the following:

Notice that $P_{X}(x,r)\geq P_{XY}(x,y,r)$ for all $r$, so the second expectation is always no larger than 1. For the first expectation, we plug in $F_{\rho_{k,1}}(r)$ and integrate over $P_{XY}(x,y,r)$, let $t=P_{XY}(x,y,r)$ and observe,

For sufficiently large $N$ and $k$, it is upper bounded by $C_{1}(1/N+1/k)$ for some constant $C_{1}>0$. Therefore,

Similarly, by using the fact that $\log(x/y)>(x-y)/x$ and Cauchy-Schwarz inequality again, we conclude that there are some constant $C_{2}>0$ such that

Therefore, by combining (33), (36) and (37), we obtain

where $C^{\prime}=\max\{C_{1},C_{2}\}$. Since (25) and (25) are symmetric, the same upper bound (38) also applies to (25). Combine (29), (30) and (38), we have

for every $(x,y)\in\Omega_{3}\setminus E$. By integration over $\Omega_{3}\setminus E$, we have

By Assumption 1, $k$ increases as $N\to\infty$. By Assumption 3, $\int_{\mathcal{X}\times\mathcal{Y}}|\log f(x,y)|dP_{XY}<+\infty$. Therefore, this quantity vanishes as $N\to\infty$. Combining with the case that $(x,y)\in E$, we have

The proof of this lemma utilizes the Lebesgue-Besicovitch differentiation theorem [11, Theorem 1.32], stated below

For our lemma, let $f=\frac{dP_{XY}}{dP_{X}P_{Y}}$ and $\mu=P_{X}P_{Y}$. Since $\mu$ is a probability measure, it is a Radon measure of Euclidean space. Also, since $\int_{\mathcal{X}\times\mathcal{Y}}|f|d\mu=1$, so $f$ is globally integrable, hence locally integrable with respect to $\mu$. So the conditions of Lebesgue-Besicovitch differentiation theorem are satisfied, so

A.2 Proof of Lemma A.2

since $\zeta\geq\min\left\{x,x_{0}\right\}=\min\left\{x,Np\right\}$, we have:

Now taking the conditional expectations from both sides, we have:

In which we used the fact that $2i^{2}\geq(i+1)(i+2)$ for $i\geq 4$, and $(N+2)p\geq 4p$ for $N\geq 2$. Plugging it into Equation 60, we have:

A.3 Proof of Lemma A.3

Now we deal with the case that $\rho_{k,1}=r>0$. Given that $(X_{1},Y_{1})=(x,y)$ and $\rho_{k,1}=r>0$, we sort the samples $\{(X_{i},Y_{i})\}_{i=2}^{N}$ by their distance to $(x,y)$ defined as $d_{i}=\max\{\|X_{i}-x\|,\|Y_{i}-y\|\}$. To avoid the case that two samples have identical distance, we introduce a set of random variables $\{Z_{i}\}_{i=2}^{N}$ i.i.d. samples from ${\rm Unif}$ and define a comparison operator $\prec$ as:

Since for any $i\neq j$, the probability that $Z_{i}=Z_{j}$ is zero, so we can have either $i\prec j$ or $i\succ j$ with probability 1. Now let $\{2,3,\dots,N\}=S\cup\{j\}\cup T$ be a partition of the indices with $\left|S\right|=k-1$ and $\left|T\right|=N-k-1$. Define an event $\mathcal{A}_{S,j,T}$ associated to the partition as:

A.4 Proof of Lemma A.4

(i) $Np\geq m$. In this case, for any $x$, by applying Taylor’s theorem around $x_{0}=Np+m$, there exists $\zeta$ between $x$ and $x_{0}$ such that

By noticing that $\zeta\geq\min\{x,x_{0}\}=\min\{x,Np+m\}$, we have

Now let $X-m$ be a ${\rm Bino}(N,p)$ random variable. By taking expectation on both sides, we have:

for $m\geq 1$ and $N\geq 4$. Plug these in (81), we have

where $1/(2Np)\leq 1/(Np+m)$ comes from the fact that $Np\geq m$.

(ii) $Np<m$. In this case, for any $x$, by applying Taylor’s theorem around $x_{0}=Np+m$, there exists $\zeta$ between $x$ and $x_{0}$ such that

By noticing that $\zeta\geq\min\{x,x_{0}\}\geq m\geq(Np+m)/2$, we have:

Similarly, by taking expectation on both sides, we have

Combining the two cases, we obtain the desired statement.

Appendix B Proof of Theorem 2

We use the Efron-Stein inequality to bound the variance of the estimator. For simplicity, let $\widehat{I}^{(N)}(Z)$ be the estimate based on original samples $\{Z_{1},Z_{2},\dots,Z_{N}\}$, where $Z_{i}=(X_{i},Y_{i})$. For the usage of Efron-Stein inequality, we consider another set of i.i.d. samples $\{Z^{\prime}_{1},Z^{\prime}_{2},\dots,Z^{\prime}_{n}\}$ drawn from $P_{XY}$. Let $\widehat{I}^{(N)}(Z^{(j)})$ be the estimate based on $\{Z_{1},\dots,Z_{j-1},Z^{\prime}_{j},Z_{j+1},\dots,Z_{N}\}$. Then Efron-Stein inequality states that

Now we will give an upper bound for the difference $|\widehat{I}^{(N)}(Z)-\widehat{I}^{(N)}(Z^{(j)})|$ for given index $j$. First of all, let $\widehat{I}^{(N)}(Z_{\setminus j})$ be the estimate based on $\{Z_{1},\dots,Z_{j-1},Z_{j+1},\dots,Z_{N}\}$, then by triangle inequality, we have:

where the last equality comes from the fact that $\{Z_{1},\dots,Z_{j-1},Z^{\prime}_{j},Z_{j+1},\dots,Z_{N}\}$ has the same joint distribution as $\{Z_{1},\dots,Z_{N}\}$. Now recall that

Now we need to upper-bound the difference $|\,\xi_{i}(Z)-\xi_{i}(Z_{\setminus j})\,|$ created by eliminating sample $Z_{j}$ for different $i$ ’s. There are three cases of $i$’s as follows,

Case I. $i=j$. Since the upper bounds $|\xi_{i}(Z)|\leq 2\log N$ and $|\xi_{i}(Z_{\setminus j})|\leq 2\log(N-1)$ always holds, so $|\,\xi_{i}(Z)-\xi_{i}(Z_{\setminus j})\,|\leq 4\log N$. The number of $i$’s in this case is only 1. So $\sum_{\textrm{Case I}}|\,\xi_{i}(Z)-\xi_{i}(Z_{\setminus j})\,|\leq 4\log N$.

The number of $i$’s in this case is the number if $i$’s such that $X_{i}=X_{j}$ but $Y_{i}\neq Y_{j}$, which is less than $n_{x,j}$. Therefore, $\sum_{\textrm{Case II.2}}|\,\xi_{i}(Z)-\xi_{i}(Z_{\setminus j})\,|\leq 2n_{x,j}/n_{x,j}\leq 2$.

Combining the four sub-cases, we conclude that $\sum_{\textrm{Case II}}|\,\xi_{i}(Z)-\xi_{i}(Z_{\setminus j})\,|\leq 13$.

Case III.1. $Z_{j}$ is in the $k$-nearest neighbors of $Z_{i}$. In this case, we don’t know how $n_{x,i}$ and $n_{y,i}$ will change by eliminating $Z_{j}$, so we just use the loosest bound $|\,\xi_{i}(Z)-\xi_{i}(Z_{\setminus j})\,|\leq 4\log N$. However, the number of $i$’s in this case is upper bounded by the following lemma.

By the first inequality in Lemma B.1, the number of $i$’s in this case is upper bounded by $k\gamma_{d}$. Therefore, $\sum_{\textrm{Case III.1}}|\,\xi_{i}(Z)-\xi_{i}(Z_{\setminus j})\,|\leq 4k\gamma_{d_{x}+d_{y}}\log N$.

Case III.2. $Z_{j}$ is not in the $k$-nearest neighbors of $Z_{i}$, but $\|X_{j}-X_{i}\|\leq\rho_{i,xy}$, i.e., $X_{j}$ is in the $n_{x,i}$-nearest neighbors of $X_{i}$. In this case, $n_{x,i}$ will decrease by 1 and $n_{y,i}$ remains the same. So

We don’t have an upper bound for the number of $i$’s in this case, but from the second inequality in Lemma B.1, we have the following upper bound, where $\mathcal{X}_{i,j}=\{X_{1},\dots,X_{i-1},X_{j},X_{i+1},\dots,X_{N}\}$:

Case III.3. $Z_{j}$ is not in the $k$-nearest neighbors of $Z_{i}$, but $\|Y_{j}-Y_{i}\|\leq\rho_{i,xy}$, i.e., $Y_{j}$ is in the $n_{y,i}$-nearest neighbors of $Y_{i}$. In this case, $n_{y,i}$ will decrease by 1 and $n_{x,i}$ remains the same. Follow the same analysis in Case III.2, we have $\sum_{\textrm{Case III.2}}|\,\xi_{i}(Z)-\xi_{i}(Z_{\setminus j})\,|\leq 2\gamma_{d_{x}+d_{y}}(\log N+1)$ as well.

Case III.4. $Z_{j}$ is not in the $k$-nearest neighbors of $Z_{i}$, and $\|X_{j}-X_{i}\|>\rho_{i,xy}$, $\|Y_{j}-Y_{i}\|>\rho_{i,xy}$. In this case, neither $n_{x,i}$ nor $n_{y,i}$ will change. Similar to Case II.4, $\sum_{\textrm{Case III.4}}|\,\xi_{i}(Z)-\xi_{i}(Z_{\setminus j})\,|\leq 1$.

Combining the four sub-cases, we conclude that $\sum_{\textrm{Case III}}|\,\xi_{i}(Z)-\xi_{i}(Z_{\setminus j})\,|\leq(4k+4)\gamma_{d_{x}+d_{y}}\log N+4\gamma_{d_{x}+d_{y}}+1$.

for $k\geq 1$, $\log N\geq 1$ and all $\{Z_{1},\dots,Z_{N}\}$. Plug it into (91), we obtain,

Plug it into Efron-Stein inequality (88), we obtain:

Since $1800\gamma^{2}_{d_{x}+d_{y}}$ is a constant independent of $N$, and $(k_{N}\log N)^{2}/N\to 0$ as $N\to\infty$ by Assumption 6, we have $\lim_{N\to\infty}{\textrm{ Var }}\left[\,\widehat{I}^{(N)}(Z)\,\right]=0$.

For the first part of the lemma, we refer to Lemma 20.6 in .

The second part of the lemma is a consequence of the first part. We reorder the indices $i$’s by $k_{i}$ and rewrite the summation as follows,

References