Next Article in Journal
Pricing Interval European Option with the Principle of Maximum Entropy
Previous Article in Journal
Entropy-Based Clustering Algorithm for Fingerprint Singular Point Detection
Previous Article in Special Issue
Gender Diversity in STEM Disciplines: A Multiple Factor Problem
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Geometric Estimation of Multivariate Dependency

by
Salimeh Yasaei Sekeh
*,† and
Alfred O. Hero
Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI 48109, USA
*
Author to whom correspondence should be addressed.
Current address: School of Computing and Information Science, University of Maine, Orono, ME 04469, USA.
Entropy 2019, 21(8), 787; https://doi.org/10.3390/e21080787
Submission received: 21 May 2019 / Revised: 2 August 2019 / Accepted: 8 August 2019 / Published: 12 August 2019
(This article belongs to the Special Issue Women in Information Theory 2018)

Abstract

:
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.

1. 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 [1,2,3,4,5,6,7,8,9,10], respectively. It is difficult to accurately estimate the mutual information in high-dimensional settings, specially where the data is multivariate with an absolutely continuous density with respect to Lebesgue measure—the setting considered in this paper. An important and regular 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 [11,12], and related to the multivariate runs test [13]. 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 the bound in [14,15] is tighter than previous divergence-type bounds such as the Bhattacharrya bound [16]. Furthermore, the authors of [17] proposed a non-parametric 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 X Y 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,
H 0 : f X Y = f X f Y , versus H 1 : f X Y f X f Y .
The GMI is defined as the Henze–Penrose divergence between f X Y and f X f Y which can be used as a dependency measure. In this paper, we prove that for large sample size n the randomized MST statistic spanning the original multivariate sample realizations and a randomly shuffled data set converges almost surely 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 optimal test of independence. 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 density is bounded and smooth. Importantly unlike other measures of MI, our proposed GMI estimator does not require explicit estimation of the joint and marginal densities.
Computational complexity is an important challenge in machine learning and data science. Most plug-in-based estimators, such as the kernel density estimator (KDE) or the K-nearest-neighbor (KNN) estimator with known convergence rate, 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) [18]. The NNR estimator is based on k-NN graph and computationally more tractable than other competing estimators with complexity O ( k n 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. [19] can be applied. This algorithm is based on the construction of Borůvka [20] and implements the Euclidean MST in approximately O ( n l o g n ) time. In this paper, we experimentally show that for large sample size the proposed GMI estimator has faster runtime than the KDE plug-in method.

1.1. Related Work

Estimation of mutual information has a rich history. The most common estimators of MI are based on plug-in density estimation, e.g., using the histogram, kernel density or kNN density estimators [21,22]. Motivated by ensemble methods applied to divergence estimation [23,24], in [22] 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 estimators of Rényi- α and KL mutual information which are motivated by the asymptotic limit of the length of the KNN graph, [25,26] when joint density is smooth. The estimator of [27] builds on KNN methods for Rényi entropy estimation. The authors of [26], showed that when MI is large the KNN and KDE approaches are ill-suited for estimating MI since the joint density may be insufficiently smooth when there are strong dependencies. To overcome this issue an assumption on the smoothness of the density is required, see [28,29], and [23,24]. For all 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 Hölder continuous classes of densities Σ d ( η , K ) , where the smoothness parameter η ( 0 , 1 ] and the Lipschitz constant  K > 0 .
A MI measure based on the Pearson chi-square divergence was considered in [31] that is computational efficient and numerically stable. The authors of [27,32] used nearest-neighbor graph and minimal spanning tree approaches, respectively, to estimate Rényi mutual information. In [22], a non-parametric 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 dependency 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) [33] which is a multivariate generalization of the univariate Smirnov maximum deviation test [34] and the Wald-Wolfowitz [35] runs test in one dimension. We also emphasize that the proposed GMI estimator does not require boundary correction, in contrast to other graph-based estimators, such as, the NNR estimator [18], scalable MI estimator [36], or cross match statistic [37].

1.2. 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), which is based on graph-based divergence estimation. The geometric mutual information is constructed using a minimal spanning tree and is a function of the Friedman–Rafsky multivariate test statistic.
(2)
We establish properties of the proposed dependency measure analogous to those of Shannon mutual information, such as, convexity, concavity, chain rule, and a type of 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 the need for boundary correction, which is required for most plug-in estimators.

1.3. Organization

The rest of the paper is organized as follows. In Section 2, we define the geometric mutual information and establish some of its mathematical properties. In Section 2.2 and Section 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 and the variance by abbreviation Var . Bold face type indicates random vectors. All densities are assumed to be absolutely continuous with respect to non-atomic Lebesgue measure.

2. The Geometric Mutual Information (GMI)

In this section, we first review the definition of the Henze–Penrose (HP) divergence measure defined by Berisha and Hero in [13,14]. The Henze–Penrose divergence between densities f and g with domain R d for parameter p ( 0 , 1 ) is defined as follows (see [13,14,15]):
D p ( f , g ) = 1 4 p q p f ( x ) q g ( x ) 2 p f ( x ) + q g ( x ) d x ( p q ) 2 ,
where q = 1 p . This functional is an f-divergence [38], equivalently, as an Ali-Silvey distance [39], i.e., it satisfies the properties of non-negativity, monotonicity, and joint convexity [15]. The measure (1) takes values in [ 0 , 1 ] and D p ( f , g ) = 0 if and only if f = g almost surely.
The mutual information measure is defined as follows. Let f X , f Y , and f X Y 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
I p ( X ; Y ) = D p ( f X Y , f X f Y ) = 1 4 p q p f X Y ( x , y ) q f X ( y ) f Y ( y ) 2 p f X Y ( x , y ) + q f X ( x ) f Y ( y ) d x d y ( p q ) 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 the merged labeled samples X Y . Thus, we call I p ( X ; Y ) the 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
I p ( X ; Y ) = 1 A p ( X ; Y ) = u p ( X ; Y ) 4 p q ( p q ) 2 4 p q ,
where
A p ( X ; Y ) = f X Y ( x , y ) f X ( x ) f Y ( y ) p f X Y ( x , y ) + q f X ( x ) f Y ( y ) d x d y = E X Y p f X Y ( X , Y ) f X ( X ) f Y ( Y ) + q 1 , and u p ( X ; Y ) = p f X Y ( x , y ) q f X ( x ) f Y ( y ) 2 p f X Y ( x , y ) + q f X ( x ) f Y ( y ) d x d y = 1 4 p q A p ( X ; Y ) .
The function A p ( X ; Y ) was defined in [13] and is called the geometric affinity between X and Y . The next subsection of the paper is dedicated to the basic inequalities and properties of the proposed GMI measure (2).

2.1. Properties of the Geometric Mutual Information

In this subsection we establish basic inequalities and properties of the GMI, I p , given in (2). 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.1.
Theorem 1.
Denote by I ˜ p ( f X Y ) the GMI I p ( X ; Y ) when X R d x and Y R d y have joint density f X Y . 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 . Then for λ 1 , λ 2 [ 0 , 1 ] , λ 1 + λ 2 = 1
I ˜ p λ 1 f Y | X g X + λ 2 f Y | X h X λ 1 I ˜ p ( f Y | X g X ) + λ 2 I ˜ p ( f Y | X h X ) .
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 λ 1 , λ 2 [ 0 , 1 ] , λ 1 + λ 2 = 1
I ˜ p λ 1 g Y | X f X + λ 2 h Y | X f X λ 1 I ˜ p ( g Y | X f X ) + λ 2 I ˜ p ( h Y | X f X ) .
The inequality is strict unless either λ 1 or λ 2 are zero or h Y | X = g Y | X .
The GMI, I p ( X ; Y ) , satisfies properties analogous to the standard chain rule and the data-processing inequality [40]. For random variables X R d x , Y R d y , and Z R d z with conditional density f X Y | Z we define the conditional GMI
I p ( X ; Y | Z ) = E Z I ˜ p ( f X Y | Z ) , where I ˜ p ( f X Y | Z ) = 1 f X Y | Z ( x , y | z ) f X | Z ( x | z ) f Y | Z ( y | z ) p f X Y | Z ( x , y | z ) + q f X | Z ( x | z ) f Y | Z ( y | z ) d x d y .
The next theorem establishes a relation between the joint and conditional GMI.
Theorem 2.
For given d-dimensional random vector X with components X 1 , X 2 , , X d and random variable Y,
I p ( X ; Y ) I p ( X 1 ; Y ) i = 1 d 1 1 I p ( X i ; Y | X i 1 ) ,
where X i : = X 1 , X 2 , , X i and the conditional GMI I p ( X i ; Y | X i 1 ) is defined in (7).
For d = 2 Theorem 2 reduces to
I p ( X 1 , X 2 ; Y ) I p ( X 1 ; Y ) 1 I p ( X 2 ; Y | X 1 ) ,
Please note that when i = 1 d 1 1 I p ( X i ; Y | X i 1 ) 1 , the inequality (8) is trivial since 0 I p ( X 1 ; Y ) 1 . The proof of Theorem 2 is given in Appendix A.2. Theorem 2 is next applied to the case where X and Y form a Markov chain. The proof of the following “leany” data-processing inequality (Proposition 1) is provided in Appendices section, Appendix A.3.
Proposition 1.
Suppose random vectors X , Y , Z form a Markov chain denoted, X Y Z , in the sense that f X Y Z = f X | Y f Y | Z f Z . Then for p ( 0 , 1 )
I p ( Y ; X ) I p ( Z ; X ) p E X Y δ X , Y + ( 1 p ) 1 ,
where
δ X , Y = f X | Y ( X | Y ) f Z | Y ( z | Y ) f X | Z ( X | z ) d z .
Furthermore, 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 (10) becomes interpretable as the standard data-processing inequality I p ( Y ; X ) I p ( Z ; X ) , when
E Z f ( Z | Y ) f ( Z | X ) = ,
since
E X Y δ X , Y = E X Y f ( X | Y ) f ( X ) E Z f ( Z | Y ) f ( Z | X ) .

