Geometric Estimation of Multivariate Dependency

This paper proposes a geometric estimator of dependency between a pair of multivariate random variables. The proposed estimator of dependency is based on a randomly permuted geometric graph (the minimal spanning tree) over the two multivariate samples. This estimator converges to a quantity that we call the geometric mutual information (GMI), which is equivalent to the Henze–Penrose divergence. between the joint distribution of the multivariate samples and the product of the marginals. The GMI has many of the same properties as standard MI but can be estimated from empirical data without density estimation; making it scalable to large datasets. The proposed empirical estimator of GMI is simple to implement, involving the construction of an minimal spanning tree (MST) spanning over both the original data and a randomly permuted version of this data. We establish asymptotic convergence of the estimator and convergence rates of the bias and variance for smooth multivariate density functions belonging to a Hölder class. We demonstrate the advantages of our proposed geometric dependency estimator in a series of experiments.


Introduction
Estimation of multivariate dependency has many applications in fields such as information theory, clustering, structure learning, data processing, feature selection, time series prediction and reinforcement learning, see [2] and [3], [4], [5], [6][7][8][9], [10], and [11], respectively.It is difficult to accurately estimate the mutual information in high dimensional settings where the data is multivariate with a Lebesgue continuous density-the setting considered in this paper.An important and ubiquitous measure of dependency is the Shannon Mutual Information (MI), which has seen extensive use across many application domains.However, the estimation of mutual information can often be challenging.In this paper we focus on a measure of MI that we call the Geometric MI (GMI).This MI measure is defined as the asymptotic large sample limit of a randomized minimal spanning tree (MST) statistic spanning the multivariate sample realizations.The GMI is related to a divergence measure called the Henze-Penrose divergence [12], [13], and related to the multivariate runs test [1].In [14], [15], it was shown that this divergence measure can be used to specify a tighter bound for the Bayes error rate for testing if a random sample comes from one of two distributions as compared to previous divergence-type bounds such as the Bhattacharrya bound.Furthermore, the authors of [16] proposed a nonparametric bound on multi-class classification Bayes error rate using a global MST graph.
Let X and Y be random variables with unknown joint density f XY and marginal densities f X and f Y , respectively, and consider two hypotheses: H 0 , X and Y are independent and H 1 , X and Y are dependent, The GMI is defined as the Henze-Penrose divergence between f XY and f X f Y which can be used as a dependency measure.In this paper we prove that for large sample size the randomized minimal spanning tree (MST) statistic spanning the original multivariate sample realizations and randomly shuffled data set in which for given the first element in each pair the second element is replaced by a randomly selected from all second elements, convergences to the GMI measure.A direct implication of [14], [15] is that the GMI provides a tighter bound on the Bayes misclassification rate for the independence test of H 0 vs H 1 .In this paper we propose an estimator based on a random permutation modification of the Friedman-Rafsky multivariate test statistic and show that under certain conditions the GMI estimator achieves the parametric mean square error (MSE) rate when the joint and marginal densities are bounded and smooth.Importantly our proposed GMI estimator does not require explicit estimation of the joint and marginal densities.
In addition to achieving the optimal theoretical MSE rate, computational complexity is an important challenge addressed by researchers in machine learning and data science.Most plug-in based estimators like kernel density estimator (KDE) or K-nearest neighbour (KNN) estimator require runtime complexity of O(n 2 ) which is not suitable for large scale applications.Noshad et al. proposed a graph theoretic direct estimation method based on nearest neighbor ratios (NNR) [17].The NNR estimator is based on k-NN graph and computationally more tractable than other competing estimators with complexity O(kn log n).The construction of the minimal spanning tree lies at the heart of the GMI estimator proposed in this paper.Since the GMI estimator is based on the Euclidean MST the dual-tree algorithm by March et al. [18] can be applied.This algorithm is based on the construction of Bor ůvka [19] and implements the Euclidean MST in approximately O(nlogn) time which is faster than the previous estimators.In this paper we experimentally show that for large number of sample sizes the proposed GMI estimator is less computationally complex than the KDE plug-in method.

