Scalable Mutual Information Estimation using Dependence Graphs
Morteza Noshad Yu Zeng Alfred O. Hero III
cs.IT math.IT stat.ML
 













 
    


            
       
       
       
      












 
 
 
































g
  #1#1



og
  #1
     #2#2 
     #2
    #1
{#2} 





g
  #1
  
  #1



lng
  #1ln ln#1



og
  #1
     #2ℛeℛe{#2} 
     #2ℛe
    #1
ℛe#2 

og
  #1
     #2ℐmℐm{#2} 
     #2ℐm
    #1
ℐm#2 



cosg
  #1coscos#1

sing
  #1sinsin#1

tang
  #1tantan#1

arccosg
  #1arccosarccos#1

arcsing
  #1arcsinarcsin#1

arctang
  #1arctanarctan#1

g
  #1cotcot#1










g
  #1trtr#1






og
  #1
     #2 
     #2 



expg
  #1#1





























g
  #1CovCov#1





















{}og#1#2ℱℱ{#2}#2ℱ
        #1
    ℱ#2og#1#2ℱ^-1ℱ^-1{#2}#2ℱ^-1
        #1
    ℱ^-1#2g#1g#1g#1[1]#1#2#3test@foreach#1#2#3,end@token
    swallow#1test@foreach#1#2ifnextcharend@token
        swallowforeach#1#2foreach#1#2#3,#4end@token
      #1#2#3test@foreach#1#2#4end@token
    theoremTheorem[section]
    lemma[theorem]Lemmaproposition[theorem]Propositioncorollary[theorem]CorollaryremarkRemarkdefinition[1][Definition]
       #1example[1][Example]
       #1ℝł|||̊|ł(()̊)ł[[]̊]𝐗𝐘Ł LMorteza Noshad, Yu Zeng, Alfred O. Hero IIIThis research was partially supported by ARO grant W911NF-15-1-0479.University of Michigan, Electrical Engineering and Computer Science, Ann Arbor, Michigan, U.S.AScalable Mutual Information Estimation using Dependence Graphs
        Christine Nattrass
        December 30, 2023
    ==============================================================
    
    
    The Mutual Information (MI) is an often used measure of dependency between two random variables utilized in information theory, statistics and machine learning. Recently several MI estimators have been proposed that can achieve parametric MSE convergence rate. However, most of the previously proposed estimators have high computational complexity of at least O(N^2).
    We propose a unified method for empirical non-parametric estimation of general MI function between random vectors in ℝ^d based on N i.i.d. samples. The reduced complexity MI estimator, called the ensemble dependency graph estimator (EDGE), combines randomized locality sensitive hashing (LSH), dependency graphs, and ensemble bias-reduction methods. We prove that EDGE achieves optimal computational complexity O(N), and can achieve the optimal parametric MSE rate of O(1/N) if the density is d times differentiable. To the best of our knowledge EDGE is the first non-parametric MI estimator that can achieve parametric MSE rates with linear time complexity. We illustrate the utility of EDGE for the analysis of the information plane (IP) in deep learning. Using EDGE we shed light on a controversy on whether or not the  compression property of information bottleneck (IB) in fact holds for ReLu and other rectification functions in deep neural networks (DNN). 
    
    
    
    § INTRODUCTION
    
    
    The Mutual Information (MI) is an often used measure of dependency between two random variables or vectors <cit.>, and it has a wide range of applications in information theory  <cit.> and machine learning <cit.>. Non-parametric MI estimation methods have been studied that use estimation strategies including KSG <cit.>, KDE <cit.> and Parzen window density estimation <cit.>. The performance of these estimators has been evaluated and compared based on both empirical studies <cit.> and asymptotic analysis <cit.>. Recently several MI estimators have been proposed that can achieve parametric MSE rate of convergence. For example, in <cit.> a KDE plug-in estimator for Rényi divergence and mutual information achieves the MSE rate of O(1/N) when the densities are at least d times differentiable. Another KDE based mutual information estimator was proposed in <cit.> that can achieve the MSE rate of O(1/N) when the densities are d/2 times differentiable. 
    Recently Moon et al <cit.> and Gao et al <cit.> respectively proposed KDE and KNN based MI estimators for random variables with mixtures of continuous and discrete components.
    Most of these estimators, however, have high  computational cost and require knowledge of the density support boundary. 
    
    In this paper we propose a reduced complexity MI estimator called the ensemble dependency graph estimator (EDGE). The estimator combines  randomized locality sensitive hashing (LSH), dependency graphs, and ensemble bias-reduction methods. A dependence graph is a bipartite directed graph consisting of two sets of nodes V and U. The data points are mapped to the sets V and U using a randomized LSH function H that depends on a hash parameter ϵ. Each node is assigned a weight that is proportional to the number of hash collisions. Likewise, each edge between the vertices v_i and u_j has a weight proportional to the number of (X_k,Y_k) pairs mapped to the node pairs (v_i,u_j).
    For a given value of the hash parameter ϵ, a base estimator of MI is proposed as a weighted average of non-linearly transformed of the edge weights. The proposed EDGE estimator of MI is obtained by applying the method of weighted ensemble bias reduction <cit.> to a set of base estimators with different hash parameters.
    This estimator is a non-trivial extension of the LSH divergence estimator defined in <cit.>. LSH-based methods have previously been used for KNN search and graph constructions problems <cit.>, and they result in fast and low complexity algorithms.
    
    Recently, Shwartz-Ziv and Tishby utilized MI to study the training process in Deep Neural Networks (DNN) <cit.>. Let X, T and Y respectively denote the input, hidden and output layers. The authors of  <cit.> introduced the information bottleneck (IB) that represents the tradeoff between two mutual information measures: I(X,T) and I(T,Y). They observed that the training process of a DNN consists of two distinct phases; 1) an initial fitting phase in which I(T,Y) increases, and 2) a subsequent compression phase in which I(X,T) decreases. Saxe  et al in <cit.> countered the claim of <cit.>, asserting that this compression property is not universal, rather it depends on the specific activation function. Specifically, they claimed that the compression property does not hold for ReLu activation functions. The authors of <cit.> challenged these claims, arguing that the authors of <cit.>  had not observed compression due to poor estimates of the MI. We use our proposed rate-optimal ensemble MI estimator to explore this controversy, observing that our estimator of MI does exhibit the compression phenomenon in the ReLU network studied by <cit.>.  
    
    
    Our contributions are as follows:
    
    
      *  
    To the best of our knowledge the proposed MI estimator is the first estimator to have 
     linear complexity and can achieve the optimal  MSE rate of O(1/N).
    
    
      *  
    The proposed MI estimator provides a simplified and unified treatment of mixed continuous-discrete variables. This is due to the hash function approach that is adopted.
    
    
    
    
      *  EDGE is applied to IB theory of deep learning, and provides evidence that the compression property does indeed occur in ReLu DNNs, contrary to the claims of <cit.>. 
    
    
    The rest of the paper is organized as follows. In Section<ref>, we introduce the general definition of MI and define the dependence graph. 
    In Section <ref>, we introduce the hash based MI estimator and give theory for the bias and variance. In section <ref> we introduce the ensemble dependence graph MI estimator (EDGE) and show how the ensemble estimation method can be used to improve the convergence rates.
    Finally, in Section <ref> we provide numerical results as well as study the IP in DNNs.
    
    
    § MUTUAL INFORMATION
    
    
    In this section, we introduce the general mutual information function based on the f-divergence measure. Then, we define a consistent estimator for the mutual information function. 
    Consider the probability measures P and Q on a Euclidean space 𝒳. Let g:(0,∞)→ℝ be a convex function with g(1)=0. The f-divergence between P and Q can be defined as follows <cit.>.
    
        D(PQ):=𝔼_Q[gdP/dQ].
    
    Let 𝒳 and 𝒴 be Euclidean spaces and let P_XY be a probability measure on the space 𝒳×𝒴. For any measurable sets A⊆𝒳 and B⊆𝒴, we define the marginal probability measures P_X(A):=P_XY(A×𝒴) and P_Y(B):=P_XY(𝒳× B). Similar to <cit.>, the general MI denoted by I(X,Y) is defined as
    
        D(P_XYP_XP_Y)=𝔼_P_XP_Y[gdP_XY/dP_XP_Y],
    
    where dP_XY/dP_XP_Y is the Radon-Nikodym derivative, and g:(0,∞)→ℝ is, as in (<ref>) a convex function with g(1)=0. Shannon mutual information is a particular cases of (<ref>) for which g(x)=xlog x. 
    
    
     §.§ Dependence Graphs
    
    
    
    
    
    Consider N i.i.d samples (X_i,Y_i), 1≤ i≤ N drawn from the probability measure P_XY, defined on the space 𝒳×𝒴. 
    Define the sets 𝐗={X_1,X_2,...,X_N} and 𝐘={Y_1,Y_2,...,Y_N}.
    The dependence graph G(X,Y) is a directed bipartite graph, consisting of two sets of nodes V and U with cardinalities denoted as |V| and |U|, and the set of edges E_G. Each point in the sets 𝐗 and 𝐘 is mapped to the nodes in the sets U and V, respectively, using the hash function H, described as follows.
    
    
    
    
    A vector valued hash function H is defined in a similar way as defined in <cit.>. First, define the vector valued hash function H_1: ℝ^d→ℤ^d as 
    
        H_1(x)=[h_1(x_1), h_1(x_2), ... ,h_1(x_d)],
    
    where x_i denotes the ith component of the vector x. In (<ref>), each scalar hash function h_1(x_i): ℝ→ℤ is given by
    
        h_1(x_i)=x_i+b/ϵ,
    
    for a fixed ϵ>0, where ⌊ y⌋ denotes the floor function (the smallest integer value less than or equal to y), and b is a fixed random variable in [0,ϵ]. Let ℱ:={1,2,..,F}, where F:=c_HN and c_H is a fixed tunable integer. We define a random hash function H_2:ℤ^d→ℱ with a uniform density on the output and consider the combined hashing function
    
        H(x):=H_2(H_1(x)),
    
    which maps the points in ℝ^d to ℱ.
    
    H(x) reveals the index of the mapped vertex in G(X,Y). The weights ω_i and ω'_j corresponding to the nodes v_i and u_j, and ω_ij, the weight of the edge (v_i,u_j), are defined as follows.
    
        ω_i=N_i/N,ω'_j=M_j/N,ω_ij=N_ijN/N_iM_j,
    
    where N_i and M_j respectively are the the number of hash collisions at the vertices v_i and u_j, and N_ij is the number of joint collisions of the nodes (X_k,Y_k) at the vertex pairs (v_i,u_j). The number of hash collisions is defined as the number of instances of the input variables map to the same output value. In particular, 
    
        N_ij:=#{(X_k,Y_k)  s.t   H(X_k)=i  and   H(Y_k)=j }.
    
    Fig. <ref> represents a sample dependence graph. Note that the nodes and edges with zero collisions do not show up in the dependence graph. 
    
    
    
    
    § THE BASE ESTIMATOR OF MI
    
    
    
     §.§ Assumptions
    
    
    The following are the assumptions we make on the probability measures and g:
    
    A1.  The support sets 𝒳 and 𝒴 are bounded.
    
    A2.  The following supremum exists and is bounded:
    
        sup_P_XP_YgdP_XY/dP_XP_Y≤ U.
    A3. Let x_D and x_C respectively denote the discrete and continuous components of the vector x. Also let f_X_C(x_C) and p_X_D(x_D) respectively denote density and pmf functions of these components associated with the  probability measure P_X. The density functions f_X_C(x_C), f_Y_C(y_C), f_X_CY_C(x_C,y_C), and the conditional densities f_X_C|X_D(x_C|x_D), f_Y_C|Y_D(y_C|y_D), f_X_CY_C|X_DY_D(x_C,y_C|x_D,y_D) are Hölder continuous.
    
    
    
    Given a support set 𝒳, a function f:𝒳→ℝ is called Hölder continuous with parameter 0<γ≤ 1, if there exists a positive constant G_f, possibly depending on f, such that for every x≠ y ∈𝒳,
    
        |f(y)-f(x)|≤ G_fy-x^γ.
    A4.  Assume that the function g in (<ref>) is Lipschitz continuous; i.e. g is Hölder continuous with γ=1.
    
    
    
     §.§ The Base Estimator of MI
    
    For a fixed value of the hash parameter ϵ, we propose the following base estimator of MI (<ref>) function based on the dependence graph:
    
        I(X,Y):=∑_e_ij∈ E_Gω_iω'_jgω_ij,
    
    where the summation is over all edges e_ij:(v_i→ u_j) of G(X,Y) having non-zero weight and g(x):=max{g(x),U}. 
    
    When X and Y are strongly dependent, each point X_k hashed into the bucket (vertex) v_i corresponds to a unique hash value for Y_k in U. Therefore, asymptotically ω_ij→ 1 and the mutual information estimation in (<ref>) takes its maximum value. On the other hand, when X and Y are independent, each point X_k hashed into the bucket (vertex) v_i may be associated with different values of Y_k, and therefore asymptotically ω_ij→ω_j and the Shannon MI estimation tends to 0.
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
     §.§ Various LSH Functions
    
    There are various types of LSH functions <cit.>, and all of them share the common property that they map similar items to the same bins with high probability.
    
    In equations (<ref>) and (<ref>) we considered a simple floor function on the scaled input, however in general, any other type of LSH might be used for our estimation method. In particular, the hash functions based on random projections can reduce the dimensionality of data. SimHash <cit.>, which is based on cosine distance, and the LSH based on p-stable distributions <cit.> are among well known LSH functions that reduce the dimension of data. For example, the LSH based on p-stable
    distribution is defined similarly to the floor hash function in (<ref>) and (<ref>), except that the input vector is projected on random hyperplanes with p-stable distributions. The formal definition is H_p-stable:ℝ^d→ℤ^r,
    
    
        H_p-stable(x)= H_1(XW),
    
    where H_1 is defined in (<ref>), and W is a d× r matrix with entries chosen independently from a stable distribution. For high-dimensional datasets one can choose r<<d in order to reduce the  dimensionality. Finally, note that for theoretical analysis, we only focus on performance of the simple floor hash function defined in (<ref>) and (<ref>). 
    
    
    
    
    
    
    
    
     §.§ Convergence Rates
    
    In the following theorems we state upper bounds on the bias and variance rates of the proposed MI estimator (<ref>). The proofs are given in appendices A and B.
    We define the notations 𝔹[T̂]=𝔼[T̂]-T for bias and 𝕍[T̂]=𝔼[T̂^2]-𝔼[T̂]^2 for variance of T̂. The following theorem states an upper bound on the bias.
    
    
    
    Let d=d_X+d_Y be the dimension of the joint random variable (X,Y). Under the aforementioned assumptions A1-A4, and assuming that the density functions in A3 have bounded derivatives up to order q≥ 0, the following upper bound on the bias of the estimator in (<ref>) holds
    
        𝔹[I(X,Y)] = 
             
               Oϵ^γ+O1/Nϵ^d,       q=0
        
        ∑_i=1^q C_iϵ^i+Oϵ^q+O1/Nϵ^d      q≥ 1,
    
    where  ϵ is the hash parameter in (<ref>), γ is the smoothness parameter in (<ref>), and C_i are real constants.
    
    
    
    In (<ref>), the hash parameter, ϵ needs to be a function of N to ensure that the bias converges to zero. For the case of q=0, the optimum bias is achieved when ϵ=1/N^γ/(γ+d). When q≥ 1, the optimum bias is achieved for ϵ=1/N^1/(1+d).
    
    
    
    Under the assumptions A1-A4 the variance of the proposed estimator can be bounded as 𝕍[I(X,Y)]≤ O1/N.
    Further, the variance of the variable ω_ij is also upper bounded by O(1/N).
    
    
     §.§ Computational Complexity
    
    
    We analyze the computational complexity of the proposed estimator. The procedure for estimation of MI equivalent to the proposed estimator in (<ref>) is given in Algorithm <ref>.
    
    
    
    
    
    We go over all of the data points and map them using the hash function H. We compute the number of marginal and joint collisions (N_i,M'_j, N_ij) and based on these we compute the vertex and edge weights. Note that computing the hashing of all of the data points takes about O(N) time. Finally in the last line we compute the MI estimate by going over all of the edges e_ij. The number the edges is upper bounded by O(N), since each edge correspond to at least one pair of (X_i,Y_i). Finally, note that for high-dimensional data sets, the computational of computing the hash function of each input may depend on d, however, is would not be greater than O(d). Hence, overall the the computational complexity of computing the proposed MI estimate is linear with respect to both N and d.
    
    
    
    
    
     §.§ Comparison to the Other Estimation Methods
    
    
    So far, various estimation methods for information measures (entropy, divergence and mutual information) have been proposed, most of which are based on kernel density estimates (KDE) <cit.>, k-nearest neighbors (KNN) <cit.> or histogram binning <cit.>. Certain KDE and KNN based estimators can achieve the optimal parametric MSE rate <cit.>, however, implementation of the KDE and KNN methods respectively have O(N^2) and O(kNlog N) computational complexity, where N is the number of samples. One could probably could approximation methods for KDE and KNN, however, there would be no theoretical guarantees for the estimation based on these approximations. Empirical histograms, on the other hand, are simpler and easier to implement, however, their convergence rate is not as good as the KNN and KDE based methods <cit.>. 
    
    LSH methods have previously been applied to the approximate nearest neighbor search, however, it has never been directly utilized for estimation of densities or information theoretic quantities.  
    The simplest LSH function considered in (<ref>) has similarities with the histogram estimator in terms of binning, however, there are also certain differences. As opposed to the LSH based method, the histogram binning requires a pre-knowledge of the support set or needs an extra computation to estimate the support set. The number of 
    the bins in the histogram estimator gets exponentially large with increasing dimension which results in a huge computational complexity for high-dimensional datasets. In addition, most of the bins would be empty. In contrast, the LSH based method results in no empty hash bins, the number of the bins is upper bounded by O(N), and the computational complexity is linear in dimension and the number of samples. The plug-in methods including histograms, KDE and KNN require the estimation of the densities p_X, p_Y and p_XY for computation of mutual information, while the proposed LSH-based method finds the mapping of  and  data points, and then estimate the mutual information, based on the hash collisions. Finally, by giving an accurate bias and variance rates for the base LSH estimator, we use an ensemble estimation technique to achieve the optimal parametric convergence rate, discussed in the following section.
    
    
    
    
    
    
    
    § ENSEMBLE DEPENDENCE GRAPH ESTIMATOR (EDGE)
    
    
    Given the expression for the bias in Theorem <ref>, the ensemble estimation technique proposed in <cit.> can be applied to improve the convergence rate of the MI estimator (<ref>). Assume that the densities in A3 have continuous bounded derivatives up to the order q, where q≥ d.
    Let 𝒯:={t_1,...,t_T} be a set of index values with t_i<c, where c>0 is a constant. Let ϵ(t):=tN^-1/2d. For a given set of weights w(t) the weighted ensemble estimator is then defined as 
    
        I_w:=∑_t∈𝒯w(t)I_ϵ(t),
    
    where I_ϵ(t) is the mutual information estimator with the parameter ϵ(t). Using (<ref>), for q>0 the bias of the weighted ensemble estimator (<ref>) takes the form
    
        𝔹(Î_w) = ∑_i=1^q Ci N^-i/2d∑_t∈𝒯 w(t) t^i +Ot^d/N^1/2+O1/Nϵ^d
    
    
    Given the form (<ref>), as long as T≥ q, we can select the  weights w(t) to force to zero the slowly decaying terms in (<ref>), i.e.  ∑_t∈τ w(t)t^i/d=0 subject to the constraint that∑_t∈τ w(t)=1. However, T should be strictly greater than q in order to control the variance, which is upper bounded by the euclidean norm squared of the weights  ω. In particular we have the following theorem (the proof is given in Appendix C):
    
    
    For T>d let  w_0 be the solution to:
    
        min_w        w_2 
        subject to       ∑_t∈𝒯w(t)=1, 
               ∑_t∈𝒯w(t)t^i=0, i∈ℕ, i≤ d.
    
    Then the MSE rate of the ensemble estimator I_w_0 is O(1/N).
    
    
    § EXPERIMENTS
    
    
    We first use simulated data to compare the proposed estimator to the competing MI estimators Ensemble KDE (EKDE)<cit.>, and generalized KSG <cit.>. 
    Both of these estimators work on mixed continuous-discrete variables. We also apply EDGE to study the information bottleneck in different networks trained on MNIST hand-written datset.
    
    
    Fig. <ref>, shows the MSE estimation rate of Shannon MI between the continuous random variables X and Y having the  relation Y=X+a N_U,
    where X is a 2D Gaussian random variable with the mean [0, 0] and covariance matrix C=I_2. Here I_d denote the d-dimensional identity matrix. N_U is a uniform random vector with the support 𝒩_U=[0,1]× [0,1].  
    We compute the MSE of each estimator for different sample sizes. The MSE rates of EDGE, EKDE and KSG are compared for a=1/5. Further, the MSE rate of EDGE is investigated for noise levels of a={1/10,1/5,1/2,1}. As the dependency between X and Y increases the MSE rate becomes slower. 
    
    
    
    
    
    
    
    Fig. <ref>, shows the MSE estimation rate of Shannon MI between a discrete random variables X and a continuous random variable Y.
    We have X∈{1,2,3,4}, and each X=x is associated with multivariate Gaussian random vector Y, with d=4, the expectation [x/2, 0, 0, 0] and covariance matrix C=I_4. 
    In general in Figures <ref> and <ref>, EDGE has better convergence rate than EKDE and KSG estimators. Fig. <ref> represents the runtime comparison for the same experiment as in Fig. <ref>. It can be seen from this graph how fast our proposed estimator performs compared to the other other methods.
    
    
    
    
    
    
    
    
    Next, we use EDGE to study the information bottleneck <cit.> in DNNs. Fig. <ref> represents the information plane of a DNN with 4 fully connected hidden layers of width 784-1024-20-20-20-10 with tanh and ReLU activations. The sequence of colored points shows different iterations of the training process. Each gray line connects the points with the same iterations for diferent layers. The left most sequence of points corresponds to the last hidden layer and the right most sequence of points corresponds to the first hidden layer.
    The network is trained with Adam optimization with a learning rate of 0.003 and cross-entropy loss functions to classify the MNIST handwritten-digits dataset. 
    We repeat the experiment for 20 iterations with different randomized initializations and take the average over all experiments. In both cases of ReLU and tanh activations we observe some degree of compression in all of the hidden layers. However, the amount of compressions is different for ReLU and tanh activations. The average test accuracy in both of these networks are around 0.98. This network is the same as the one studied in <cit.>, for which it is claimed that no compression happens with a ReLU activation. The base estimator used in <cit.> provides KDE-based lower and upper bounds on the true MI <cit.>. According to our experiments (not shown) the upper bound is in some cases twice as large as the lower bound. In contrast, our proposed ensemble method estimates the exact mutual information with significantly higher accuracy. 
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    Fig. <ref> represents the information plane for another network with 4 fully connected hidden layers of width 784-200-100-60-30-10 with ReLU activation. The network is trained with Adam optimization with a learning rate of 0.003 and cross-entropy loss functions to classify the MNIST handwritten-digits dataset. Again, we observe compression for this network with ReLU activation.
    
    
    
    
    
    
    
    
    
    
    
    
    
    Finally, we study the information plane curves in a CNN with three convolutioal ReLU layers and a dense ReLU layer. The convolutional layers respectively have depths of 4,8,16 and  the dense layer has the dimension 256. 
    Max-pooling functions are used in the second and third layers. Note although for a certain initialization of the weights this model can achieve the test accuracy of 0.99, the average test accuracy (over different weight initializations) is around 0.95. That's why the converged point of the last layer has smaller I(T,Y) compared to the examples in Fig. <ref>, which achieves the average test accuracy of 0.98. Another interesting point about the information plane in CNN is that the convolutional layers have larger I(T,Y) compared to the hidden layers in the fully connected models in <ref> and <ref>, which implies that the convolutional layers can extract almost all of the useful information about the labels after small number of iterations.
    
    
    
    
    
    
    § CONCLUSION
    
    
    In this paper we proposed a fast non-parametric estimation method for MI based on random hashing, dependence graphs, and ensemble estimation. Remarkably, the proposed estimator has linear computational complexity and attains optimal (parametric) rates of MSE convergence. We provided bias and variance convergence rate, and validated our results by numerical experiments.
    
    
    
    
    ieeetr
    
    § A. BIAS PROOF
    
    We first prove a theorem that establishes an upper bound on the number of vertices in V and U.
    
    
    Cardinality of the sets U and V are upper bounded as |V|≤ Oϵ^-d and |U|≤ Oϵ^-d, respectively.
    
    
    Let {X_i}_i=1^L_X and {Y_i}_i=1^L_Y
     respectively denote distinct outputs of H_1 with the N i.i.d points X_k and Y_k as input. Then according to <cit.> (Lemma 4.1), we have
    
        L_X≤ Oϵ^-d,  L_Y≤ Oϵ^-d.
    
    
    Simply, because of the deterministic feature of H_2, the number of its distinct inputs is greater than or equal to the number of its outputs. So, |V|≤ L_X and |U|≤ L_Y. Using the bounds in (<ref>) completes the proof.  
    
    
    The bias proof is based on analyzing the hash function defined in (<ref>). The proof consists of two main steps: 1) Finding the expectation of hash collisions of H_1; and 2) Analyzing the collision error of H_2. An important point about H_1 and H_2 is that collision of H_1 plays a crucial role in our estimator, while the collision of H_2 adds extra bias to the estimator. We introduce the following events to formally define these two biases:
    
    
        E_ij:    The event that there is an edge between the vertices v_i and u_j.
         
        E_ℰ:    The event that ℰ is the set of all edges in G, i.e. ℰ=E_G.
        
        E_v_i^>0:    The event that there is at least one vector from {X_i}_i=1^L_X that maps to v_i using H_2.
        
        E_v_i^=1:    The event that there is exactly one vector from {X_i}_i=1^L_X that maps to v_i using H_2.
        
        E_v_i^>1:    The event that there are at least two vectors from {X_i}_i=1^L_X that map to v_i using H_2.
    E_u_i^>0, E_u_i^=1 and E_u_i^>1 are defined similarly. Further, let for any event E, E denote the complementary event. Let E_ij^=1:= E_v_i^=1∩ E_u_i^=1. Finally, we define E^=1:=∩_i=1^L_XE_v_i^=1∩∩_j=1^L_YE_u_j^=1, which represent the event of no collision.
    
    
    Consider the notation I(X,Y):=∑_e_ij∈ E_Gω_iω'_jgω_ij (Notice the difference from the definition in (<ref>)). We can derive its expectation as
    
    
        I(X,Y)   =∑_e_ij∈ E_Gω_iω'_jgω_ij| E_G
           =∑_e_ij∈ E_Gω_iω'_jgω_ij| E_ij
           =∑_e_ij∈ E_GP(E_ij^=1| E_ij) ω_iω'_jgω_ij| E_ij^=1,E_ij
           +∑_e_ij∈ E_GPE_ij^=1| E_ijω_iω'_jgω_ij|E_ij^=1,E_ij.
    
    Note that the second term in (<ref>) is the bias due to collision of H_2 and we denote this term by 𝔹_H. 
    
    
    
    
     §.§ Bias Due to Collision
    
    
    The following lemma states an upper bound on the bias error caused by H_2.
    
    The bias error due to collision of H_2 is upper bounded as
    
        𝔹_H ≤ O1/ϵ^dN.
    
    
    
    Before proving this lemma, we provide the following lemma.
    
    P(E_ij^=1| E_ij) is given by
    
        P(E_ij^=1| E_ij)= 1-O1/ϵ^dN.
    
    
    Let 𝐗=𝐱 and 𝐘=𝐲 respectively abbreviate the equations 
    X_1=x_1,...,X_L_X=x_L_X and Y_1=y_1,...,Y_L_Y=y_L_Y. Let 𝐱:={x_1,x_2,...,x_L_X} and 𝐲:={y_1,y_2,...,y_L_Y}. Define 𝐳:=𝐱∪𝐲 and L_Z:=|𝐳|.
    
    
    
        P(E_ij^=1| E_ij)   =∑_𝐱,𝐲 P𝐗=𝐱,𝐘=𝐲| E_ijP(E_ij^=1| E_ij, 𝐗=𝐱,𝐘=𝐲).
    
    
    Define a=2 for the case i≠ j and a=1 for the case i=j. Then we have
    
    
        P(E_ij^=1| E_ij)   =∑_𝐱,𝐲 P𝐗=𝐱,𝐘=𝐲| E_ijOF-a/F^L_Z-a
           =∑_𝐱,𝐲 P𝐗=𝐱,𝐘=𝐲| E_ij1-OL_Z/F
           ≤∑_𝐱,𝐲 P𝐗=𝐱,𝐘=𝐲| E_ij1-OL_X+L_Y/F
           =∑_𝐱,𝐲 P𝐗=𝐱,𝐘=𝐲| E_ij1-O1/ϵ^dN
           =1-O1/ϵ^dN∑_𝐱,𝐲 P𝐗=𝐱,𝐘=𝐲| E_ij
           =1-O1/ϵ^dN,
    
    where in the fourth line we have used (<ref>).
    N'_i and M'_j respectively are defined as the number of the input points 𝐗 and 𝐘 mapped to the buckets X_i and Y_j using H_1. Define 𝒜_i:= {j: H_2(X_j)=i} and ℬ_i:= {j: H_2(Y_j)=i}. For each i we can rewrite N_i and M_i as
    
        N_i=∑_j=1^L_X1_𝒜_i(j)N_j',  M_i=∑_j=1^L_Y1_ℬ_i(j)M_j'.
     
    
    
    Thus, 
    
    
        𝔹_H    ≤∑_i,j∈ℱ PE_ij^> 11_E_ijω_iω'_jgω_ij| E_ij^> 1
           = ∑_i,j∈ℱ PE_ij^> 1PE_ij| E_ij^> 1ω_iω'_jgω_ij| E_ij^> 1,E_ij+
        PE_ij| E_ij^> 1ω_iω'_jgω_ij| E_ij^> 1,E_ij
           = ∑_i,j∈ℱ PE_ijPE_ij^> 1| E_ijω_iω'_jgω_ij| E_ij^> 1,E_ij
           ≤ OU /ϵ^dN∑_i,j∈ℱ PE_ijω_iω'_j| E_ij^>1,E_ij
           = OU /ϵ^dN^3∑_i,j∈ℱ PE_ijN_iM_j| E_ij^>1,E_ij
           = OU /ϵ^dN^3∑_𝐱,𝐲 p_𝐗,𝐘𝐱,𝐲∑_i,j∈ℱ PE_ijN_iM_j| E_ij^>1,E_ij, 𝐗=𝐱,𝐘=𝐲
           = OU /ϵ^dN^3∑_𝐱,𝐲p_𝐗,𝐘∑_i,j∈ℱ PE_ij∑_r=1^L_X1_𝒜_i(r)N_r'∑_s=1^L_Y1_ℬ_j(s)M_s'| E_ij^>1,E_ij, 𝐗=𝐱,𝐘=𝐲
           = OU /ϵ^dN^3∑_𝐱,𝐲 p_𝐗,𝐘∑_i,j∈ℱ PE_ij∑_r=1^L_X∑_s=1^L_Y1_𝒜_i(r)1_ℬ_j(s)| E_ij^>1,E_ij, 𝐗=𝐱,𝐘=𝐲N_r'M_s'| E_ij, 𝐗=𝐱,𝐘=𝐲
           = OU /ϵ^dN^3∑_𝐱,𝐲 p_𝐗,𝐘∑_i,j∈ℱ PE_ij∑_r=1^L_X∑_s=1^L_YPr∈𝒜_i, s∈ℬ_j| E_ij^>1,E_ij, 𝐗=𝐱,𝐘=𝐲N_r'M_s'| E_ij, 𝐗=𝐱,𝐘=𝐲,
    
    where in (<ref>) we have used the Bayes rule, and the fact that gω_ij=0 conditioned on the event E_ij. In (<ref>) we have used the bound in Lemma <ref>, and the upper bound on g(ω_ij) . Equation (<ref>) is due to (<ref>).
    Now we simplify Pr∈𝒜_i, s∈ℬ_j| E_ij^>1,E_ij, 𝐗=𝐱,𝐘=𝐲 in (<ref>) as follows. First assume that X̃_r≠Ỹ_s.
    
    
        Pr∈𝒜_i, s∈ℬ_j| E_ij^>1,E_ij, 𝐗=𝐱,𝐘=𝐲   ≤ Pr∈𝒜_i, s∈ℬ_j| E_v_i^>1, E_u_j^>1, 𝐗=𝐱,𝐘=𝐲
           = Pr∈𝒜_i| E_v_i^>1, 𝐗=𝐱 Ps∈ℬ_j| E_u_j^>1,𝐘=𝐲,
    
    where the second line is because the hash function H_2 is random and independent for different inputs.
    Pr∈𝒜_i| E_v_i^>1,𝐗=𝐱 in (<ref>) can be written as
    
        Pr∈𝒜_i| E_v_i^>1,𝐗=𝐱   =Pr∈𝒜_i,E_v_i^>1|𝐗=𝐱/PE_v_i^>1|𝐗=𝐱.
    
    We first find PE_v_i^>1|𝐗=𝐱:
    
    
        PE_v_i^>1|𝐗=
        𝐱   =1-PE_v_i^=0|𝐗=
        𝐱-PE_v_i^=1|𝐗=
        𝐱
           = 1-F-1/F^L_X-L_X/FF-1/F^L_x-1
           =L_X^2/2F^2+oL_X^2/2F^2.
    
    Next, we find Pr∈𝒜_i,E_v_i^>1|𝐗=𝐱 in (<ref>) as follows.
    
        Pr∈𝒜_i,E_v_i^>1|𝐗=𝐱   =PE_v_i^>1| r∈𝒜_i,𝐗=𝐱
        Pr∈𝒜_i|𝐗=𝐱
           =1-F-1/F^L_X-11/F=OL_X/F^2
    
    Thus, using (<ref>) and (<ref>) yields 
    
    
        Pr∈𝒜_i| E_v_i^>1,𝐗=𝐱
        =O1/L_X.
    
    Similarly, we have  
    
        Ps∈ℬ_j| E_u_j^>1,𝐘=𝐲
        =O1/L_Y.
    
    
    Now assume the case X̃_r=Ỹ_s. Then since H_2(X̃_r)=H_2(Ỹ_s), we can simplify Pr∈𝒜_i, s∈ℬ_j| E_ij^>1,E_ij, 𝐗=𝐱,𝐘=𝐲 in (<ref>) as
    
    
        Pr∈𝒜_i, s∈ℬ_j| E_ij^>1,E_ij, 𝐗=𝐱,𝐘=𝐲=δ_ij Pr∈𝒜_i| E_v_i^>1, 𝐗=𝐱,𝐘=𝐲.
    
    Recalling the definition 𝐳:=𝐱∪𝐲 and L_Z:=|𝐳|, similar to  
    
    
        Pr∈𝒜_i| E_v_i^>1, 𝐗=𝐱,𝐘=𝐲
        =O1/L_Z.
    
    
    By using equations (<ref>), (<ref>), (<ref>) and (<ref>) in (<ref>), we can write the following upper bound for the bias estimator due to collision. 
    
    
        𝔹_H    ≤ OU /ϵ^dN^3∑_𝐱,𝐲 p_𝐗,𝐘∑_i,j∈ℱ PE_ij∑_r=1^L_X∑_s=1^L_YN_r'M_s'| E_ij, 𝐗=𝐱,𝐘=𝐲O1/L_XL_Y+δ_ijO1/L_Z
           = OU /ϵ^dN^3∑_𝐱,𝐲 p_𝐗,𝐘∑_i,j∈ℱ PE_ij∑_r=1^L_XN_r'∑_s=1^L_YM_s'| E_ij, 𝐗=𝐱,𝐘=𝐲O1/L_XL_Y+δ_ijO1/L_Z
           = OU /ϵ^dN^3∑_𝐱,𝐲 p_𝐗,𝐘∑_i,j∈ℱ PE_ij N^2O1/L_XL_Y+δ_ijO1/L_Z
           = OU /ϵ^dN^3∑_𝐱,𝐲 p_𝐗,𝐘ON^2/L_XL_Y+ON/L_Z∑_i,j∈ℱ PE_ij
           = OU /ϵ^dN^3∑_𝐱,𝐲 p_𝐗,𝐘ON^2/L_XL_Y+ON/L_Z∑_i,j∈ℱ1_E_ij
           ≤ OU /ϵ^dN^3∑_𝐱,𝐲 p_𝐗,𝐘ON^2/L_XL_Y+ON/L_ZL_XL_Y
           ≤ O1/ϵ^dN.
    
    
     §.§ Bias without Collision
    
    
    A key idea in proving Theorem <ref> is to show that the expectation of the edge weights ω_ij are proportional to the Radon-Nikodym derivative dP_XY/dP_XP_Y at the points that correspond to the vertices v_i and u_j. This fact is stated in the following lemma:
    
    
    Under the assumptions A1-A4, and assuming that the density functions in A3 have bounded derivatives up to order q≥ 0 we have:
    
        ω_ij=dP_XY/dP_XP_Y+𝔹(N,ϵ,q,γ),
    
    where
    
        𝔹(N,ϵ,q,γ):= 
             
               Oϵ^γ+O1/Nϵ^d,       q=0
        
        ∑_i=1^q C_iϵ^i+Oϵ^q+O1/Nϵ^d,       q≥ 1,
    
    and C_i are real constants.
    
    
    
    
    Note that since ω_ij=N_ijN/N_iM_j, and N_ij, N_i and N_j are not independent variables, deriving the expectation is not trivial.
     In the following we give a lemma that provides conditions under which the expectation of a function of random variables is close to the function of expectations of the random variables. We will use the following lemma to simplify ω_ij.
    
    
    Assume that g(Z_1,Z_2,...,Z_k): 𝒵_1×...×𝒵_k → R is a Lipschitz continuous function with constant H_g> 0 with respect to each of variables Z_i, 1≤ i≤ k. Let   𝕍[Z_i] and 𝕍[Z_i|X] respectively denote the variance and the conditional variance of each variable Z_i for a given variable X. Then we have
    
    
        a) gZ_1,Z_2,...,Z_k-gZ_1,Z_2,...,Z_k≤ H_g∑_i=1^k√(𝕍[Z_i]),
           b) gZ_1,Z_2,...,Z_k|X-gZ_1|X,Z_2|X,...,Z_k|X≤ H_g∑_i=1^k√(𝕍[Z_i|X]).
    
        gZ_1,Z_2,...,Z_k-gZ_1,Z_2,...,Z_k   = gZ_1,Z_2,...,Z_k-gZ_1,Z_2,...,Z_k
           ≤gZ_1,Z_2,...,Z_k-
        gZ_1,Z_2,...,Z_k
           ≤𝔼[| gZ_1,Z_2,...,Z_k-gZ_1,Z_2,...,Z_k+
             +gZ_1,Z_2,...,Z_k-gZ_1,Z_2,...,Z_k
             +...
             +gZ_1,Z_2,...,Z_k-1,Z_k-gZ_1,Z_2,...,Z_k|]
           ≤| gZ_1,Z_2,...,Z_k-gZ_1,Z_2,...,Z_k|
             +| gZ_1,Z_2,...,Z_k-gZ_1,Z_2,...,Z_k|
             +...
             +| gZ_1,...,Z_k-1,Z_k-gZ_1,...,Z_k|
           ≤ H_gZ_1-Z_1+H_gZ_2-Z_2+...+H_gZ_k-Z_k
           ≤ H_g∑_i=1^k√(𝕍[Z_i]).
    
    In (<ref>) and (<ref>) we have used triangle inequalities. In (<ref>) we have applied Lipschitz condition, and finally in (<ref>) we have used Cauchy Schwarz inequality. Since the proofs of parts (a) and (b) are similar, we omit the proof of part (b).
    
    Define ν_ij=N_ij/N, and recall the definitions ω_ij=N_ijN/N_iN_j, ω_i=N_i/N, and ω_j'=N_j/N. Then we can write
    
        ω_ij=ν_ij/ω_iω_j'+O√(1/N)
    
    The proof follows by Lemma <ref> and the fact that 𝕍[ω_ij]≤ O1/N (proved in Lemma <ref>).
    
    
    
    Let x_D and x_C respectively denote the discrete and continuous components of the vector x, with dimensions d_D and d_C. Also let f_X_C(x_C) and p_X_D(x_D) respectively denote density and pmf functions of these components associated with the  probability measure P_X. Let S(x,r) be the set of all points that are within the distance r/2 of x in each dimension i, i.e. 
    
        S(x,r):{x|∀ i≤ d,|X_i-x_i|<r/2}.
    
    
    Denote P_r(x):=P(x∈ S(x,r)). Then we have the following lemma.
    
    
    Let r<s_𝒳, where s_𝒳 is the smallest possible distance in the discrete components of the support set, 𝒳.  Under the assumption A3, and assuming that the density functions in A3 have bounded derivatives up to the order q≥ 0, we have
    
        P_r(x)=P(X_D=x_D)r^d_Cf(x_C|x_D)+μ(r,γ,q,𝐂_X),
    
    where 
    
        μ(r,γ,q,𝐂_X):= 
             
               Or^γ,       q=0
        
        ∑_i=1^q C_i r^i+Or^q,       q≥ 1.
    
    
    In the above equation, 𝐂_X:=(C_1,C_2,...,C_q), and C_i are real constants depending on the probability measure P_X.
    
    
    The proof is straightforward by using (<ref>) for the case q=0 (similar to (27)-(29) in <cit.>), and using the Taylor expansion of f(x_C|x_D) for the case q≥ 1 (similar to (36)-(37) in <cit.>). 
    
    Let H(x)=i, H(y)=j. Under the assumptions A1-A3, and assuming that the density functions in A3 have bounded derivatives up to the order q≥ 0, we have
    
    
        ω_ij|E_ij^≤ 1=dP_XY/dP_XP_Y(x,y)+μ(ϵ,γ,q,𝐂_XY')+O1/√(N),
    
    where μ(ϵ,γ,q,𝐂'_XY) is defined in (<ref>).
    
    Define ν_ij=N_ij/N, and recall the definitions ω_ij=N_ijN/N_iN_j, ω_i=N_i/N, and ω_j'=N_j/N. Using Lemma <ref> we have
    
    
        ω_ij|E_ij^≤ 1=ν_ij|E_ij^≤ 1/ω_i|E_ij^≤ 1ω_j'|E_ij^≤ 1+O1/√(N)
    
    
    
    
    Assume that H(x)=i. Let 𝒳 have d_C and d_D continuous and discrete components, respectively. Also let 𝒴 have d'_C and d'_D continuous and discrete components, respectively. Then we can write
    
    
        ω_i|E_ij^≤ 1   =1/NN_i|E_ij^≤ 1
           =P(X∈ S(x,ϵ))
           =P(X_D=x_D)ϵ^d_Cf(x_C|x_D)+μ(ϵ,γ,q,𝐂_X),
    
    where in the third line we have used Lemma <ref>. Similarly we can write 
    
        ω'_j|E_ij^≤ 1   =P(Y_D=y_D)ϵ^d'_Cf(y_C|y_D)+μ(ϵ,γ,q,𝐂_X), 
        ν_ij|E_ij^≤ 1   =P(X_D=x_D,Y_D=y_D)ϵ^(d_C+d'_C)f(x_C,y_C|x_D,y_D)+μ(ϵ,γ,q,𝐂_XY).
    
    
    Using (<ref>) and (<ref>) in (<ref>) results in
    
    
        ω_ij|E_ij^≤ 1=P(X_D=x_D)P(Y_D=y_D)f(x_C|x_D)f(y_C|y_D)/P(X_D=x_D,Y_D=y_D)f(x_C,y_C|x_D,y_D)+μ(ϵ,γ,q,𝐂'_XY)+O1/√(N),
    
    where 𝐂'_XY depends only on P_XY. 
    Now note that using Lemma <ref>, dP_XY/dP_XP_Y(x,y) can be simplified as 
    
    
        dP_XY/dP_XP_Y(x,y)=dP_XY,r/dr(x,y)/dP_X,rP_Y,r/dr(x,y)=P(X_D=x_D)P(Y_D=y_D)f(x_C|x_D)f(y_C|y_D)/P(X_D=x_D,Y_D=y_D)f(x_C,y_C|x_D,y_D)+μ(ϵ,γ,q,𝐂”_XY).
    
    Finally, using (<ref>) in (<ref>) gives
    
    
        ω_ij|E_ij^≤ 1=dP_XY/dP_XP_Y(x,y)+μ(ϵ,γ,q,𝐂_XY)+O1/√(N),
    
    where H(x)=i,H(y)=j.
    
    Lemma <ref> is a simple consequence of Lemma <ref>. We have
    
    
        ω_ij   = PE_ij^≤ 1ω_ij|E_ij^≤ 1+PE_ij^> 1ω_ij|E_ij^> 1.
    
    
    Recall the definitions 𝐗:=X_1,X_2,...,X_L_X  and 𝐘:=Y_1,Y_2,...,Y_L_Y as the mapped 𝐗 and 𝐘 points through H_1. Let 𝐙:=𝐗∪𝐘 and L_Z:=|𝐙|. We first find PE_ij^≤ 1 as follows. For a fixed set 𝐙 we have
    
    
        PE_ij^≤ 1   =PE_v_i^=0∩ E_u_j^=0+PE_v_i^=0∩ E_u_j^=1+PE_v_i^=1∩ E_u_j^=0+PE_v_i^=1∩ E_u_j^=1
           =(F-2)^L_Z/F^L_Z+L_Y(F-2)^L_Z-1/F^L_Z+L_X(F-2)^L_Z-1/F^L_Z+L_YL_X(F-2)^L_Z-2/F^L_Z
           =1-OL_Z/F
           ≤ 1-OL_X+L_Y/F
           =1-O1/ϵ^dN.
    
    Now note that the second term in (<ref>) is the bias due to collision of H_2, and similar to (<ref>) it is upper bounded by O1/ϵ^dN. Thus, (<ref>) and (<ref>) give rise to 
    
        ω_ij   = dP_XY/dP_XP_Y(x,y)+μ(ϵ,γ,q,𝐂_XY)+O1/√(N)+ O1/ϵ^dN.
    
    which completes the proof.
    
    
    
    
    
    In the following lemma we make a relation between the bias of an estimator and the bias of a function of that estimator.
    
    
    Assume that g(x): 𝒳→ℝ is infinitely differentiable. If Z is a random variable estimating a constant Z with the bias 𝔹[Z] and the variance 𝕍[Z], then the bias of g(Z) can be written as
    
        g(Z)-g(Z) = ∑_i=1^∞ξ_i𝔹[Z]^i+O√(𝕍[Z]),
    
    where ξ_i are real constants.
    
        gZ-g(Z)    = gZ-g(Z) +gZ-gZ
            = ∑_i=1^∞Z-Z^ig^(i)(Z)/i! + OgZ-gZ
           =∑_i=1^∞ξ_i𝔹[Z]^i+O√(𝕍[Z]).
    
    
    In the second line we have used Taylor expansion for the first term, and triangle inequality for the second term. In the third line we have used the definition ξ_i:=g^(i)(Z)/i!, and the Cauchy Schwarz inequality for the second term.
    
    
    
    In the following we compute the expectation of the first term in (<ref>) and prove Theorem <ref>.
    
     
    
    Recall that N'_i and M'_j respectively are defined as the number of the input points 𝐗 and 𝐘 mapped to the buckets X_i and Y_j using H_1. Similarly, N'_ij is defined as the number of input pairs 𝐗,𝐘 mapped to the bucket pair X_i,Y_j using H_1. Define the notations r(i):=H_2^-1(i) for i∈ℱ and s(x):=H_1(x) for x∈𝒳∪𝒴. Then from (<ref>) since there is no collision of mapping with H_2 into v_i and u_j we have
    
    
        N'_s(x)s(y)N/N'_s(x)N'_s(y)=dP_XY/dP_XP_Y(x,y)+μ(ϵ,γ,q,𝐂_XY)+O1/√(N),
    
    
    
    By using (<ref>) and defining h̃(x)=g̃(x)/x we can simplify the first term of (<ref>) as
    
    
        ∑_i,j∈ℱ PE_ij^≤ 11_E_ijω_iω'_jgω_ij| E_ij^≤ 1   = 1-O1/ϵ^dN∑_i,j∈ℱ1_E_ijω_iω'_jgω_ij| E_ij^≤ 1
           = ∑_i,j∈ℱ1_E_ijN_iM_j/N^2gN_ijN/N_iM_j| E_ij^≤ 1+ O1/ϵ^dN
           = ∑_i,j∈ℱ1_E_ijN'_r(i)M'_r(j)/N^2gN'_r(i)r(j)N/N'_r(i)M'_r(j)+ O1/ϵ^dN
           = ∑_i,j∈ℱ1_E_ijN'_r(i)r(j)/NhN'_r(i)r(j)N/N'_r(i)M'_r(j)+ O1/ϵ^dN
           = 1/N∑_i,j∈ℱN'_r(i)r(j)hN'_r(i)r(j)N/N'_r(i)M'_r(j)+ O1/ϵ^dN
           = 1/N∑_i,j∈ℱ N'_r(i)r(j)hN'_r(i)r(j)N/N'_r(i)M'_r(j)+ O1/ϵ^dN
           = 1/N∑_i=1^N hN'_s(X)s(Y)N/N'_s(X)M'_s(Y)+ O1/ϵ^dN
           = 1/N∑_i=1^N hN'_s(X)s(Y)N/N'_s(X)M'_s(Y)+ O1/ϵ^dN
           = 𝔼_(X,Y)∼ P_XY[𝔼[hN'_s(X)s(Y)N/N'_s(X)M'_s(Y)| X=x,Y=y]]+ O1/ϵ^dN
           = 𝔼_(X,Y)∼ P_XY[dP_XY/dP_XP_Y]+μ(ϵ,γ,q,𝐂_XY)+O1/√(N)+O1/ϵ^dN.
    (<ref>) is due to the fact that N'_r(i)r(j)=0 if there is no edge between v_i and u_j. Also, (<ref>) is due to (<ref>).
    
    From (<ref>) and (<ref>) we obtain
    
    
        I(X,Y)   =∑_e_ij∈ E_Gω_iω'_jgω_ij= 𝔼_(X,Y)∼ P_XY[dP_XY/dP_XP_Y]+μ(ϵ,γ,q,𝐂_XY)+O1/√(N)+O1/ϵ^dN.
    
    Finally using Lemma <ref> results in (<ref>).
    
    
    
    § B. VARIANCE PROOF
    
    
    
    In this section we first prove bounds on the variances of the edge and vertex weights and then we provide the proof of Theorem <ref>.  
    
    
    Under the assumptions A1-A4, the following variance bounds hold true.
    
    
        𝕍[ω_i]   ≤ O1/N,    𝕍[ω'_j]≤ O1/N,    𝕍[ω_ij]≤ O1/N,    𝕍[ν_ij]≤ O1/N.
    
    Here we only provide the variance proof of ω_i. The variance bounds of ω'_j, ω_ij and ν_ij can be proved in the same way. The proof is based on Efron-Stein inequality.
    Define Z_i:=(X_i,Y_i). For using the Efron Stein inequality on 𝐙:=(Z_1,...,Z_N), we consider another independent copy of 𝐙 as  𝐙':=(Z'_1,...,Z'_N) and define 𝐙^(i):=(Z_1,...,Z_i-1,Z'_i,Z_i+1,...,Z_N). Define ω_i(𝐙) as the weight of vertex v_i in the dependence graph constructed by the set 𝐙. By applying Efron Stein inequality <cit.> we have
    
    
        𝕍[ω_i]    ≤1/2∑_i=1^N 𝔼[(ω_i𝐙-ω_i𝐙^(j))^2]
           = 1/2N^2∑_i=1^N 𝔼[(N_i𝐙-N_i𝐙^(j))^2]
           ≤1/2N^2 ON
           ≤ O1/N.
    
    
    In the third line we have used the fact that the absolute value of N_i𝐙-N_i𝐙^(j) is at most 1.
    
    
    
    
    We follow similar steps as the proof of Lemma <ref>. Define I_g(𝐙) as the mutual information estimation using the set 𝐙. By applying Efron Stein inequality we have
    
    
        𝕍[I(X,Y)]    ≤1/2∑_k=1^N 𝔼[(I(𝐙)-I(𝐙^(k)))^2]
           ≤N/2𝔼[∑_e_ij∈ E_Gω_i𝐙ω'_j𝐙gω_ij𝐙-∑_e_ij∈ E_Gω_i𝐙^(k)ω'_j𝐙^(k)gω_ij𝐙^(k)^2]  
           = N/2N^4𝔼[∑_e_ij∈ E_G N_i𝐙M_j𝐙gN_ij𝐙N/N_i𝐙M_j𝐙-∑_e_ij∈ E_G N_i𝐙^(k)M_j𝐙^(k)gN_ij𝐙^(k)N/N_i𝐙^(k)M_j𝐙^(k)^2]  
           ≤1/2N^3𝔼[Σ_n_1+Σ_n_2+Σ_m_1+Σ_m_1+D_n_1m_1+D_n_2m_2^2].
     
    
    Note that in equation (<ref>), when (X_k,Y_k) is resampled, at most two of N_i for i∈ℱ are changed exactly by one (one decrease and the other increase). The same statement holds true for M_j. Let these vertices be v_n_1, v_n_2, v_m_1 and v_m_2. Also the pair collision counts N_ij are fixed except possibly N_n_1m_1 and N_n_2m_2 that may change by one. So, in the fourth line Σ_n_1 and Σ_n_2 account for the changes in MI estimation due to the changes in N_n_1 and N_n_2, and Σ_m_1 and Σ_m_2 account for the changes in M_m_1 and v_m_2, respectively. Finally D_n_1m_1 and D_n_2m_2 account for the changes in MI estimation due to the changes in N_n_1m_1 and N_n_2m_2. For example, Σ_n_1 is precisely defined as follows:
    
        Σ_n_1:=∑_j:e_mj∈ E_G N_mM_jgN_mjN/N_mN_j-(N_m+1)M_jgN_mjN/(N_m+1)M_j
    
    where we have used the notations N_i and N_i^(k) instead of N_i(𝐙) and N_i(𝐙^(k)) for simplicity. Now note that by assumption A4 we have
    
    
        gN_mjN/N_mM_j-gN_mjN/(N_m+1)M_j   ≤ G_g N_mjN/N_mM_j-N_mjN/(N_m+1)M_j
           ≤ ON_mjN/N^2_mM_j.
    
    
    Thus, using (<ref>), Σ_n_1 can be upper bounded as follows
    
    
        Σ_n_1   ≤∑_j:e_mj∈ E_G ON_mjN/N_m^2 =  ON/N_m≤ O(N).
    
    
    It can similarly be shown that N_n_2, Σ_m_1, Σ_m_2, D_n_1m_1 and D_n_2m_2 are upper bounded by O(N). Thus, (<ref>) simplifies as follows
    
    
        𝕍[I(X,Y)]    ≤36O(N^2)/2N^3= O(1/N).
    
    
    § C. OPTIMUM MSE RATES OF EDGE
    
    In this short section we prove Theorem <ref>. 
    
    
    
    
    The proof simply follows by using the ensemble theorem in (<cit.>, Theorem 4) with the parameters ψ_i(t)=t^i and ϕ_i,d(N)=N^-i/2d for the bias result in Theorem <ref>. Thus, the following weighted ensemble estimator (EDGE) can achieve the optimum parametric MSE convergence rate of O(1/N) for q≥ d.
    
    
        I_w:=∑_t∈𝒯w(t)I_ϵ(t),