2.2. The Friedman–Rafsky Estimator

Let a random sample { x i , y i } i = 1 n from f X Y ( 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 [33] 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 } i = 1 n . 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
( x i k , y j k ) , k = 1 , , n , selected at random from Z n :
Entropy 21 00787 i001
This means that for each z i k = ( x i k , y i k ) Z n given the first element x i k the second element y i k is replaced by a randomly selected y { y j k } j = 1 n . This results in a random shuffling of the binary relation relating y i k in y j k . The estimator of I p ( X ; Y ) is derived based on the Friedman–Rafsky (FR) multivariate runs test statistic [33] 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 merged 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 ) . Please 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.
Algorithm 1: MST-based estimator of GMI
Input: Data set Z n : = ( x i , y i ) i = 1 n
   1: Find α ˜ using arguments in Section 2.4
   2: n α ˜ n , n ( 1 α ˜ ) n
   3: Divide Z n into two subsets Z n and Z n
   4: Z ˜ n { ( x i k , y j k ) k = 1 n : 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 2 n n
Output: I ^ p , where p = α ˜
Theorem 3 shows that the output in Algorithm 1 estimates the GMI with parameter p = α . The proof is provided in Appendix A.4.
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 / n β = 1 α . Then
1 R n , n ( Z n , Z ˜ n ) n 2 n n I α ( X ; Y ) , a . s .
Please note that the asymptotic limit in (11) depends on the proportionality parameter α . Later in Section 2.4, we discuss the choice of an optimal parameter α ˜ . In Figure 1, we illustrate the MST constructed over merged independent ( ρ = 0 ) and highly dependent ( ρ = 0.9 ) data sets drawn from two-dimensional normal distributions with correlation coefficients ρ . Notice that the edges of the MST connecting samples with different colors, corresponding to independent and dependent samples, respectively, are indicated in green. The total number of green edges is the FR test statistic R n , n ( Z n , Z ˜ n ) .

2.3. Convergence Rates

In this subsection we characterize the MSE convergence rates of the GMI estimator of Section 2.2 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 [22,27]. The rates are typically derived in terms of a smoothness condition on the densities, such as the Hölder condition [41]. Here we assume f X , f Y and f X Y have support sets S X , S Y , and S X Y : = S X × S Y , respectively, and are smooth in the sense that they belong to Hölder continuous classes of densities Σ d s ( η , K), 0 < η 1 [42,43]:
Definition 1.
(Hölder class): Let X R d be a compact space. The Hölder class of functions Σ d ( η , K ) , with Hölder parameters η and K, consists of functions g that satisfy
g : g ( z ) p x η ( z ) d K x z d η , x , z X ,
where p x k ( 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 η.
To explore the optimal choice of parameter α we require bounds on the bias and variance bounds, provided in Appendix A.5. To obtain such bounds, we will make several assumptions on the absolutely continuous densities f X , f Y , f X Y and support sets S X , S Y , S X Y :
(A.1)
Each of the densities belong to Σ d ( η , K ) with smoothness parameters η and Lipschitz constant K.
(A.2)
The volumes of the support sets are finite, i.e., 0 < V ( S X ) < , 0 < V ( S Y ) < .
(A.3)
All densities are bounded i.e., there exist two sets of constants C X L , C Y L , C X Y L and C X U , C Y U , C X Y U such that 0 < C X L f X C X U < , 0 < C Y L f Y C Y U < and 0 < C X Y L f X Y C X Y U < .
The following theorem on the bias follows under assumptions (A.1) and (A.3):
Theorem 4.
For given α ( 0 , 1 ) , β = 1 α , d 2 , and 0 < η 1 the bias of the R n , n : = R n , n ( Z n , Z ˜ n ) satisfies
| E R n , n n 2 α β f X Y ( x , y ) f X ( x ) f Y ( y ) α f X Y ( x , y ) + β f X ( x ) f Y ( y ) d x d y | O max n η 2 / ( d ( 1 + η ) ) , ( β n ) η / ( 1 + η ) , c d n 1 ,
where c d is the largest possible degree of any vertex of MST on Z n Z ˜ n . The explicit form of (13) is provided in Appendix A.5.
Please note that according to Theorem 13 in [44], the constant c d is lower bounded by Ω d 2 n ( 1 H ( γ ) ) , γ = 2 d and H ( γ ) is the binary entropy i.e.,
H ( γ ) = γ log γ ( 1 γ ) log ( 1 γ ) .
A proof of Theorem 4 is given in Appendix A.5. 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 Efron–Stein inequality [45]). It is similar to arguments in ([46], Appendix A.3), and is omitted. In Theorem 5 we assume that the densities f X , f Y , and f X Y are absolutely continuous and bounded (A.3).
Theorem 5.
Given α ( 0 , 1 ) , the variance of the estimator R n , n : = R n , n ( Z n , Z ˜ n ) is bounded by
V a r R n , n n ( 1 α ) c d n , α = n / n ,
where c d is a constant depending only on the dimension d.

2.4. Minimax Parameter α

Recall assumptions (A.1), (A.2), and (A.3) in Section 2.3. The constant α can be chosen to minimize the maximum the MSE converges rate where the maximum is taken over the space of Hölder smooth joint densities f X Y .
Throughout this subsection we use the following notations:
  • ϵ X Y : = f X Y / f X f Y ,
  • C ϵ L : = C X Y L / C X U C Y U and C ϵ U : = C X Y U / C X L C Y L ,
  • C n : = C X Y L n / 2 ,
  • α 0 L : = 2 C n and α 0 U : = min 1 4 , 1 + 1 / C n 4 + 2 C ϵ U , 1 n η / d 1 , where η is the smoothness parameter,
  • l n : = n η / ( d 2 ( 1 + η ) ) .
Now define G ˜ ϵ X Y , n α , β ( x , y ) by
( ϵ X Y ( x , y ) + 1 / ( β C n ) ) ( 1 + ϵ X Y ( x y ) + 1 / ( β C n ) ) ( α + β ϵ X Y ( x , y ) ) 2 , β = 1 α .
Consider the following optimization problem:
min α max ϵ X Y Δ ˜ ( α , ϵ X Y ) + c d ( 1 α ) n 1 subject to C ϵ L ϵ X Y C ϵ U , α 0 L α α 0 U ,
where
Δ ˜ ( α , ϵ X Y ) : = D ( n , l n , d , η ) + D ˜ ( n , l n , d ) C X Y U S X Y G ˜ ϵ X Y , n α , β ( x , y ) d x d y ,
and
D ( n , l n , d , η ) = c 2 l n d n 1 + c d 2 d n 1 + c l n d n η / d + c l n d n 1 / d + 2 c 1 l n d 1 n 1 / d 1 + c 3 l n d η ,
D ˜ ( n , l n , d ) = 2 + n 1 2 c i = 1 M l n l n d a i 1 + n 3 / 2 2 c 1 i = 1 M l n l n d / 2 b i a i 2 + n 1 i = 1 M 2 n 3 / 2 l n d / 2 b i a i 2 n a i l n d + n 2 a i 2 1 / 2 n b i l n d + n 2 b i 2 1 / 2 .
Please note that in (18), c , c , c 1 , c 2 are constants, and c d only depends on the dimension d. Also, in (19), a i and b i are constants. Let ϵ X Y * be the optimal ϵ X Y i.e., ϵ X Y * be the solution of the optimization problem (16). Set
Ξ ( α ) : = d d α Δ ˜ ( α , ϵ X Y * ) + c d ( 1 α ) n 1 ,
such that Δ ˜ ( α , ϵ X Y * ) is (17) when ϵ X Y = ϵ X Y * . For α [ α 0 L , α 0 U ] , the optimal choice of ϵ X Y in terms of maximizing the MSE is ϵ X Y * = C ϵ U and the saddle point for the parameter α , denoted by α ˜ , is given as follows:
  • α ˜ = α 0 U , if Ξ ( α 0 U ) < 0 .
  • α ˜ = α 0 L , if Ξ ( α 0 L ) > 0 .
  • α ˜ = Ξ 1 ( 0 ) , if α 0 L Ξ 1 ( 0 ) α 0 U .
Further details are given in Appendix A.6.

3. 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 GMI 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 density plug-in estimators of the GMI, in particular the standard KDE density plug-in estimator of [22], where the convergence rates of Theorems 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 the proposed GMI estimator approach with KDE method.
Here we perform four sets of experiments to illustrate the proposed GMI estimator. For the first set of experiments the MSE of the GMI estimator in Algorithm 1 is shown in Figure 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 ] . Please 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 Figure 2-left validates the theoretical MSE growth rates derived from (13) and (14), i.e., decreasing sub-linearly in n and increasing exponentially in d.
In Figure 2-right, we compare the proposed MST-based GMI estimator with the KDE-GMI estimator [22]. For the KDE approach, we estimated the joint and marginal densities and then plugged them into the proposed expression (2). The bandwidth h used for the KDE plug-in estimator was set as h = n 1 / ( d + 1 ) . The choice of h minimizes the bound on the MSE of the plug-in estimator. We generated data from the two-dimensional normal distribution with zero mean and covariance matrix
1 ρ ρ 1 .
The coefficient ρ is varied in range [ 0.1 , 0.9 ] . The true GMI was computed by the Monte Carlo approximation to the integral (2). Please note that as ρ increases, the MST-GMI outperforms the KDE-GMI approach. In this set of experiments α = 0.6 .
Figure 3 again compares the MST-GMI estimator with the KDE-GMI estimator. samples are drawn 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 Figure 3 show the MSE (100 trials) of the GMI estimator implemented with an KDE estimator (with bandwidth as in Figure 2 i.e., h = n 1 / ( d + 1 ) ) for dimensions d = 4 , 12 and various sample sizes. For all dimensions and sample sizes the MST-GMI estimator also outperforms the plug-in KDE-GMI estimator based on the estimated log–log MSE slope given in Figure 3 (left plots). The right plots in Figure 3 compares the MST-GMI with the KDE-GMI. In this experiment, the error bars denote standard deviations with 100 trials. We observe that for higher dimension d = 12 and larger sample size n, the KDE-GMI approaches the true GMI at a slower rate than the MST-GMI estimator. This reflects the power of the proposed graph-based approach to estimating GMI.
The comparison between MSEs for various dimension d is shown in Figure 4 (left). This experiment highlights the impact of higher dimension on the GMI estimators. 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 Figure 4 (left) we observe that for larger sample size, MSE curves are ordered based on their corresponding dimensions. Results in Section 2.4 strongly depend on the lower bounds C X L , C Y L , C X Y L and upper bounds C X U , C Y U , C X Y U and provide optimal parameter α in the range [ α 0 L , α 0 U ] , therefore in the experiment section we only analyze one case where the lower bounds C X L , C Y L , C X Y L and upper bounds C X U , C Y U , C X Y U are known and the optimal α becomes α 0 L . Figure 4 (right) illustrates the MSE vs proportion parameter α when n = 500 , 10 4 samples are generated from truncated normal distribution with ρ = 0.7 , 0.5 . First, following Section 2.4, we compute the bound [ α 0 L , α 0 U ] and then derive the optimal α in this range. Therefore, each experiment with different sample size and ρ provides different range [ α 0 L , α 0 U ] . We observe that the MSE does not appeared a monotonic function in α and its behavior strongly depends on sample size n, d, and density functions’ bounds. Additional study of the dependency is described in Appendix A.6. In this set of experiments Ξ ( α 0 L ) > 0 , therefore following the results in Section 2.4, we have α ˜ = α 0 L . In this experiment the optimal value of α is always the lower bound α 0 L and indicated in the Figure 4 (right).
The parameter α is studied further for three scenarios where the lower bounds C X L , C Y L , C X Y L and upper bounds C X U , C Y U , C X Y U are assumed unknown, therefore results in Section 2.4 are not applicable. In this set of experiments, we varied α in the range ( 0 , 1 ) to divide our original sample. We generated sample from an isotropic multivariate standard normal distribution ( ρ = 0 ) in all three scenarios (all features are independent). Therefore, the true GMI is zero and in all scenarios the GMI column, corresponding to the MST-GMI, 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 with α = 0.5 .
Finally, Figure 5 shows the runtime 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-GMI method is slightly faster but as n becomes large we see significant relative speedup of the proposed MST-GMI method.