Related work
Estimation of mutual information is related to estimation of the information histogram.The most common estimators of MI are based on plug-in density estimation, e.g., using kernel density or kNN density estimators [20], [21].Motivated by ensemble methods applied to divergence estimation [22]- [23], in [21] an ensemble method for combining multiple KDE bandwidths was proposed for estimating MI.Under certain smoothness conditions this ensemble MI estimator was shown to achieve parametric convergence rates.
Another class of estimators of multivariate dependency bypasses the difficult density estimation task.This class includes the statistically consistent estimator of Rényi-α and KL mutual information which is based on the asymptotic of the length of the KNN graph, [24] and [25].The estimator of [26] builds on the KNN methods for Rényi entropy estimation.As MI increases, the dependencies between random variables increase which results in less smooth densities.In [27], it has been shown that for large MI the KNN or KDE approaches are ill-suited candidates for estimating MI since the assumption of local uniformity can be violated when they are strong dependencies.To overcome this issue an assumption on the smoothness of the density is required, see [28], [29], and [22]- [23].For all of these methods, the optimal parametric rate of MSE convergence is achieved when the densities are either d, (d + 1)/2 or d/2 times differentiable [30].In this paper, we assume that joint and marginal densities are smooth in the sense that they belong to a strong Hölder continuous classes of densities Σ s d (η, K), such that the smoothness parameter η ranges in (0, 1]. A MI measure based on the Pearson divergence was considered in [31] that is computational efficient and numerically stable.The authors of [32] and [33] used minimal spanning tree generalized nearest-neighbor graph approaches, respectively, to estimate Rényi mutual information.In [21], a nonparametric mutual information estimator was proposed using a weighted ensemble method with O(1/n) parametric convergence rate.This estimator was based on plug-in density estimation, which is challenging in high dimension.
Our proposed estimator differs from previous methods in the following ways.First, it estimates a different measure of mutual information, the GMI.Second, instead of using the KNN graph the estimator of GMI uses a randomized minimal spanning tree that spans the multivariate realizations.The proposed GMI estimator is motivated by the multivariate runs test of Friedman and Rafsky (FR) [34] which is a multivariate generalization of the univariate Smirnov maximum deviation test [35] and the Wald-Wolfowitz [36] runs test in one dimension.We also emphasize that our proposed GMI estimator does not require boundary correction, in contrast to other graph-based estimators, such as, the NNR estimator [17], scalable MI estimator [37], or cross match statistic [38].

Contribution
The contribution of this paper has three components (1) We propose a novel non-parametric multivariate dependency measure, referred to as Geometric Mutual Information (GMI), that is based on graph-based divergence estimation.The geometric mutual information is constructed using a minimal spanning tree and is a function of Friedman-Rafsky multivariate test statistic.(2) We establish properties of the proposed dependency measure analogous to those of standard mutual information, such as, convexity, concavity, chain rule, and the data processing inequality.(3) We derive a bound on the MSE rate for the proposed geometric estimator.An advantage of the estimator is that it achieves the optimal MSE rate without explicitly performing boundary correction, which is required for most plug-in estimators.

Organization
The paper is organized as follows.In Section 2, we define the geometric mutual information and establish some of its mathematical properties.In Subsection 2.2 and 2.3, we introduce a statistically consistent GMI estimator and derive a bound on its mean square error convergence rate.In Section 3 we verify the theory through experiments.
Throughout the paper, we denote statistical expectation by E, the variance by abbreviation Var.Bold face type indicates random vectors.All densities are assumed to be Lebesgue continuous with no atoms.

The geometric mutual information (GMI)
In this section, we first review the definition of Henze-Penrose (HP) divergence measure defined by Berisha and Hero in [14] and [1].The Henze-Penrose divergence between densities f and g with domain R d for parameter p ∈ (0, 1) is defined as follows (see [1], [39], and [40]): where q = 1 − p.This functional is an f -divergence [41], also known as an Ali-Silvey distance [42].
That is, it satisfies the properties of non-negativity, monotonicity, and joint convexity [15].This measure takes values in [0, 1] and D p = 0 if and only if f = g almost surely.
The geometric mutual information measure is defined as follows.Let f X , f Y , and f XY be the marginal and joint distributions, respectively, of random vectors X ∈ R d x , Y ∈ R d y where d x and d y are positive integers.Then by using (1), a Henze-Penrose generalization of the mutual information between X and Y, is defined by (2) We will show below that I p (X; Y) has a geometric interpretation in terms of the large sample limit of a minimal spanning tree spanning n sample realizations of X ∪ Y. Thus we call I p (X; Y) the geometric mutual information (GMI) between X and Y.The GMI satisfies similar properties to other definitions of mutual information, such as Shannon and Rényi mutual information.Recalling (3) in [14], an alternative form of I p is given by where , and The affinity A p (X; Y) is called the geometric affinity between X and Y, and it was originally defined in [1].The next subsection of the paper is dedicated to the basic inequalities and properties of the proposed GMI measure (2).

Properties of the geometric mutual information
In this subsection we establish basic inequalities and properties of the GMI, I p , given in (2).