4. Conclusions

In this paper, we have proposed a new measure of mutual information, called Geometric MI (GMI), which 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, which depends on the MST over all data points. We established some properties of the GMI in terms of convexity/concavity, chain rule, and a type of data-processing inequality. A direct estimator of the GMI, called the MST-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.

Author Contributions

S.Y.S. wrote this article primarily under the supervision of A.O.H. as principle investigator (PI), and A.O.H. edited the paper. S.Y.S. provided the primary contributions for the proofs of all theorems and performed all experiments.

Funding

The work presented in this paper was partially supported by ARO grant W911NF-15-1-0479 and DOE grant DE-NA0002534.

Acknowledgments

The authors would like to thank Brandon Oselio for the helpful comments.

Conflicts of Interest

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Appendix A

We organize the appendices as follows: Theorem 1 which establishes convexity/concavity is proved in Appendix A.1. Appendix A.2 and Appendix A.3 establish the inequality (9) and (10) for given p ( 0 , 1 ) , respectively. In Appendix A.4, we first prove that the set Z ˜ n which is randomly generated from original dependent data, contains asymptotically independent samples. Later by using the generated independent sample Z ˜ n we show that for given α the FR estimator of the GMI given in Algorithm 1 tends to I α . Appendix A.5 dedicates proof of Theorem 4. The proportionality parameter ( α ) optimization strategy is presented in Appendix A.6.

Appendix A.1. Theorem 1

Proof. 
The proof is similar to the result for standard (Shannon) mutual information. However, we require the following lemma, proven in analogous manner to the log-sum inequality:
Lemma A1.
For non-negative real numbers α 1 , , α n and β 1 , , β n , given p ( 0 , 1 ) , q = 1 p ,
i = 1 n α i p β i α i + q 1 i = 1 n α i p i = 1 n β i i = 1 n α i + q 1 .
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 x y for d x , d y and d x d y , respectively. To prove part (i) of Theorem 1, we represent the LHS of (5) as:
I ˜ p λ 1 f Y | X g X + λ 2 f Y | X h X = 1 x y λ 1 f Y | X g X + λ 2 f Y | X h X × p λ 1 f Y | X g X + λ 2 f Y | X h X x λ 1 f Y | X g X + λ 2 f Y | X h X y λ 1 f Y | X g X + λ 2 f Y | X h X + q 1 = 1 x y λ 1 f Y | X g X + λ 2 f Y | X h X p f Y | X x λ 1 f Y | X g X + λ 2 f Y | X h X + q 1 .
Furthermore, the RHS of (5) can be rewritten as
λ 1 I ˜ p ( f Y | X g X ) + λ 2 I ˜ p ( f Y | X h X ) = 1 x y λ 1 f Y | X g X p f Y | X g X x f Y | X g X y f Y | X g X + q 1 + λ 2 f Y | X h X p f Y | X h X x f Y | X h X y f Y | X h X + q 1 = 1 x y λ 1 f Y | X g X p f Y | X x f Y | X g X + q 1 + λ 2 f Y | X h X p f Y | X x f Y | X h X + q 1 .
Thus, to prove LHS RHS , we use the inequality below:
λ 1 f Y | X g X + λ 2 f Y | X h X p f Y | X x λ 1 f Y | X g X + λ 2 f Y | X h X + q 1 λ 1 f Y | X g X p π x f Y | X g X + q 1 + λ 2 f Y | X h X p π x f Y | X h X + q 1 .
In Lemma A1, let
α 1 = λ 1 x f Y | X g X λ 1 f Y | X g X + λ 2 f Y | X h X x λ 1 f Y | X g X + λ 2 f Y | X h X , α 2 = λ 2 x f Y | X h X λ 1 f Y | X g X + λ 2 f Y | X h X x λ 1 f Y | X g X + λ 2 f Y | X h X ,
and for i = 1 , 2 ,
β i = λ i f Y | X λ 1 f Y | X g X + λ 2 f Y | X h X x λ 1 f Y | X g X + λ 2 f Y | X h X .
Then the claimed assertion (i) is obtained. Part (ii) follows by convexity of D p and the following expression:
I ˜ p λ 1 g Y | X f X + λ 2 h Y | X f X = D p λ 1 f X g Y | X + λ 2 f X h Y | X , x λ 1 f X g Y | X + λ 2 f X h Y | X y λ 1 ϕ π 1 + λ 2 ϕ π 2 = D p ( λ 1 f X g Y | X + λ 2 f X h Y | X , f X x λ 1 f X g Y | X + λ 2 f X h Y | X = D p λ 1 f X g Y | X + λ 2 f X h Y | X , λ 1 x f X g Y | X y f X g Y | X + λ 2 x f X h Y | X y f X h Y | X .
Therefore, the claim in (6) is proved. □

Appendix A.2. Theorem 2

Proof. 
We start with (9). Given p ( 0 , 1 ) and q = 1 p , we can easily check that for positive t > q , s > q , such that t , s 1 :
( t + s ) ( t s ) + q ( 1 t s ) p ( t s ) 0 .
This implies
p t q p s q p + q t s t + s .
By substituting
f X 1 Y ( x 1 , y ) f X 1 ( x 1 ) f Y ( y ) = t q p , f X 2 Y | X 1 ( x 2 , y | x 1 ) f X 2 | X 1 ( x 2 | x 1 ) f Y | X 1 ( y | x 1 ) = s q p ,
we get
p f X 1 X 2 Y ( x 1 , x 2 , y ) f X 1 X 2 ( x 1 , x 2 ) f Y ( y ) + q 1 p f X 1 Y ( x 1 , y ) f X 1 ( x 1 ) f Y ( y ) + q 1 + p f X 2 Y | X 1 ( x 2 , y | x 1 ) f X 2 | X 1 ( x 2 | x 1 ) f Y | X 1 ( y | x 1 ) + q 1 .
Consequently
I p ( X 1 , X 2 ; Y ) I p ( X 1 ; Y ) E f p f X 2 Y | X 1 ( x 2 , y | x 1 ) f X 2 | X 1 ( x 2 | x 1 ) f Y | X 1 ( y | x 1 ) + q 1 .
Here f is the joint PDF of random vector ( X 1 , X 2 , Y ) . From the conditional GMI definition in (7) the expectation term in (A2) is equivalent to 1 I p ( X 2 ; Y | X 1 ) . This completes the proof. □

Appendix A.3. Proposition 1

Proof. 
Recall the Theorem 2, part (i). First from X Y Z we have f X Y Z = f X Y f Z | Y and then by applying the Jensen inequality,
I p ( X ; Y ) = I p ( X ; Y , Z ) and I p ( X ; Y , Z ) I p ( Z ; X ) E p π ( X , Y , Z ) + q 1 I p ( Z ; X ) p E π ( X , Y , Z ) + q 1 ,
where
π ( x , y , z ) = f Y X | Z ( y , x | z ) f Y | Z ( y | z ) f X | Z ( x | z ) .
Now by Markovian property we can immediately simplify the RHS in (A3) to the RHS in (10).
Furthermore, we can easily show that if X Z Y , we have f X Y Z = f Z X f Y | Z and therefore I p ( Z ; X ) = I p ( X ; Y , Z ) . This together with (A3) 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 A.4. Theorem 3

Proof. 
We first derive two required Lemmas A2 and A3 below:
Lemma A2.
Consider random vector Z = ( X , Y ) with joint probability density function (pdf) f X Y . Let Z n = { z 1 , , z n } = { ( x i , y i ) } i = 1 n be a set of samples with pdf f X Y . 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 mutually independent.
Proof. Consider two subsets A , B R n , then we have
P ( X ˜ A , Y ˜ B ) = E I A ( X ˜ ) . I B ( Y ˜ ) = E i , j I A ( X i ) . I B ( Y j ) . P ( X ˜ , Y ˜ ) = ( X i , Y j ) | Z n .
Here I A stands for the indicator function. Please note that
P ( X ˜ , Y ˜ ) = ( X i , Y j ) | Z n = 1 n 2 ,
and X i and Y j , i j are independent, therefore
P ( X ˜ A , Y ˜ B ) = 1 n 2 i j P ( X i A ) P ( Y j B ) + 1 n 2 i = 1 n P ( X i A , Y i B ) = P ( X i A ) P ( Y j B ) + 1 n P ( X i A , Y i B ) P ( X i A ) P ( Y i B ) ,
this implies that
| P ( X ˜ A , Y ˜ B ) P ( X ˜ A ) P ( Y ˜ B ) | 1 n | f X Y ( x , y ) f X ( x ) f Y ( y ) | d x d y .
On the other hand, we know that n = β n , so we get
| P ( X ˜ A , Y ˜ B ) P ( X ˜ A ) P ( Y ˜ B ) | 1 β n | f X Y ( x , y ) f X ( x ) f Y ( y ) | d x d y .
From (A5), 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 so that β = constant then the RHS in (A5) tends to zero. This implies that X ˜ and Y ˜ become independent when n . □
An immediate result of Lemma A2 is the following:
Lemma A3.
For given random vector Z n = ( X n , Y n ) from joint density function f X Y and with marginal density functions f X and f Y let Z ˜ n = z ˜ 1 , , z ˜ n be realization of random vector Z ˜ as in Lemma A2 with parameter β = n / n . Then for given points of Z ˜ n at z ˜ = ( x ˜ , y ˜ ) , we have
| f Z ˜ ( x ˜ , y ˜ ) f X ( x ˜ ) f Y ( y ˜ ) | = O 1 β n .
Now, we want to provide a proof of assertion (11). Consider two subsets Z n and Z ˜ n as described in Section 2.2. Assume that the components of sample Z ˜ n follow density function f ˜ X ˜ Y ˜ . Therefore by owing to Lemmas A2 and A3, 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 , and denote the FR statistic R n , n on these processes. Following the arguments in [13,46] we shall prove the following:
E R n , n n + n 2 α β f X , Y ( x , y ) f X ( x ) f Y ( y ) α f X Y ( x , y ) + β f X ( x ) f Y ( y ) d x d y .
This follows due to | R n , n R n , n | K d | M n n | + | N n n | , where K d is a constant defined in Lemma 1, [13] and n + n = n . Thus, ( n + n ) 1 E | R n , n R n , n | 0 as n . Let W 1 n , n , W 2 n , n , be independent variables with common density
ϕ n , n ( x , y ) = n f X Y ( x , y ) + n f ˜ X ˜ Y ˜ ( x , y ) / ( n + n ) ,
for ( x , y ) R d × R d . Let L n , n be an independent Poisson variable with mean n + n . Let F n , n = W 1 n , n , , W L n , n n , n a non-homogeneous Poisson process of rate n f X Y + n f ˜ X ˜ Y ˜ . Assign mark 1 to a point in F n , n with probability
n f X Y ( x , y ) / n f X Y ( x , y ) + n f ˜ X ˜ Y ˜ ( x , y ) ,
and mark 2 otherwise. By the marking theorem [13,47], the FR test statistic R ˜ n , n has the 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 (A7).
g n , n ( z , z ) = n f X Y ( x , y ) n f ˜ X ˜ , Y ˜ ( x , y ) + n f ˜ X ˜ , Y ˜ ( x , y ) n f X Y ( x , y ) n f X Y ( x , y ) + n f ˜ X ˜ , Y ˜ ( x , y ) n f X Y ( x , y ) + n f ˜ X ˜ , Y ˜ ( x , y ) ,
define
g ( z , z ) = α β f X Y ( x , y ) f X ( x ) f Y ( y ) + f X ( x ) f Y ( y ) f X Y ( x , y ) α f X Y ( x , y ) + β f X ( x ) f Y ( y ) α f X Y ( x , y ) + β f X ( x ) f Y ( y ) ,
then
E R ˜ n , n | F n , n = i < j L n , n g n ( W i n , n , W j n , n ) I F n , n ( W i n , n , W j n , n ) .
Now recall (A8). We observe that g n , n ( z , z ) g ( z , z ) . Going back to (A9), we can write
E R ˜ n , n = i < j L n , n g n , n ( W i n , n , W j n , n ) I F n , n ( W i n , n , W j n , n ) + o ( n + n ) .
For fixed n , n consider the collection:
F n , n = W 1 n , n , , W n , n n + n .
By the fact that E M n + N n ( n + n ) = o ( n + n ) , we have
E R ˜ n , n = i < j n + n g n , n ( W i n , n , W j n , n ) I F n , n ( W i n , W j n ) + o ( n + n ) .
Introduce
ϕ ( x , y ) = α f X Y ( x , y ) + β f X ( x ) f Y ( y ) .
Then ϕ n , n ( x , y ) ϕ ( x , y ) uniformly as n / n α and n / n β . Thus, using Proposition 1 in [13], we get
E R ˜ n , n n g ( z , z ) ϕ ( z ) d z = 2 α β f X Y ( x , y ) f X ( x ) f Y ( y ) α f X Y ( x , y ) + β f X ( x ) f Y ( y ) d x d y .
 □

Appendix A.5. Theorem 4

Proof. 
We begin by providing a family of bias rate bounds for the FR test statistic R n , n in terms of a parameter l. Assume f X Y , f X , and f Y are in Σ d ( η , K ) . Then by plugging the optimal l, we prove the bias rate bound given in (13).
Theorem A1.
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
| E R n , n n 2 α β f X Y ( x , y ) f X ( x ) f Y ( y ) α f X Y ( x , y ) + β f X ( x ) f Y ( y ) d x | O l d ( n ) η / d + O l d η + O l d β 1 n 1 + O c d n 1 ,
where 0 < η 1 is the Hölder smoothness parameter and c d is the largest possible degree of any vertex of MST.
Set
α i = α n a i l d 1 a i l d + α n 2 a i 2 , β i = β n b i l d 1 b i l d + β n 2 b i 2 .
and
A f , n β , α ( x , y ) = 2 f X Y ( x , y ) f X ( x ) f Y ( y ) + δ f / ( β n ) f X Y ( x , y ) α + ( f X ( x ) f Y ( y ) + δ f / ( β n ) ) β a i 2 l d α f X Y ( x , y ) + β f X ( x ) f Y ( y ) + δ f / ( β n ) 2 ,
where
δ f = | f X Y ( x , y ) f X ( x ) f Y ( y ) | d x d y ,
A more explicit form for the bound on the RHS is given below:
Δ ( α , f X Y , f X f Y ) : = c 2 l d ( n ) 1 + c d 2 d ( n ) 1 + O l d ( n ) η / d + O l d ( n ) 1 / 2 + O c d ( n ) 1 / 2 + 2 c 1 l d 1 ( n ) ( 1 / d ) 1 + δ f ( ( β n ) 1 ) 2 α β f X Y ( x , y ) α f X Y ( x , y ) + β f X ( x ) f Y ( y ) d x d y + ( n ) 1 i = 1 M 2 f X Y ( x , y ) f X ( x ) f Y ( y ) + δ f / ( β n ) ( α i β i ( α n a i l d f X Y 2 ( x , y ) + β n b i l d f X ( x ) f Y ( y ) + δ f / ( β n ) 2 ) ) 1 / 2 / α n a i f X Y ( x , y ) + β n b i f X ( x ) f Y ( y ) 2 d x d y + ( n ) 1 i = 1 M O ( l ) l d ( a i ) 1 2 f X Y ( x , y ) f X ( x ) f Y ( y ) + δ f / ( β n ) α f X Y ( x , y ) + β f X ( x ) f Y ( y ) d x d y + O l d η + ( n ) 3 / 2 i = 1 M O ( l ) l d / 2 b i A f , n β , α ( x , y ) d x d y .
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 { Z ˜ 1 , , Z ˜ N n } . Likewise Appendix A.4, set R n , n = R ( Z n , Z ˜ n ) . Applying Lemma 1, and (12) in [13], we can write
| R n , n R n , n | c d ( | M n n | + | N n n | ) .
Here c d denotes the largest possible degree of any vertex of the MST in R d . Following the arguments in [46], we have E | M n n | = O n 1 / 2 and E | N n n | = O n 1 / 2 . Hence
E R n , n n + n = E R n , n n + n + O c d ( n + n ) 1 / 2 .
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 i = 1 l d a i l d = i = 1 l d b i l d = 1 . Therefore, using the subadditivity property in Lemma 2.2, [46], we can write
E R n , n i = 1 M E E R n i , n i | n i , n i + 2 c 1 l d 1 ( n + n ) 1 / d ,
where M = l d , and η > 0 is the Hölder smoothness parameter. Furthermore, for given n i , n i , let W 1 n i , n i , W 2 n i , n i , be independent variables with common densities for ( x , y ) R d × R d :
g n i , n i ( x , y ) = n i f X Y ( x , y ) + n i f ˜ X ˜ Y ˜ ( x , y ) / ( n i + n i ) .
Denote L n i , n i be an independent Poisson variable with mean n i + n i and F n i , n i = W 1 n i , n i , , W L n i . n i n i , n i a non-homogeneous Poisson of rate n i f X Y + n i f ˜ X ˜ Y ˜ . Let F n i , n i be the non-Poisson point process W 1 n i , n i , W n i + n i n i , n i . Assign a mark from the set { 1 , 2 } to each point of F n i , n i . Let Z ˜ n i be the sets of points marked 1 with each probability n i f X Y ( x , y ) / n i f X Y ( x , y ) + n i f ˜ X ˜ Y ˜ ( x , y ) and let Z ˜ n i be the set points with mark 2. Please note that owing to the marking theorem [47], 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 as FR test statistic on nodes in Z ˜ n i Z ˜ n i , we have
E R n i , n i | n i , n i = E R ˜ n i , n i | n i , n i .
By the fact that E | M n + N n n n | = O ( ( n + n ) 1 / 2 ) , we have
E R ˜ n i , n i | n i , n i = E E R ˜ n i , n i | F n i , n i = E s < j < n i + n i P n i , n i ( W s n i , n i , W j n i , n i ) 1 ( W s n i , n i , W j n i , n i ) F n i , n i + O ( n i + n i ) 1 / 2 ) .
Here z = ( x , y ) , z = ( x , y ) , and P n i , n i ( z , z ) is given in below:
P n i , n i ( z , z ) : = P r mark z mark z , ( z , z ) F n i , n i = n i f X Y ( x , y ) n i f ˜ X ˜ Y ˜ ( x , y ) + n i f ˜ X ˜ Y ˜ ( x , y ) n i f X , Y ( x , y ) n i f X Y ( x , y ) + n i f ˜ X ˜ Y ˜ ( x , y ) n 1 f X Y ( x , y ) + n i f ˜ X ˜ Y ˜ ( x , y ) .
Next set
α i = n a i l d ( 1 a i l d ) + n 2 a i 2 , β i = n b i l d ( 1 b i l d ) + n 2 b i 2 .
By owing to the Lemma B.6 in [46] and applying analogous arguments, we can rewrite the expression in (A20):
E R n , n i = 1 M a i b i l d 2 n n f X Y ( x , y ) f X ˜ Y ˜ ( x , y ) n a i f X , Y ( x , y ) + n b i f X ˜ Y ˜ ( x , y ) d x d y + 2 c 1 l d 1 ( n + n ) 1 / d + i = 1 M 2 f X Y ( x , y ) f X ˜ Y ˜ ( x , y ) α i β i n a i l d f X Y 2 ( x , y ) + n b i l d f X ˜ Y ˜ 2 ( x , y ) 1 / 2 n a i f X Y ( x , y ) + n b i f X ˜ Y ˜ ( x , y ) 2 d x d y + i = 1 M E n i , n i ( n i + n i ) ς η ( l , n i , n i ) + O l d ( n + n ) 1 η / d + O l d ( n + n ) 1 / 2 ,
where
ς η ( l , n i , n i ) = O l n i + n i 2 l d n i + n i g n i , n i ( z ) P n i , n i ( z , z ) d z + O ( l d η ) .
Going back to Lemma A3, we know that
f X ˜ Y ˜ ( x , y ) = f X ( x ) f Y ( y ) + O 1 β n .
Therefore, the first term on the RHS of (A20) is less and equal to
i = 1 M a i b i l d 2 n n f X Y ( x , y ) f X ( x ) f Y ( y ) n a i f X , Y ( x , y ) + n b i f X ( x ) f Y ( y ) d x d y + δ f β n i = 1 M a i b i l d 2 n n f X Y ( x , y ) n a i f X Y ( x , y ) + n b i f X ( x ) f Y ( y ) d x d y ,
and the second term is less and equal to
i = 1 M 2 f X Y ( x , y ) f X ( x ) f Y ( y ) + δ f / ( β n ) ( α i β i ( n a i l d f X Y 2 ( x , y ) + n b i l d ( f X ( x ) f Y ( y ) + δ f / ( β n ) ) 2 ) ) 1 / 2 / n a i f X Y ( x , y ) + n b i f X ( x ) f Y ( y ) 2 d x d y ,
where
δ f = | f X Y ( x , y ) f X ( x ) f Y ( y ) | d x d y .
Recall the definition of the dual MST and FR statistic denoted by R n , n * following [46]:
Definition A1.
(Dual MST, MST * and dual FR statistic R m , n * ) Let F i be the set of corner points associated with a particular subsection Q i , 1 i l d of [ 0 , 1 ] d . Define the dual MST * ( X m Y n Q i ) as the boundary MST graph in partition Q i [48], which contains X m and Y n points falling inside partition cell Q i and those corner points in F i which minimize total MST length. Please note that it is allowed to connect the MSTs in Q i and Q j through points strictly contained in Q i and Q j and therefore corner points. Thus, 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 by passing through the corner the corner points which are all connected in order to minimize the total weights. To clarify, 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 closer to a corner point otherwise the dual MST connects them to each other.
Furthermore, R m , n * ( X m , Y n Q i ) is defined as the number of edges in the MST * graph connecting nodes from different samples and number of edges connecting to the corner points. Please note that the edges connected to the corner nodes (regardless of the type of points) are always counted in the dual FR test statistic R m , n * .
Similarly, consider the Poisson processes samples and the FR test statistic over these samples, denoted by R n , n * . By superadditivity of the dual R n , n * [46], we have
E R n , n * i = 1 M a i l d 2 n n f X Y ( x , y ) f X ( x ) f Y ( y ) δ f / ( β n ) n f X Y ( x , y ) + n f X ( x ) f Y ( y ) δ f / ( β n ) d x d y i = 1 M E n i , n i ( n i + n i ) ς η ( l , n i , n i ) O l d ( n + n ) 1 η / d O l d ( n + n ) 1 / 2 c 2 l d .
The first term of RHS in (A21) is greater or equal to
2 n n f X Y ( x , y ) f X ( x ) f Y ( y ) n f X Y ( x , y ) + n f X ( x ) f Y ( y ) d x d y δ f β n 2 n n f X Y ( x , y ) n f X Y ( x , y ) + n f X ( x ) f Y ( y ) d x d y .
Furthermore,
E R n , n n + c d 2 d n E R n , n * n ,
where c d is the largest possible degree of any vertex of the MST in R d , as before. Consequently, we have
| E R n , n n 2 α β f X Y ( x , y ) f X ( x ) f Y ( y ) α f X Y ( x , y ) + β f X ( x ) f Y ( y ) d x d y | B ( α , f X Y , f X f Y ) ,
where B is defined in (A16) and A f , n β , α ( x , y ) has been introduced in (A14). The last line in (A22) follows from the fact that
i = 1 M E n i , n i ( n i + n i ) ς η ( l , n i , n i ) i = 1 M O ( l ) l d / 2 b i A f , n β , n / n ( x , y ) d x d y + i = 1 M O ( l ) l d ( a i ) 1 2 f X Y ( x , y ) f X ( x ) f Y ( y ) + O ( δ f / ( β n ) ) n f X Y ( x , y ) + n f X ( x ) f Y ( y ) d x d y .
Here A f , n β , n / n ( x , y ) is given as (A14) by substituting n / n in α such that β = 1 α . Hence, the proof of Theorem A1 is completed. □
Going back to the proof of (13), without loss of generality assume that ( n ) l d > 1 , for 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:
l ( n , β ) = a r g min l max l d ( n ) η / d , l η d , l d β 1 n 1 , c d 2 d n 1 .
The solution l = l ( n , β ) is obtained when l d ( n ) η / d = l η d , or equivalently l = ( n ) η / ( d 2 ( η + 1 ) ) or when l d β 1 n 1 = l η d , which implies l = ( β n ) 1 / ( d ( 1 + η ) ) . Substitute this l in the bound (A13) to obtain the RHS expression in (13) for d 2 . □