Basic inequalities
The following theorem shows that I p (X; Y) is a concave function in f X and a convex function in f Y|X .The proof is given in Appendix A, Subsection 5.1.
Theorem 1. Denote by Ĩp ( f XY ) the GMI I p (X; Y) when X ∈ R d x and Y ∈ R d y have joint density f XY .Then the GMI satisfies (i) Concavity in f X : Let f Y|X be conditional density of Y given X and let g X and h X be densities on R d x .For p ∈ (0, 1) Ĩp The inequality is strict unless either λ 1 or λ 2 are zero or h X = g X .
(ii) Convexity in f Y|X : Let g Y|X and h Y|X be conditional densities of Y given X and let f X be marginal density.Then for p ∈ (0, 1) The equality occurs when either λ 1 or λ 2 are zero or h Y|X = g Y|X .
The GMI, I p (X; Y), satisfies a property analogous to the standard chain rule and the data processing inequality [43].For random variables X ∈ R d x , Y ∈ R d y , and Z ∈ R d z we define the conditional GMI measure by I p (X; Y|Z) = E Z I p (X; Y|Z = z) , where The next theorem establishes a relation between the joint and conditional geometric mutual information.
Theorem 2. For given d-dimensional random vector X with components X 1 , X 2 , . . ., X d and random variable Y, Here f is the joint PDF of random vector (X 1 , X 2 , Y) and note that 0 ≤ δ X 2 ,Y|X 1 ≤ 1.Then we have Furthermore . ., X i .Define a general form of δ in (i) by where f is the join PDF of random vector (X 1 , . . ., X d , Y).A general form of ( 8) is given by Note that in the special case where Further, if both X → Y → Z and X → Z → Y together hold true, we have I p (Y; X) = I p (Z; X).The inequality in (11) becomes a tighter inequality interpretable as a standard data-processing inequality I p (Y; X) ≥ I p (Z; X), when

The Friedman-Rafsky Estimator
Let a random sample {x i , y i } n i=1 from f XY (x, y) be available.Here we show that the GMI I p (X; Y) can be directly estimated without estimating the densities.The estimator is inspired by the MST construction of [34] that provides a consistent estimate of the Henze-Penrose divergence [14], [15].We denote by z i the i-th joint sample x i , y i and by Z n the sample set {z i } n i=1 .Divide the sample set Z n into two subsets Z n and Z n with the proportion α = n /n and β = n /n, where α + β = 1.
Denote by Z n the set This means that for each z ik = (x ik , y ik ) ∈ Z n given the first element x ik the second element y ik is replaced by a randomly selected y ∈ {y jk } n j=1 .This results in a shuffling of the binary relation relating y ik in y jk .The estimator of I p (X; Y) is derived based on the Friedman-Rafsky (FR) multivariate runs test statistic [34] on the concatenated data set , Z n ∪ Z n .The FR test statistic is defined as the number of edges in the MST spanning the concatenated data set that connect a point in Z n to a point in Z n .This test statistic is denoted by R n ,n := R n ,n (Z n , Z n ).Note that since the MST is unique with probability one (under the assumption that all density functions are Lebesgue continuous) then all inter point distances between nodes are distinct.This estimator converges to I p (X; Y) almost surely as n → ∞.The procedure is summarized in Algorithm 1. Theorem 3 shows that the output in Algorithm Algorithm 1 FR estimator of GMI Theorem 3.For given proportionality parameter α ∈ (0, 1), choose n , n such that n + n = n and as n → ∞, we have n /n → α and n Note that the asymptotic limit in (12) depends on the proportionality parameter α.Later in Subsection 2.4, we discuss on choosing an optimal parameter α.In Fig. 1, we show a visualization of the MST constructed over merged independent (ρ = 0) and highly dependent (ρ = 0.9) data sets drawn from two dimensional standard Normal distribution with correlation coefficient ρ.Notice that the MST connects different colored samples, corresponding to independent and dependent samples, respectively, indicated in green.The total number of green edges is the FR test statistic.
The MST and FR statistic of spanning the merged set of Normal points when X and Y are independent (denoted in blue points) and when X and Y are highly dependent (denoted in red points).
The FR test statistic is the number of edges in the MST that connect samples from different color nodes (denoted in green) and it is used to estimate the GMI I p .

Convergence Rates
In this subsection we provide the MSE convergence rates in the form of upper bounds on the bias and the variance.This MSE bound is given in terms of the sample size n, the dimension d, and the proportionality parameter α.Deriving convergence rates for mutual information estimators has been of interest in information theory and machine learning [26], [21].The rates are typically derived in terms of a smoothness condition on the densities, such as the Hölder condition [44].Here we assume f X , f Y and f XY with support sets S X , S Y , and S XY := S X × S Y , respectively, are smooth in the sense that they belong to strong Hölder continuous classes of densities Σ s d (η, K), 0 < η ≤ 1 [45], [44]: (Strong Hölder class): Let X ⊂ R d be a compact space.The strong Hölder class Σ s d (η, K), with η-Hölder parameter and constant K, of functions on L d -norm, consists of the functions g that satisfy where p k x (z) is the Taylor polynomial (multinomial) of g of order k expanded about the point x and η is defined as the greatest integer strictly less than η.Note that for the standard Hölder class the factor g(x) on the RHS of ( 13) is omitted.
To find the optimal parameter α we require explicit expressions for the bias and variance.Bounds are provided in Appendix E, Subsection 5.5.To obtain such expressions, we require several assumptions on the absolutely continuous densities f X , f Y , f XY and support sets S X , S Y , S XY : (A.1)Each of the densities belong to Σ s d (η, K) class in their support sets with smoothness parameters η and smoothness constants K.
(A.3)All densities are bounded on their support sets, i.e. there exist two sets of constants

Remarks:
1.The assumption (A.3) implies that the strong Hölder class, (A.1), is a subset of standard Hölder class: XY are equal i.e. the density functions are uniform, g(x) = K , and if K = K.K then two classes ( 13) and ( 14) are equal.3.By using the Strong Hölder class (13) our proposed bias bound in Theorem 4 below, becomes tighter.Another bound can be derived for the standard Hölder class ( 14) but the proposed bias bound becomes weaker, and so in the result for bounding convergence rate.
Note that according to Theorem 13 in [46], the constant c d is lower bounded by is the binary entropy i.e.
A proof of Theorem 4 is given in Appendix E. The next theorem gives an upper bound on the variance of the FR estimator R n ,n .The proof of the variance result requires a different approach than the bias bound.The proof uses the Efron-Stein inequality [47] and is similar to arguments in [48, Appendix C], and is omitted.In Theorem 5 we assume that the densities f X , f Y , and f XY are absolutely continuous and bounded.Note that we can obtain a weaker bound than (15) when we relax (A.1) to the standard Hölder class with a constant coefficient K instead of Kg(x).
where c d is a constant depending only on the dimension d.

Minimax Parameter α
A question in Algorithm 1 is the existance of a best proportion α that partitions sample.In this subsection we obtain on an upper bound of the MSE rate and characterize the parameter α that minimizes the upper bound.This optimal α depends on maximize densities, however.
Recall assumptions (A.1), (A.2), and (A.3) in Subsection 2.3.The constant α can be chosen to ensure that the MSE converges rate obtained from the bias and variance rates ( 15) and ( 5) by selecting α to minimize the maximum MSE, where the maximum is taken over the space of Hölder smooth joint densities f XY .
Throughout this subsection we use the following notations: , where η is the smoothness parameter, Consider the following optimization problem: where and Note that in (20), c, c , c 1 , c 2 are constants, and c d only depends on the dimension d.Also, in (21), a i and b i are constants.Let * XY be the optimal XY i.e. * XY be the solution of the optimization problem (18).Set , the optimal choice of XY in terms of maximizing the MSE is * XY = C U and the saddle point for the parameter α, denoted by α, is given as follows: Further details are given in Appendix F.

Simulation Study
In this section numerical simulations are presented that illustrate the theory in Section 2. We perform multiple experiments to demonstrate the utility of the proposed direct estimator of the HP-divergence in terms of the dimension d and the sample size n.Our proposed MST-based estimator of the GMI is compared to other plug-in GMI estimators, in particular the standard KDE estimator of [21], where the convergence rates of Theorem 4 and 5 are validated.We use multivariate Normal simulated data in the experiments.In this section, we also discuss the choice of the proportionality parameter α and compare runtime of our proposed FR estimator approach with KDE method.
Here we perform four sets of experiments to illustrate the estimator and the theory derived above.For the first set of experiments the MSE of the FR estimator proposed in Algorithm 1 is shown in Fig. 2-(left).The samples were drawn from d-dimensional Normal distribution, with various sample sizes and dimensions d = 6, 10, 12.We selected the proportionality parameter α = 0.3 and computed the MSE in terms of the sample size n.We show the log-log plot of MSE when n varies in [100,1500].Note that the empirically optimal proportion α depends on n, so to avoid the computational complexity we fixed α for this experiment.The experimental result shown in Fig. 2-(left) validates the theoretical MSE growth rates derived from ( 15) and ( 16), i.e., decreasing sublinearly in n and increasing exponentially in d.In Fig. 2-(right), we compare our estimator with the true GMI and the standard KDE [21].For the KDE approach, we estimated the joint and marginal densities and then plugged them into the proposed expression (3) in [21].The bandwidth used for the the KDE plug-in estimator was selected by setting h = n −1/(d+1) to minimize the MSE of the plug in estimator.We generated data from the two dimensional Normal distribution with zero mean and covariance matrix We varied ρ in range [0.1, 0.9] and computed the proposed estimated GMI and the KDE.We used the Monte Carlo method to approximate the integral (2) and compute the GMI measure (labeled True GMI).We see that as ρ increases the estimated GMI using the FR test statistic outperforms the KDE approach (FR curve approaches to True curve and KDE curve diverges from True curve).In this set of experiments α = 0.6.Fig. 3 again compares the FR estimator with the standard KDE estimator.In this setting, we draw samples from the multivariate standard Normal distribution with dimensions d = 4 and d = 12.In both cases the proportionality parameter α = 0.5.The left plots in Fig. 3 show the MSE (100 trials) of the GMI estimator implemented with an KDE estimator (with bandwidth as in Fig. 2 i.e.In both cases the proportionality parameter α is 0.5.The FR estimator in both plots for sample size n in [700, 1600] outperforms the standard kernel plug-in estimator, especially for larger dimensions.h = n −1/(d+1) ) for dimensions d = 4, 12 and various sample sizes.For all dimensions and sample sizes the FR estimator also outperforms the plug-in KDE estimator based on the estimated log-log MSE slope given in Fig. 3 (left plots).The right plots in Fig. 3 show the geometric mutual information estimated by KED and FR approaches.The iteration in this experiment is 100 and error bars are standard deviations.We observe that for higher dimension d = 12 and larger sample size n, the KDE estimator approaches to the true GMI slower than the FR estimator.This reflects the power of graph-based (direct) estimators, and in particular our proposed FR estimator.
The comparison between MSEs for various dimension d is shown in Fig. 4 (left).This experiment highlights higher dimension's impact on our proposed FR estimator for the GMI measure.As expected, for larger sample size n, MSE decreases while for higher dimension it increases.In this setting, we have generated samples from standard Normal distribution with size n ∈ [10 2 , 4 × 10 3 ] and α = 0.5.From Fig. 4 (left) we observe that for larger sample size, MSE curves are ordered based on their corresponding dimensions.Fig. 4 (right) illustrates the MSE vs proportion parameter α when n = 500, 10 4 samples are generated from Normal distribution with ρ = 0.7, 0.5.First, following subsection 2.4, we compute the bound [α L 0 , α U 0 ] and then derive optimal α in this range.Therefore, each experiment with different sample size and ρ provides different range [α L 0 , α U 0 ].We observe that the MSE is not necessarily a monotonic function in α and its behavior strongly depends on sample size n, d, and density functions' bounds.This is studied extensively in Appendix F. In this set of experiments Ξ(α L 0 ) > 0, therefore following the results in subsection 2.4, we have α = α L 0 .The optimal α is indicated in the Fig. 4 (right).Next, we consider three scenarios to analyze parameter α.In these scenarios the lower bounds C L X , C L Y , C L XY and upper bounds C U X , C U Y , C U XY are unknown, therefore results in Section 2.4 are not applicable.In these set of experiments we varied α in the range (0, 1) to divide our original sample.We generated sample from the multivariate standard Normal distribution in all three scenarios (all features are independent).Therefore the true GMI is zero and in all scenarios the GMI column is compared with zero.In each scenario we fixed dimension d and sample size n and varied α = 0.2, 0.5, 0.8.The dimension and sample size in Scenarios 1,2, and 3 are d = 6, 8, 10 and n = 1000, 1500, 2000, respectively.In Table 1 the last column (α) stars the parameter α ∈ {0.2, 0.5, 0.8} with the minimum MSE and GMI (I α ) in each scenario.Table 1 shows that in these sets of experiments when α = 0.5, the GMI estimator has less MSE (i.e. is more accurate) than when α = 0.2 or α = 0.8.This experimentally demonstrates that if we split our training data, the proposed Algorithm 1 performs better than α = 0.2 or 0.8 to estimate the GMI measure.We applied the FR approach to estimate the GMI (I α ) with α = 0.2, 0.5, 0.8.We varied dimension d = 6, 8, 10 and sample size n = 1000, 1500, 2000 in each scenario.We observe that for α = {0.2,0.5, 0.8}, the GMI estimator provides less MSE when α = 0.5.

KDE FR Estimator
Figure 5. Runtime of KDE approach and proposed FR estimator of GMI measure vs sample size.The proposed GMI estimator achieves significant speedup, while for small sample size, the KDE method becomes overly fast.Note that in this experiment we have generated our sample from standard Gaussian distribution.
Finally Fig. 5 shows the runtime (KDE -FR methods) as a function of sample size n.We vary sample size in the range [10 3 , 10 4 ].Observe that for smaller number of samples the KDE method is slightly faster but as n becomes large we see significant relative speedup of the proposed FR method.

Conclusion
In this paper we have proposed a new measure of mutual information, called Geometric MI (GMI), that is related to the Henze-Penrose divergence.The GMI can be viewed as dependency measure that is the limit of the Friedman-Rafsky test statistic that depends on the minimal spanning tree over all data points.We established some properties of the GMI in terms of convexity/concavity, chain rule, and analogous to a data processing inequality.A direct estimator of the GMI was introduced that uses random permutations of observed relationships between variables in the multivariate samples.An explicit form for the MSE convergence rate bound was derived that depends on a free parameter called the proportionality parameter.An asymptotically optimal form for this free parameter was given that minimizes the MSE convergence rate.Simulation studies were performed that illustrate and verify the theory.

Appendices
We organize the appendices as the following: Theorem 1 which establishes convexity/concavity property is proved in Appendix A. Appendix B and C are arranged to establish the inequality ( 8) and (11) for given p ∈ (0, 1), respectively.In Appendix D, we first prove that the set Z n which is randomly generated from original dependent data, contains independent samples asymptotically.Later by using the generated independent sample Z n we show that for given α the FR estimator of the GMI given in Algorithm 1 intends to I α approximately.Appendix E is dedicated to the Theorem 4. A full discussion on the proportionality parameter (α) optimization strategy is provided in Appendix F.
Notice this follows by using the convex function u(y) = y 2 (p + q y) for any p ∈ (0, 1), q = 1 − p, and the Jensen inequality.

Define the shorthand
x , y , and xy for dx, dy and dxdy, respectively.To prove part (i) of the Theorem 1, we represent the LHS of ( 5) as: Furthermore, the RHS of ( 5) can be rewritten as Thus, in order to prove LHS ≥ RHS, we use the inequality below: In Lemma 1, let and for i = 1, 2, Then the claimed assertion (i) is obtained.Part (ii) follows by convexity of D p and the following expression: Therefore the claim in (6) is proved.

Appendix B: Theorem 2
Proof.We prove part (i) and the second part (ii) is shown by repetition method.To show (8), we use the inequality below.Given p ∈ (0, 1) and q = 1 − p, we can easily check that for positive t > q, s > q, such that t, s = 1: This implies Consequently Here f is the joint PDF of random vector (X 1 , X 2 , Y).This completes the proof of part (i).

Appendix C: Proposition 1
Proof.Recall the Theorem 2, part (i).First from X → Y → Z we have f XYZ = f XY f Z|Y and then by applying the Jensen inequality, we can write where (x|z) .
Now by Markovian property we can immediately simplify the last line of ( 26) to the RHS in (11).Furthermore we can easily show that if X → Z → Y, we have f XYZ = f ZX f Y|Z and therefore I p (Z; X) = I p (X; Y, Z).This together with (26) proves that under both conditions X → Y → Z and X → Z → Y, the equality I p (X; Y) = I p (Z; X) holds true.

Appendix D: Theorem 3
Proof.In this appendix, we first derive two required Lemmas 2 and 3 below: Lemma 2. Consider random vector Z = (X, Y) with joint probability density function (pdf) f XY .Let Z n = {z 1 , . . ., z n } = {(x i , y i )} n i=1 be a set of samples with pdf f XY .Let Z n and Z n be two distinct subsets of Z n such that n + n = n and sample proportion is α = n /n and β = 1 − α.Next, let Z n = { z 1 , . . ., z n } be a set of pairs such that z k = (x i k , y j k ), k = 1, . . ., n are selected at random from Z n .Denote Z = ( X, Y) as the random vector corresponding to samples in Z n .Then as n → ∞ such that n also grows in a linked manner that β = 0 then the distribution of Z convergences to f X × f Y i.e. random vectors X and Y become independent.
Proof.Consider two subsets A, B ⊂ R n , then we have Here I A stands for the indicator function.Note that and X i and Y j , i = j are independent, therefore On the other hand, we know that n = β n, so we get From ( 28), we observe that when β takes larger values the bound becomes tighter.So if n → ∞ such that n also becomes large enough in a linked manner that β = 0 then the RHS in (28) tends to zero.This implies that X and Y become independent when n → ∞.
An immediate result of Lemma 2 is the following: Lemma 3.For given random vector Z n = (X n , Y n ) from joint density function f XY and with marginal density functions f X and f Y set Z n = z 1 , . . ., z n be realization of random vector Z as in Lemma 2 with parameter β = n /n.Then for given points of Z n at z = ( x, y), we have Now, let us get back to our main goal which is the proof of assertion (12).Consider two subsets Z n and Z n as described in Subsection 2.2 .Assume that the components of sample Z n follow density function f X Y .Therefore by owing to Lemma 2 and 3, when n → ∞ then f X Y → f X f Y .Let M n and N n be Poisson variables with mean n and n independent of one another and {Z i } and { Z j }.Assume two Poisson processes Z n = Z 1 , . . ., Z M n and Z n = Z 1 , . . ., Z N n , with the FR statistic R n ,n .So by owing to the arguments in [1], [48] we shall prove the following: , . . .be independent variables with common density  [49], [1], has a same distribution as R n ,n .Given points of F n ,n at z = (x , y ) and z = (x , y ), the probability that they have different marks is given by (30). Set then Now recall (31).We observe that g n ,n (z , z ) → g(z , z ).Going back to (32), we can write Consider the non-Poisson process So by the fact that Introduce Then φ n ,n (x, y) → φ(x, y) uniformly because of n /n → α and n /n → β.Thus by using Proposition 1 in [1], we get So, we conclude the proof.

Appendix E: Theorem 4
Proof.We begin by providing a family of bias rate bound for the FR estimator R n ,n in terms of a parameter l.Then by plugging the optimal l, we prove the bias rate bound given in (15).Theorem 6.Let R n ,n := R(Z n , Z n ) be the FR test statistic.Then a bound on the bias rate of the R n ,n estimator for 0 < η ≤ 1, d ≥ 2 is given by where 0 < η ≤ 1 is the Hölder smoothness parameter and c d is a constant depending only on d.Set where A more explicit form for the bound on the RHS is given below: Proof.Consider two Poisson variables M n and N n with mean n and n respectively and independent of one another and {Z i } and { Z j }.Let Z n and Z n be the Poisson processes {Z 1 , . . ., Z M n } and Applying Lemma 1, and ( 12) in [1], we can write Here c d denotes the largest possible degree of any vertex of the MST in R d .Following the arguments in [48], we have Next let n i and n i be independent binomial random variables with marginal densities B(n , a i l −d ) and B(n , b i l −d ) such that a i , b i are non-negative constants a i ≤ b i and Therefore using subadditivity property in Lemma 2.2, [48], we can write where M = l d , and η > 0 stands with Hölder smoothness parameter.Further, for given n i , n i , let , . . .be independent variables with common densities for (x, y) ∈ R d × R d : Here z = (x , y ), z = (x , y ), and P n i ,n i (z , z ) is given in below: By owing to the Lemma B.6 in [48] and applying the analogous arguments, we can write the expression in ( 43): where Going back to Lemma 3, we know that Therefore the first term in RHS of (43) turns to be less and equal than dx dy, and the second term is less and equal than where Recall the definition of the dual MST and FR statistic denoted by R * n ,n from [48]: Definition: (Dual MST, MST * and dual FR statistic R * m,n ) Let F i be the set of corner points of the subsection Q i for 1 ≤ i ≤ l d .Then we define MST * (X m ∪ Y n ∩ Q i ) as the boundary MST graph of partition Q i [50], which contains X m and Y n points falling inside the section Q i and those corner points in F i which minimize total MST length.Notice it is allowed to connect the MSTs in Q i and Q j through points strictly contained in Q i and Q j and corner points are taking into account under condition of minimizing total MST length.Another word, the dual MST can connect the points in Q i ∪ Q j by direct edges to pair to another point in Q i ∪ Q j or the corner the corner points (we assume that all corner points are connected) in order to minimize the total length.To clarify this, assume that there are two points in Q i ∪ Q j , then the dual MST consists of the two edges connecting these points to the corner if they are closed to a corner point otherwise dual MST consists of an edge connecting one to another.Further, R * m,n (X m , Y n ∩ Q i ) is defined as the number of edges in MST * graph connecting nodes from different samples and number of edges connecting to the corner points.Note that the edges connected to the corner nodes (regardless of the type of points) are always counted in dual FR test statistic R * m,n .See [48] for more detail about the dual MST and dual FR test statistic.Similar discussion as above and in [48] The first term of RHS in (44)  Here A β,n /n f ,n (x, y) is given as (37) by substituting n /n in α such that β = 1 − α.Hence, the proof of Theorem 6 is completed.
Going back to the proof of (15), without loss of generality assume that (n)l −d > 1.In the range d ≥ 2 and 0 < η ≤ 1.We select l as a function of n and β to be the sequence increasing in n which minimizes the maximum of these rates: The solution l = l(n, β) occurs when l d (n) −η/d = l −ηd , or equivalently l = (n) η/(d 2 (η+1)) and also l d β −1 n −1 = l −ηd which implies l = (β n) 1/(d(1+η)) .Substitute this in l in the bound (36), the RHS expression in (15) for d ≥ 2 is derived.

i=1 1 :
Find α using arguments in Subsection 2.4 2: n ← αn, n ← (1 − α)n 3: Divide Z n into two subsets Z n and Z n 4: Z n ← (x ik , y jk ) n k=1 : shuffle first and second elements of pairs in Z n 5: Z ← Z n ∪ Z n 6: Construct MST on Z 7: R n ,n ← # edges connecting a node in Z n to a node of Z n 8: I p ← 1 − R n ,n n + n 2n n Output: I p , where p = α 1 estimates the GMI with parameter p = α.The proof is provided in Appendix D, Subsection 5.4.

Figure 2 .
Figure 2. (left) Log-log plot of theoretical and experimental MSE of the FR estimator of the GMI as a function of sample size n for d = 6, 10, 12 and fixed smoothness parameter η. (right) The GMI estimator using two approaches, FR test statistic and KDE method along with True GMI.In this experiment, we generated data from the two dimensional Normal distribution with zero mean and covariance matrix (23) for various ρ ∈ [0.1, 0.9].

Figure 3 .
Figure 3. MSE log-log plots as a function of sample size n (left) for the proposed Friedman-Rafsky estimator ("Estimated GMI") and the standard KDE plug-in estimator ("KDE").The two right plots each corresponds to the geometric mutual information with the dimension d = 4 (top) and d = 12 (bottom).In both cases the proportionality parameter α is 0.5.The FR estimator in both plots for sample size n in [700, 1600] outperforms the standard kernel plug-in estimator, especially for larger dimensions.

5 Figure 4 .
Figure 4. MSE log-log plots as a function of sample size n for the proposed FR estimator.We compare the MSE of our proposed FR estimator for various dimensions d = 15, 20, 50 (left).As d increases, the blue curve takes larger values than green and orange curves i.e.MSE increases as d grows.However, this is more evidential for large sample size n.The second experiment (right) focuses on optimal proportion α for n = 500, 10 4 and ρ = 0.7, 0.5.α is the optimal α for α ∈ [α L 0 , α U 0 ].
, consider the Poisson processes samples and the FR test statistic under the union of samples, denoted by R * n ,n , and superadditivity of dual R * n ,n , we have−d 2n n f XY (x, y) f X (x) f Y (y) − δ f /(βn) n f XY (x, y) + n f X (x) f Y (y) − δ f /(βn)

Table 1 :
Comparison between different scenarios of various dimensions and sample sizes in terms of parameter α.
Let L n ,n be an independent Poisson variable with mean n + n .Let F n ,n = + n f X Y (x, y) , and mark 2 otherwise, having the FR test statistic R n ,n which by the marking theorem Denote L n i ,n i be an independent Poisson variable with mean n i + n i and F n i ,n i = W Poisson of rate n i f XY + n i f X Y .Let F n i ,n i be the non-Posisson point process .Assign a mark from the set {1, 2} to each points of F n i ,n i .Let Z n i be the sets of points marked 1 with each probability n i f XY (x, y) n i f XY (x, y) + n i f X Y (x, y) and let Z n i be the set points with mark 2. Note that owing to the marking theorem [49], Z n i and Z n i are independent Poisson processes with the same distribution as Z n i and Z n i , respectively.Considering R n i ,n i is greater and equal than2n n f XY (x, y) f X (x) f Y (y) n f XY (x, y) + n f X (x) f Y (y) dxdy − δ f βn 2n n f XY (x, y) n f XY (x, y) + n f X (x) f Y (y)dxdy. is the largest possible degree of any vertex of the MST in R d as before.Consequently, we haveE R n ,n n − 2αβ f XY (x, y) f X (x) f Y (y) α f XY (x, y) + β f X (x) f Y (y) dxdy ≤ B(α, f XY , f X f Y ), E n i ,n i (n i + n i ) ς η (l, n i , n i ) ≤ (a i ) −1 2 f XY (x, y) f X (x) f Y (y) + O(δ f (βn))n f XY (x, y) + n f X (x) f Y (y) dxdy.
(45)e B is defined in(39)and A β,α f ,n (x, y) has been introduced in (37).The last line in(45)is implied from the fact that M ∑ i=1