Appendix A.6.

Our main goal in Section 2.4 was to find proportion α such that the parametric MSE rate depending on the joint density f X Y and marginal densities f X , f Y is minimized. Recalling the explicit bias bound in (A16), it can be seen that this function will be a complicated function of f X Y , f X f Y and α . By rearrangement of terms in (A13), we first find an upper bound for Δ in (A16), denoted by Δ ¯ , as follows:
Δ ¯ ( α , f X Y , f X f Y ) = D ( n , l n , d , η ) + D ˜ ( n , l n , d ) E X Y G f , n α , β ( X , Y ) ,
where l n : = n η / ( d 2 ( 1 + η ) ) . From Appendix A.5 we know that optimal l is given by (A23). One can check that for α 1 n ( η / d ) 1 , the optimal l = n η / ( d 2 ( 1 + η ) ) provides a tighter bound. In (A24), the constants D and D ¯ are
D ( n , l n , d , η ) = c 2 l n d n 1 + c d 2 d n 1 + c l n d n η / d + c l n d n 1 / d + 2 c 1 l n d 1 n 1 / d 1 + c 3 l n d η ,
D ˜ ( n , l n , d ) = 2 + n 1 2 c i = 1 M l n l n d a i 1 + n 3 / 2 2 c 1 i = 1 M l n l n d / 2 b i a i 2 + n 1 i = 1 M 2 n 3 / 2 l n d / 2 b i a i 2 n a i l n d + n 2 a i 2 1 / 2 n b i l n d + n 2 b i 2 1 / 2 .
And the function G f , n α , β ( x , y ) is given as the following:
G f , n α , β ( x , y ) = f X ( x ) f Y ( y ) + δ f / ( n β ) ( α f X Y ( x , y ) + β f X ( x ) f Y ( y ) + δ f / ( β n ) / α f X Y ( x , y ) + β f X ( x ) f Y ( y ) 2 ,
where δ f is given in (A15). Next After all still the expression (A27) is complicated to optimize therefore we use the fact that 0 α , β 1 to bound the function G f , n α , β ( x , y ) . Define the set Γ
Γ : = ϵ X Y : | ϵ X Y ( t ) ϵ X Y ( t ) | K ¯ t t d η ,
where
K ¯ = C ϵ U K C X Y L + C X L + C Y L C X L C X U .
Here K is the smoothness constant in the Hölder class. Notice that set Γ is a convex set. We bound Δ ¯ by
Δ ˜ ( α , ϵ X Y ) = D ( n , l n , d , η ) + D ˜ ( n , l n , d ) C X Y U S X Y G ˜ ϵ X Y , n α , β ( x , y ) d x d y .
Set C n = C X Y L n / 2 ,
G ˜ n α , β ( ϵ X Y ) = ϵ X Y 1 ( x , y ) + ( β C n ) 1 ) ( 1 + ϵ X Y 1 ( x y ) + ( β C n ) 1 α + β ϵ X Y 1 ( x , y ) 2 .
This simplifies to
G ˜ n α , β ( ϵ X Y ) = 1 + ( β C n ) 1 ϵ X Y 1 + ϵ X Y + ( β C n ) 1 ϵ X Y α ϵ X Y + β 2 .
Under the condition
2 C n α min 1 2 + 1 2 C n , 1 3 + 2 3 C n ,
G ˜ n α , β ( ϵ X Y ) is an increasing function in ϵ . Furthermore, for α 1 4 and
C ϵ L ϵ X Y min C ϵ U , θ U ( α ) , where θ U ( α ) = 1 4 α + 1 / C n 2 α ,
the function G ˜ n α , β ( ϵ X Y ) is strictly concave. Next, to find an optimal α we consider the following optimization problem:
min α max ϵ X Y Γ Δ ˜ ( α , ϵ X Y ) + c d ( 1 α ) n 1 subject to C ϵ L ϵ X Y C ϵ U ,
here ϵ X Y = f X Y / f X f Y , C ϵ U = C X Y U / C X L C Y L and C ϵ L = C X Y L / C X U C Y U , such that C ϵ L 1 . We know that under conditions (A31) and (A32), the function G ˜ n α , β is strictly concave and increasing in ϵ X Y . We first solve the optimization problem:
max ϵ X Y Γ S X Y G ˜ n α , β ϵ X Y ( x , y ) d x d y subject to θ ϵ L ( α ) V ( S X Y ) S X Y ϵ X Y ( x , y ) d x d y θ ϵ U ( α ) V ( S X Y ) ,
where
θ ϵ L ( α ) : = C ϵ L , θ ϵ U ( α ) : = min { C ϵ U , θ U ( α ) } .
The Lagrangian for this problem is
L ( ϵ X Y , λ 1 , λ 2 ) = S X Y G ˜ n α , β ϵ X Y ( x , y ) d x d y λ 1 S X Y ϵ X Y ( x , y ) d x d y θ ϵ U ( α ) V ( S X Y ) λ 2 θ ϵ L ( α ) V ( S X Y ) S X Y ϵ X Y ( x , y ) d x d y .
In this case, the optimum ϵ X Y * is bounded, θ ϵ L ( α ) ϵ X Y * θ ϵ U ( α ) , and the Lagrangian multiplier λ 1 * , λ 2 * 0 is such that
min λ 1 , λ 2 0 max ϵ X Y Γ L ( ϵ X Y , λ 1 , λ 2 ) = L ( ϵ X Y * , λ 1 * , λ 2 * ) .
Set G n ( ϵ X Y ) = d d ϵ X Y G ˜ n α , β ( ϵ X Y ) . In view of the concavity of G ˜ ϵ X Y , n α , β and Lemma 1, page 227 in [49], maximizing L ( ϵ X Y , λ 1 * , λ 2 * ) over ϵ X Y is equivalent to
S X Y G n ϵ X Y * ( x , y ) ( λ 1 * λ 2 * ) ϵ X Y ( x , y ) d x d y 0 ,
for all θ ϵ L ( α ) ϵ X Y * θ ϵ U ( α ) , and
S X Y G n ϵ X Y * ( x , y ) ( λ 1 * λ 2 * ) ϵ X Y * ( x , y ) d x d y = 0 .
Denote G n 1 the inverse function of G n . Since G n is strictly decreasing in ϵ X Y * (this is because G ˜ n α , β ( ϵ X Y ) is strictly concave, so that G n 1 is continuous and strictly decreasing in ϵ X Y * ). From (A36) and (A37), we see immediately that on any interval θ ϵ L ( α ) ϵ X Y * θ ϵ U ( α ) , we have ϵ X Y * = G n 1 ( λ 1 * λ 2 * ) . We can write then
G n θ ϵ U ( α ) λ 1 * λ 2 * G n θ ϵ L ( α ) ,
and λ 1 * , λ 2 * 0 . Next, we find the solution of
min λ 1 , λ 2 0 G ¯ n α , β ( λ 1 , λ 2 ) , where
G ¯ n α , β ( λ 1 , λ 2 ) = V ( S X Y ) G ˜ n α , β G n 1 ( λ 1 λ 2 ) ( λ 1 λ 2 ) G n 1 ( λ 1 λ 2 ) + λ 1 θ ϵ U ( α ) λ 2 θ ϵ L ( α ) .
The function G ¯ n α , β ( λ 1 , λ 2 ) is increasing in λ 1 and λ 2 , and therefore it takes its minimum at ( λ 1 * , λ 2 * ) = ( G n θ ϵ U ( α ) , 0 ) . This implies that ϵ X Y * = θ ϵ U ( α ) . Returning to our primary minimization over α :
min α Δ ˜ ( α , ϵ X Y * ) + c d ( 1 α ) n 1 subject to α 0 L α α 0 U ,
where α 0 L = 2 C n and α 0 U = min 1 4 , 1 n η / d 1 . We know that 1 4 1 3 + 2 3 C n and 1 4 1 2 + 1 2 C n , therefore the condition below
2 C n α min 1 4 , 1 n η / d 1 ,
implies the constraint α 0 L α α 0 U . Since the objective function (A38) is a complicated function in α , it is not feasible to determine whether it is a convex function in α . For this reason, let us solve the optimization problem in (A38) in a special case when C ϵ U θ U ( α ) . This implies ϵ X Y * = C ϵ U . Under assumption C ϵ U the objective function in (A38) is convex in α . Also, the case C ϵ U θ U ( α ) is equivalent to α 1 + 1 / C n 4 + 2 C ϵ U . Therefore, in the optimization problem we have constraint
2 C n α min 1 4 , 1 + 1 / C n 4 + 2 C ϵ U , 1 n η / d 1 .
We know that Δ ˜ ( α , ϵ X Y * ) + c d ( 1 α ) n 1 is convex over α [ α 0 L , α 0 U ] . Therefore, the problem becomes ordinary convex optimization problem. Let α ˜ , λ ˜ 1 and λ ˜ 2 be any points that satisfy the KKT conditions for this problem:
α 0 L α ˜ 0 , α ˜ α 0 U 0 , λ ˜ 1 , λ ˜ 2 0 , λ ˜ 1 ( α 0 L α ˜ ) = 0 , λ ˜ 2 ( α ˜ α 0 U ) = 0 , d d α Δ ˜ ( α ˜ , ϵ X Y * ) + c d ( 1 α ˜ ) n 1 λ ˜ 1 + λ ˜ 2 = 0 .
Recall Ξ ( α ) from (20):
Ξ ( α ) = d d α Δ ˜ ( α , ϵ X Y * ) + c d ( 1 α ) n 1 ,
where Δ ˜ is given in (A28). So, the last condition in (A39) becomes Ξ ( α ˜ ) = λ ˜ 1 λ ˜ 2 . We then have
α 0 L Ξ 1 ( λ ˜ 1 λ ˜ 2 ) α 0 U ,
where Ξ 1 is inverse function of Ξ . Since α 0 L α 0 U , at least one of λ ˜ 1 or λ ˜ 2 should be zero:
  • λ ˜ 1 = 0 , λ ˜ 2 0 . Then α ˜ = α 0 U and implies λ ˜ 2 = Ξ ( α 0 U ) . Since λ ˜ 2 > 0 , so this leads to Ξ ( α 0 U ) < 0 .
  • λ ˜ 2 = 0 , λ ˜ 1 0 . Then α ˜ = α 0 L and implies λ ˜ 1 = Ξ ( α 0 L ) . We know that λ ˜ 1 > 0 , hence Ξ ( α 0 L ) > 0 .
  • λ ˜ 1 = 0 , λ ˜ 2 = 0 . Then α ˜ = Ξ 1 ( 0 ) and so α 0 L Ξ 1 ( 0 ) α 0 U .
Consequently, by following the behavior of Ξ ( α ) with respect to α 0 L and α 0 U , we can often find the optimal α ˜ , λ ˜ 1 and λ ˜ 2 . For instance, if Ξ ( α ) is positive for all α [ α 0 L , α 0 U ] then we conclude that α ˜ = α 0 L .

References

  1. Lewi, J.; Butera, R.; Paninski, L. Real-time adaptive information theoretic optimization of neurophysiology experiments. In Advances in Neural Information Processing Systems; The MIT Press: Cambridge, MA, USA, 2006; pp. 857–864. [Google Scholar]
  2. Peng, H.C.; Herskovits, E.H.; Davatzikos, C. Bayesian Clustering Methods for Morphological Analysis of MR Images. In Proceedings of the IEEE International Symposium on Biomedical Imaging, Washington, DC, USA, 7–10 July 2002; pp. 485–488. [Google Scholar]
  3. Moon, K.R.; Noshad, M.; Yasaei Sekeh, S.; Hero, A.O. Information theoretic structure learning with confidence. In Proceedings of the 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), New Orleans, LA, USA, 5–9 March 2017. [Google Scholar]
  4. Brillinger, D.R. Some data analyses using mutual information. Braz. J. Probab. Stat. 2004, 18, 163–183. [Google Scholar]
  5. Torkkola, K. Feature extraction by non parametric mutual information maximization. J. Mach. Learn. Res. 2003, 3, 1415–1438. [Google Scholar]
  6. Vergara, J.R.; Estévez, P.A. A review of feature selection methods based on mutual information. Neural Comput. Appl. 2014, 24, 175–186. [Google Scholar] [CrossRef]
  7. Peng, H.; Long, F.; Ding, C. Feature selection based on mutual information criteria of max-dependency, max-relevance. IEEE Trans. Pattern Anal. Mach. Intell. 2005, 27, 1226–1238. [Google Scholar] [CrossRef] [PubMed]
  8. Peng, H.; Long, F.; Ding, C. Evaluation, Application, and Small Sample Performance. IEEE Trans. Pattern Anal. Mach. Intell. 1997, 19, 153–158. [Google Scholar]
  9. Sorjamaa, A.; Hao, J.; Lendasse, A. Mutual information and knearest neighbors approximator for time series prediction. Lecture Notes Comput. Sci. 2005, 3697, 553–558. [Google Scholar]
  10. Mohamed, S.; Rezende, D.J. Variational information maximization for intrinsically motivated reinforcement learning. In Advances in Neural Information Processing Systems; The MIT Press: Cambridge, MA, USA, 2015; pp. 2116–2124. [Google Scholar]
  11. Neemuchwala, H.; Hero, A.O. Entropic graphs for registration. In Multi-Sensor Image Fusion and its Applications; CRC Press Book: Boca Raton, FL, USA, 2005; pp. 185–235. [Google Scholar]
  12. Neemuchwala, H.; Hero, A.O.; Zabuawala, S.; Carson, P. Image registration methods in high-dimensional space. Int. J. Imaging Syst. Technol. 2006, 16, 130–145. [Google Scholar] [CrossRef] [Green Version]
  13. Henze, N.; Penrose, M.D. On the multivarite runs test. Ann. Stat. 1999, 27, 290–298. [Google Scholar]
  14. Berisha, V.; Hero, A.O. Empirical non-parametric estimation of the Fisher information. IEEE Signal Process. Lett. 2015, 22, 988–992. [Google Scholar] [CrossRef]
  15. Berisha, V.; Wisler, A.; Hero, A.O.; Spanias, A. Empirically estimable classification bounds based on a nonparametric divergence measure. IEEE Trans. Signal Process. 2016, 64, 580–591. [Google Scholar] [CrossRef]
  16. Kailath, T. The divergence and Bhattacharyya distance measures in signal selection. IEEE Trans. Commun. Technol. 1967, 15, 52–60. [Google Scholar] [CrossRef]
  17. Yasaei Sekeh, S.; Oselio, B.; Hero, A.O. Multi-class Bayes error estimation with a global minimal spanning tree. In Proceedings of the 2018 56th Annual Allerton Conference on Communication, Control, and Computing (Allerton), Monticello, IL, USA, 2–5 October 2018. [Google Scholar]
  18. Noshad, M.; Moon, K.R.; Yasaei Sekeh, S.; Hero, A.O. Direct Estimation of Information Divergence Using Nearest Neighbor Ratios. In Proceedings of the IEEE International Symposium on Information Theory (ISIT), Aachen, Germany, 25–30 June 2017. [Google Scholar]
  19. March, W.; Ram, P.; Gray, A. Fast Euclidean minimum spanning tree: Algorithm, analysis, and applications. In Proceedings of the 16th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, Boston, MA, USA, 20–23 August 2010; pp. 603–612. [Google Scholar]
  20. Borůvka, O. O jistém problému minimálním. Práce Moravské Pridovedecké Spolecnosti 1926, 3, 37–58. [Google Scholar]
  21. Kraskov, A.; Stögbauer, H.; Grassberger, P. Estimating mutual information. Phys. Rev. E 2004, 69, 066–138. [Google Scholar] [CrossRef] [PubMed]
  22. Moon, K.R.; Sricharan, K.; Hero, A.O. Ensemble Estimation of Mutual Information. In Proceedings of the IEEE International Symposium on Information Theory (ISIT), Aachen, Germany, 25–30 June 2017; pp. 3030–3034. [Google Scholar]
  23. Moon, K.R.; Hero, A.O. Multivariate f-divergence estimation with confidence. In Proceedings of the Advances in Neural Information Processing Systems 27 (NIPS 2014), Montreal, QC, Canada, 8–13 December 2014; pp. 2420–2428. [Google Scholar]
  24. Moon, K.R.; Sricharan, K.; Greenewald, K.; Hero, A.O. Improving convergence of divergence functional ensemble estimators. In Proceedings of the IEEE International Symposium on Information Theory (ISIT), Barcelona, Spain, 10–15 July 2016. [Google Scholar]
  25. Leonenko, N.; Pronzato, L.; Savani, V. A class of Rényi information estimators for multidimensional densities. Ann. Stat. 2008, 36, 2153–2182. [Google Scholar] [CrossRef]
  26. Gao, S.; Ver Steeg, G.; Galstyan, A. Efficient estimation of mutual information for strongly dependent variables. In Proceedings of the Eighteenth International Conference on Artificial Intelligence and Statistics, San Diego, CA, USA, 9–12 May 2015; pp. 277–286. [Google Scholar]
  27. Pál, D.; Póczos, B.; Szapesvári, C. Estimation of Rényi entropy and mutual information based on generalized nearest-neighbor graphs. In Proceedings of the 24th Annual Conference on Neural Information Processing Systems 2010, Vancouver, BC, Canada, 6–9 December 2010. [Google Scholar]
  28. Krishnamurthy, A.; Kandasamy, K.; Póczos, B.; Wasserman, L. Nonparametric estimation of Rényi divergence and friends. In Proceedings of the 31st International Conference on Machine Learning, Beijing, China, 21–26 June 2014; pp. 919–927. [Google Scholar]
  29. Kandasamy, K.; Krishnamurthy, A.; Póczos, B.; Wasserman, L.; Robins, J. Nonparametric von mises estimators for entropies, divergences and mutual informations. In Proceedings of the Advances in Neural Information Processing Systems, Montreal, QC, Canada, 7–12 December 2015; pp. 397–405. [Google Scholar]
  30. Singh, S.; Póczos, B. Analysis of k nearest neighbor distances with application to entropy estimation. In Proceedings of the Advances in Neural Information Processing Systems, Barcelona, Spain, 5–10 December 2016; pp. 1217–1225. [Google Scholar]
  31. Sugiyama, M. Machine learning with squared-loss mutual information. Entropy 2012, 15, 80–112. [Google Scholar] [CrossRef]
  32. Costa, A.; Hero, A.O. Geodesic entropic graphs for dimension and entropy estimation in manifold learning. IEEE Trans. Signal Process. 2004, 52, 2210–2221. [Google Scholar] [CrossRef]
  33. Friedman, J.H.; Rafsky, L.C. Multivariate generalizations of the Wald-Wolfowitz and Smirnov two-sample tests. Ann. Stat. 1979, 7, 697–717. [Google Scholar] [CrossRef]
  34. Smirnov, N.V. On the estimation of the discrepancy between empirical curves of distribution for two independent samples. Bull. Moscow Univ. 1939, 2, 3–6. [Google Scholar]
  35. Wald, A.; Wolfowitz, J. On a test whether two samples are from the same population. Ann. Math. Stat. 1940, 11, 147–162. [Google Scholar] [CrossRef]
  36. Noshad, M.; Hero, A.O. Scalable mutual information estimation using dependence graphs. In Proceedings of the International Conference on Artificial Intelligence and Statistics (AISTAT), Brighton, UK, 12–17 May 2019. [Google Scholar]
  37. Yasaei Sekeh, S.; Oselio, B.; Hero, A.O. A Dimension-Independent discriminant between distributions. In Proceedings of the 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Calgary, AB, Canada, 15–20 April 2018. [Google Scholar]
  38. Csiszár, I. Information-type measures of difference of probability distributions and indirect observations. Studia Sci. Math. Hungar. 1967, 2, 299–318. [Google Scholar]
  39. Ali, S.; Silvey, S.D. A general class of coefficients of divergence of one distribution from another. J. R. Stat. Soc. Ser. B (Methodol.) 1966, 28, 131–142. [Google Scholar] [CrossRef]
  40. Cover, T.; Thomas, J. Elements of Information Theory, 1st ed.; John Wiley & Sons: Chichester, UK, 1991. [Google Scholar]
  41. Härdle, W. Applied Nonparametric Regression; Cambridge University Press: Cambridge, UK, 1991. [Google Scholar]
  42. Lorentz, G.G. Approximation of Functions; Holt, Rinehart and Winston: New York, NY, USA; Chicago, IL, USA; Toronto, ON, Canada, 1966. [Google Scholar]
  43. Andersson, P. Characterization of pointwise Hölder regularity. Appl. Comput. Harmon. Anal. 1997, 4, 429–443. [Google Scholar] [CrossRef]
  44. Robins, G.; Salowe, J.S. On the maximum degree of minimum spanning trees. In Proceedings of the SCG 94 Tenth Annual Symposium on Computational Geometry, Stony Brook, NY, USA, 6–8 June 1994; pp. 250–258. [Google Scholar]
  45. Efron, B.; Stein, C. The jackknife estimate of variance. Ann. Stat. 1981, 9, 586–596. [Google Scholar] [CrossRef]
  46. Yasaei Sekeh, S.; Noshad, M.; Moon, K.R.; Hero, A.O. Convergence Rates for Empirical Estimation of Binary Classification Bounds. arXiv 2018, arXiv:1810.01015. [Google Scholar]
  47. Kingman, J.F.C. Poisson Processes; Oxford Univ. Press: Oxford, UK, 1993. [Google Scholar]
  48. Yukich, J.E. Probability Theory of Classical Euclidean Optimization; Vol. 1675 of Lecture Notes in Mathematics; Springer: Berlin, Germany, 1998. [Google Scholar]
  49. Luenberger, D.G. Optimization by Vector Space Methods; Wiley Professional Paperback Series; Wiley-Interscience: Hoboken, NJ, USA, 1969. [Google Scholar]
Figure 1. 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 .
Figure 1. 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 .
Entropy 21 00787 g001
Figure 2. (left) Log–log plot of theoretical and experimental MSE of the proposed MST-based GMI estimator as a function of sample size n for d = 6 , 10 , 12 and fixed smoothness parameter η . (right) The GMI estimator was implemented using two approaches, Algorithm 1 and KDE method where the KDE-GMI used KDE density estimators in the formula (2). In this experiment, samples are generated from the two-dimensional normal distribution with zero mean and covariance matrix (21) for various value of ρ [ 0.1 , 0.9 ] .
Figure 2. (left) Log–log plot of theoretical and experimental MSE of the proposed MST-based GMI estimator as a function of sample size n for d = 6 , 10 , 12 and fixed smoothness parameter η . (right) The GMI estimator was implemented using two approaches, Algorithm 1 and KDE method where the KDE-GMI used KDE density estimators in the formula (2). In this experiment, samples are generated from the two-dimensional normal distribution with zero mean and covariance matrix (21) for various value of ρ [ 0.1 , 0.9 ] .
Entropy 21 00787 g002
Figure 3. MSE log–log plots as a function of sample size n (left) for the proposed MST-GMI estimator (“Estimated GMI”) and the standard KDE-GMI plug-in estimator of GMI. The right column of plots correspond to the GMI estimated for dimension d = 4 (top) and d = 12 (bottom). In both cases the proportionality parameter α is 0.5 . The MST-GMI estimator in both plots for sample size n in [ 700 , 1600 ] outperforms the KDE-GMI estimator, especially for larger dimensions.
Figure 3. MSE log–log plots as a function of sample size n (left) for the proposed MST-GMI estimator (“Estimated GMI”) and the standard KDE-GMI plug-in estimator of GMI. The right column of plots correspond to the GMI estimated for dimension d = 4 (top) and d = 12 (bottom). In both cases the proportionality parameter α is 0.5 . The MST-GMI estimator in both plots for sample size n in [ 700 , 1600 ] outperforms the KDE-GMI estimator, especially for larger dimensions.
Entropy 21 00787 g003
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 α [ α 0 L , α 0 U ] .
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 α [ α 0 L , α 0 U ] .
Entropy 21 00787 g004
Figure 5. Runtime of KDE approach and proposed MST-based estimator of GMI vs sample size. The proposed GMI estimator achieves significant speedup, while for small sample size, the KDE method becomes overly fast. Please note that in this experiment the sample is generated from the Gaussian distribution in dimension d = 2 .
Figure 5. Runtime of KDE approach and proposed MST-based estimator of GMI vs sample size. The proposed GMI estimator achieves significant speedup, while for small sample size, the KDE method becomes overly fast. Please note that in this experiment the sample is generated from the Gaussian distribution in dimension d = 2 .
Entropy 21 00787 g005
Table 1. Comparison between different scenarios of various dimensions and sample sizes in terms of parameter α . We applied the MST-GMI estimator 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 MST-GMI estimator provides lowest MSE when α = 0.5 indicated by star (*).
Table 1. Comparison between different scenarios of various dimensions and sample sizes in terms of parameter α . We applied the MST-GMI estimator 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 MST-GMI estimator provides lowest MSE when α = 0.5 indicated by star (*).
Overview Table for Different d, n, and α
ExperimentsDimension ( d )Sample Size ( n ) GMI ( I α )MSE ( × 10 4 )Parameter ( α )
Scenario 1–1610000.0229120.2
Scenario 1–2610000.01434.79440.5 *
Scenario 1–3610000.01766.38670.8
Scenario 2–1815000.0246110.2
Scenario 2–2815000.00741.60530.5 *
Scenario 2–3815000.01375.38630.8
Scenario 3–11020000.00742.36040.2
Scenario 3–21020000.00290.541800.5 *
Scenario 3–31020000.0262110.8

Share and Cite

MDPI and ACS Style

Yasaei Sekeh, S.; Hero, A.O. Geometric Estimation of Multivariate Dependency. Entropy 2019, 21, 787. https://doi.org/10.3390/e21080787

AMA Style

Yasaei Sekeh S, Hero AO. Geometric Estimation of Multivariate Dependency. Entropy. 2019; 21(8):787. https://doi.org/10.3390/e21080787

Chicago/Turabian Style

Yasaei Sekeh, Salimeh, and Alfred O. Hero. 2019. "Geometric Estimation of Multivariate Dependency" Entropy 21, no. 8: 787. https://doi.org/10.3390/e21080787

APA Style

Yasaei Sekeh, S., & Hero, A. O. (2019). Geometric Estimation of Multivariate Dependency. Entropy, 21(8), 787. https://doi.org/10.3390/e21080787

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop