\section{Introduction}

Information theoretic quantities, such as mutual information and its generalizations, play an important role in various settings in machine learning and statistical estimation and inference. Here we summarize briefly the role of some generalizations of mutual information in learning (cf. Sec.~\ref{sec:math-inf-quantities} for a mathematical definition of these quantities).

\begin{enumerate}[leftmargin=*]
\item {\bf Conditional mutual information} measures the amount of information between two variables $X$ and $Y$ given a third variable $Z$ and is zero iff $X$ is independent of $Y$ given $Z$. CMI finds a wide range of applications in learning including causality testing \cite{dawid1979conditional, zhang2012kernel},   structure inference in graphical models \cite{whittaker2009graphical}, feature selection \cite{fleuret2004fast} as well as in providing privacy guarantees \cite{cuff2016differential}.
\item {\bf Total correlation} measures the degree to which a set of $N$ variables are independent of each other, and appears as a natural metric of interest in several machine learning problems, for example, in independent component analysis, the objective is to maximize the independence of the variables quantified through total correlation \cite{hyvarinen2000independent}. In feature selection, ensuring the independence of selected features is one goal, pursued using total correlation in \cite{vergara2014review,meyer2008information}.
\item {\bf Multivariate mutual information} measures the amount of information shared between multiple variables \cite{mcgill1954multivariate,chan2015multivariate} and is useful in feature selection \cite{lee2013feature,brown2009new} and clustering \cite{chan2016info}.
\vspace{-0.05in}\item {\bf Directed information} measures the amount of information between two random processes \cite{watanabe1960information,permuter2011interpretations} and is shown as the correct metric in identifying time-series graphical models \cite{quinn2015directed,sun2015causal,hlavavckova2007causality, amblard2011directed, rahimzamani2016network, rahimzamani2017potential}.
\end{enumerate}

Estimation of these information-theoretic quantities from observed samples is a non-trivial problem that needs to be solved in order to utilize these quantities in the aforementioned applications. While there has been long history in estimation of entropy \cite{kozachenko1987sample, beirlant1997nonparametric, wieczorkowski1999entropy, miller2003new}, and renewed recent interest \cite{lee2010sample,sricharan2013ensemble,han2015adaptive}, much less effort has been spent on the multivariate versions. A standard approach to estimating general information theoretic quantities is to write them out as a sum or difference of entropy (denoted $H$ usually) terms which are then directly estimated; we term such a paradigm as $\Sigma H$ paradigm. However, the $\Sigma H$ paradigm is applicable only when the variables involved are all discrete or there is a joint density on the space of all variables (in which case, differential entropy $h$ can be utilized). The underlying information measures themselves are well defined in very general probability spaces, for example, involving mixtures of discrete and continuous variables; however, no known estimators exist.

We motivate the requirement of estimators in general probability spaces by some examples in contemporary machine learning and statistical inference.
\begin{enumerate}[leftmargin=*]
\item It is common place in machine learning to have data-sets where {\bf some variables are discrete, and some are continuous}. For example, in recent work on utilizing information bottleneck to understand deep learning \cite{tishby2015deep}, an important step is to quantify the mutual information between the training samples (which are discrete) and the layer output (which is continuous). The employed methodology was to quantize the continuous variables; this is common practice, even though highly sub-optimal.
\item Some variables involved in the calculation {\bf may be mixtures of discrete and continuous variables}. For example, the output of ReLU neuron will not have a density even when the input data has a density. Instead, the neuron will have a discrete mass at $0$ (or wherever the ReLU breakpoint is) but will have a continuous distribution on the positive values. This is also the case in gene expression data, where a gene may have a discrete mass at expression $0$ due to an effect called drop-out \cite{liu2016single} but have a continuous distribution elsewhere.
\item The variables involved may have a joint density {\bf only on a low dimensional manifold}. For example, when calculating mutual information between input and output of a neural network, some of the neurons are deterministic functions of the input variables and hence they will have a joint density supported on a low-dimensional manifold rather than the entire space.
\end{enumerate}

In the aforementioned cases,  no existing estimators are known to work. It is not merely a matter of having provable guarantees either. When we plug in estimators that assume a joint density into data that does not, the estimated information measure can be strongly negative.

We summarize our main contributions below:
\begin{enumerate}[leftmargin=*]
\item {\bf General paradigm} (Section \ref{sec:gdm}): We define a general paradigm of graph divergence measures which captures the aforementioned generalizations of mutual information as special cases. Given a directed acyclic graph (DAG) between $n$ variables, the graph divergence is defined as the Kullback-Leibler (KL) divergence between the true data distribution $\mathbb{P}_X$ and a restricted distribution $\gls{graph_dist}_X$ defined on the Bayesian network and can be thought of as a measure of incompatibility with the given graphical model structure. These graph divergence measures are defined using the Radon Nikodym derivatives which are well-defined for general probability spaces.

\item {\bf Novel estimators} (Section \ref{sec:estimators}): We propose novel estimators for these graph divergence measures, which directly estimate the corresponding Radon-Nikodym derivatives. To the best of our knowledge, these are the first family of estimators that are well defined for general probability spaces (breaking the $\Sigma H$ paradigm).

\item {\bf Consistency proofs} (Section \ref{sec:consistency}): We prove that the proposed estimators converge to the true value of the corresponding graph divergence measures as the number of observed samples increases in a general setting which includes several cases: (1) the data has some discrete and some continuous valued components (2) some (or all) of the components themselves are discrete-continuous mixtures (3) the data is real-valued but does not have a joint density on the entire space but is supported on a low-dimensional manifold.

\item {\bf Numerical results} (Section \ref{sec:experiments}): Extensive numerical results demonstrate that (1) existing algorithms have severe failure modes in general probability spaces (strongly negative values, for example), and (2) our proposed estimator achieves consistency as well as significantly better finite-sample performance.
\end{enumerate}

\begin{figure}
  \begin{subfigure}[t]{0.3\textwidth}
  \centering
    \includegraphics[width=0.5\textwidth]{bayesian-graph.png}
    \caption{}
    \label{fig:bayesian-graph}
  \end{subfigure}
  \begin{subfigure}[t]{0.3\textwidth}
  \centering
    \includegraphics[width=0.5\textwidth]{markov.png}
      \caption{}
    \label{fig:markov}
  \end{subfigure}
  \begin{subfigure}[t]{0.3\textwidth}
  \centering
    \includegraphics[width=0.5\textwidth]{independent.png}
      \caption{}
    \label{fig:independent}
  \end{subfigure}
  \caption{(a) An example of Bayesian Network $\mathcal{G}$ with $\gls{graph_dist}_{X}$ as the induced distribution $\mathbb{P}_{X_1} \mathbb{P}_{X_2}\mathbb{P}_{X_3|X_1}$ $\mathbb{P}_{X_4|{X_1, X_2}}\mathbb{P}_{X_5|X_4}\mathbb{P}_{X_6|X_4}$. (b) A Bayesian Network $\mathcal{G}$ inducing a Markov chain $\mathbb{P}_{X_3} \mathbb{P}_{X_1|X_3}\mathbb{P}_{X_2|X_3}$. (c) A Bayesian Network $\mathcal{G}$ with $\gls{graph_dist}_{X}$ as the induced distribution $\mathbb{P}_{X_1} \mathbb{P}_{X_2}\cdots \mathbb{P}_{X_d}$. }
\end{figure}

\section{Graph Divergence Measure}
\label{sec:gdm}

In this section, we define the family of graph divergence measures. To begin with, we first define some notational preliminaries. We denote any random variable by an uppercase letter such as $X$. The sample space of the variable $X$ is denoted by $\mathcal{X}$ and any value in $\mathcal{X}$ is denoted by the lowercase letter $x$. For any subset $A \subseteq \mathcal{X}$, the probability of $A$ for a given distribution function $\mathbb{P}_X(.)$ over $X$ is denoted by $\mathbb{P}_X(A)$. We note that the random variable $X$ can be a $d$-dimensional \textit{vector} of random variables, i.e. $X \equiv (X_1, \ldots, X_d)$. The $N$ observed samples drawn from the distribution $\mathbb{P}_X$ are denoted by $x^{(1)}, x^{(2)}, \ldots, x^{(N)}$, i.e. $x^{(i)}$ is the $i^{th}$ observed sample.

Sometimes we might be interested in a subset of components of a random variable, $S \subseteq \{ X_1, \ldots, X_d \}$ instead of the entire vector $X$. Accordingly, the sample space of the variable $S$ is denoted by $\mathcal{S}$. For instance, $X=(X_1, X_2, X_3, X_4)$ and $S=(X_1,X_2)$. Throughout the entire paper, unless otherwise stated, there is a one-to-one correspondence between the notations of $X$ and any $S$. For example for any value $x \in \mathcal{X}$, the  corresponding value in $\mathcal{S}$ is simply denoted by $s$. Further,  $s^{(i)} \in \mathcal{S}$ represents the lower-dimensional sample corresponding to the $i$th observed  sample $x^{(i)} \in \mathcal{X}$. Furthermore, any marginal distribution defined over $S$ with respect to $\mathbb{P}_X$ is denoted by $\mathbb{P}_S$.

Consider a directed acyclic graph (DAG) $\mathcal{G}$ defined over $d$ nodes (corresponding to the $d$ components of the random variable $X$). A probability measure $Q$ over $X$ is said to be compatible with the graph $\mathcal{G}$ if it is a Bayesian network on $\mathcal{G}$. Given a graph $\mathcal{G}$ and a distribution $\mathbb{P}_X$, there is a natural measure $\gls{graph_dist}_X(.)$ which is compatible with the graph and is defined as follows:
\begin{equation} \gls{graph_dist}_X = \prod_{l=1}^d \mathbb{P}_{ X_l | \gls{pa_xl_var} } \end{equation}
 where $\gls{pa_xl_var} \subset X$ is the set of the \textit{parent} nodes of the random variable $X_l$, with the sample space denoted by $\gls{pa_xl_space}$, and the sample values $\gls{pa_xl_val}$ corresponding to $x$.  The distribution $\mathbb{P}_{ X_l | \gls{pa_xl_var}} $ is the \textit{conditional} distribution of $X_l$ given $\gls{pa_xl_var}$. Throughout the paper, whenever mentioning the variable $X_l$ with its own parents $\gls{pa_xl_var}$ we indicate it by $\gls{pa_n_xl_var}$, i.e. $\gls{pa_n_xl_var} \equiv \Big( X_l,\gls{pa_xl_var} \Big)$. An example is shown in Fig. \ref{fig:bayesian-graph}.

We note the fact that $ \mathbb{P}_{S | X \setminus S}$ is well defined for any subset of variables $S \subset X$. Therefore if we let $S = X \setminus \gls{pa_xl_var}$, then $\mathbb{P}_{X \setminus \gls{pa_xl_var} | \gls{pa_xl_var}}$ is well defined for any $l\in \{ 1, \ldots, d\}$. By marginalizing over $X \setminus \gls{pa_n_xl_var}$ we see that $\mathbb{P}_{ X_l | \gls{pa_xl_var} }$ and thus the distribution $\gls{graph_dist}_X$ is well defined.

The graph divergence measure is now defined as a function of the probability measure $\mathbb{P}_X$ and the graph $\mathcal{G}$.  In this work we will focus only on the KL Divergence as being the distance metric, hence unless otherwise stated $D(\cdot\parallel\cdot)=D_{KL}(\cdot\parallel\cdot)$. Let's first consider the case where $\mathbb{P}_X$ is absolutely continuous with respect to $\gls{graph_dist}_X$ and hence the \textit{Radon-Nikodym} derivative  $d\mathbb{P}_X/d\gls{graph_dist}_X$ exists. Therefore for a given set of random variables $X$ and a Bayesian Network $\mathcal{G}$, we define \textbf{Graph Divergence Measure ($\gls{graph_distance_measure}$)} as :
\begin{equation}
\gls{graph_distance_measure}(X,\mathcal{G}) = D(\mathbb{P}_X \| \gls{graph_dist}_X) = \int_{\mathcal{X}} \log \frac{d\mathbb{P}_X}{d\gls{graph_dist}_X} d\mathbb{P}_X
\label{eq:measure}\end{equation}
Here we implicitly assume that $\log \big( d\mathbb{P}_X / d\gls{graph_dist}_X \big)$ is measurable and integrable with respect to the measure $\mathbb{P}_X$. The $\gls{graph_distance_measure}$ is set to infinity wherever Radon-Nikodym derivative does not exist. It is clear that $\gls{graph_distance_measure}(X,\mathcal{G})=0$ if and only if the data distribution is compatible with the graphical model, thus the  $\gls{graph_distance_measure}$ can be thought of as a {\em measure of incompatibility with the given graphical model structure}.

We now have relevant variational characterization as below on our graph divergence measure, which can be harnessed to compute upper and lower bounds (More details in \textbf{supplementary material}):
\begin{proposition}
For a random variable $X$, a DAG $\mathcal{G}$, let $\Pi(\mathcal{G})$ be the set of measures $\mathbb{Q}_X$ defined on the Bayesian Network $\mathcal{G}$, then $\gls{graph_distance_measure}(X,\mathcal{G})=\inf_{\mathbb{Q}_X\in\Pi(\mathcal{G})}D(\mathbb{P}_X \| \mathbb{Q}_X) $.

Furthermore, let $\mathcal{C}$ denote the set of functions $h:\mathcal{X}\rightarrow \mathbb{R}$ such that $\mathbb{E}_{\mathbb{Q}_X}\left[\exp(h(X))\right]<\infty$. If $\mathbb{GDM}(X,\mathcal{G})<\infty$, then for every $h\in\mathcal{C}$, $\mathbb{E}_{\mathbb{P}_X}\left[h(X)\right]$ exists and:
\begin{align}
\mathbb{GDM}(X,\mathcal{G})=\sup_{h\in\mathcal{C}} \mathbb{E}_{\mathbb{P}_X}\left[h(X)\right] - \log \mathbb{E}_{\mathbb{Q}_X}\left[\exp(h(X))\right].
\end{align}
\end{proposition}

\vspace{-.15in}
\subsection{Special cases}
\label{sec:math-inf-quantities}

For specific choices of $X$ and Bayesian Network, $\mathcal{G}$, the Equation \ref{eq:measure} is reduced to the well-known information measures. Some examples of these measures are as follows:\\
\textbf{Mutual Information (MI):} $X=(X_1, X_2)$ and $\mathcal{G}$ has no directed edge between $X_1$ and $X_2$. Thus $\gls{graph_dist}_X = \mathbb{P}_{X_1}.\mathbb{P}_{X_2}$, and we get, $\gls{graph_distance_measure}(X,\mathcal{G}) = I(X_1;X_2)  = D(\mathbb{P}_{X_1X_2} \| \mathbb{P}_{X_1}\mathbb{P}_{X_2})$.

\textbf{Conditional Mutual Information (CMI):} We recover the conditional mutual information of $X_1$ and $X_2$ given $X_3$ by constraining $\mathcal{G}$ to be the one in Fig. \ref{fig:markov}, since $\gls{graph_dist}_X = \mathbb{P}_{X_3}.\mathbb{P}_{X_2|X_3}.\mathbb{P}_{X_1|X_3}$, i.e., $\gls{graph_distance_measure}(X,\mathcal{G}) = I(X_1;X_2|X_3)  = D(\mathbb{P}_{X_1X_2X_3} \| \mathbb{P}_{X_1|X_3}\mathbb{P}_{X_2|X_3}\mathbb{P}_{X_3})$.\\

\textbf{Total Correlation (TC):} When $X=(X_1, \cdots, X_d)$, and $\mathcal{G}$ is the graph with no edges (as in Fig. \ref{fig:independent}, we recover the total correlation of the random variables $X_1,\ldots,X_d$ since $\gls{graph_dist}_X = \mathbb{P}_{X_1}\ldots \mathbb{P}_{X_d}$, i.e., $\gls{graph_distance_measure}(X,\mathcal{G}_{dc}) =TC(X_1,\ldots,X_d)  = D(\mathbb{P}_{X_1 \ldots X_d} \| \mathbb{P}_{X_1}\ldots \mathbb{P}_{X_d})$
\\

\textbf{Multivariate Mutual Information (MMI) :} While the MMI defined by \cite{mcgill1954multivariate} is not positive in general,there is a different definition by \cite{chan2015multivariate} which is both non-negative and has an operational interpretation. Since MMI can be defined as the optimal total correlation after clustering, we can utilize our definition to define MMI (cf. \textbf{supplementary material}).

\textbf{Directed Information :} Suppose there are two stationary random processes $X$ and $Y$, the directed information rate from $X$ to $Y$ as first introduced by Massey \cite{massey1990directedinformation} is defined as:
$$ I( X \rightarrow Y ) = \frac{1}{T}\sum_{t=1}^T I \left( X^t ; Y_t \middle| Y^{t-1} \right)$$ \\
It can be seen that the directed information can be written as:
$$ I( X \rightarrow Y ) = \gls{graph_distance_measure}\Big( (X^T,Y^T), \mathcal{G}_I \Big) - \gls{graph_distance_measure}\Big( (X^T,Y^T), \mathcal{G}_C \Big) $$ \\
where the graphical model $\mathcal{G}_I$ correponds to the \textit{independent} distribution between $X^T$ and $Y^T$, and $\mathcal{G}_C$ corresponds to the \textit{causal} distribution from $X$ to $Y$ (more details provided in \textbf{supplementary material}).

\section{Estimators}
\label{sec:estimators}

\subsection{Prior Art}

Estimators for entropy date back to Shannon, who guesstimated the entropy rate of English \cite{shannon1951prediction}. Discrete entropy estimation is a well-studied topic and minimax rate of this problem is well-understood as a function of the alphabet size \cite{paninski2003estimation,jiao2015minimax,wu2016minimax}. The estimation of differential entropy is considerably  harder and also studied extensively in literature \cite{beirlant1997nonparametric, nemenman2002entropy, miller2003new, lee2010sample, lesniewicz2014expected, sricharan2012estimation, singh2014exponential} and can be broadly divided into two groups; based on either Kernel density estimates  \cite{kandasamy2015nonparametric,gao2016breaking} or based on k-nearest-neighbor estimation \cite{sricharan2013ensemble,jiao2017nearest, pal2010estimation}. In a remarkable work, Kozachenko and Leonenko suggested a nearest neighbor method for entropy estimation \cite{kozachenko1987sample} which was then generalized to a $k$th nearest neighbor approach \cite{singh2003nearest}. In this method, the distance to the $k$th nearest neighbor (KNN) is measured for each data-point, and based on this the probability density around each data point is estimated and substituted into the entropy expression. When $k$ is fixed, each density estimate is noisy and the estimator accrues a bias and a remarkable result is that the bias is distribution-independent and can be subtracted out \cite{singh2016finite}.

While the entropy estimation problem is well-studied, mutual information and its generalizations are typically estimated using a sum of signed entropy ($H$) terms, which are estimated first; we term such estimators as $\Sigma H$ estimators. In the discrete alphabet case, this principle has been shown to be worst-case optimal \cite{han2015adaptive}. In the case of distributions with a joint density, an estimator that breaks the $\Sigma H$ principle is the KSG estimator \cite{kraskov2004estimating}, which builds on the KNN estimation paradigm but couples the estimates in order to reduce the bias. This estimator is widely used and has excellent practical performance. The original paper did not have any consistency guarantees and its convergence rates were recently established \cite{pramod_demystify_2018}.  There have been some extensions to the KSG estimator for other information measures such as conditional mutual information \cite{runge2017conditional, frenzel2007partial}, directed information \cite{vejmelka2008inferring} but none of them show theoretical guarantees on consistency of the estimators, furthermore they fail completely in mixture distributions.

When the data distribution is neither discrete nor admits a joint density, the $\Sigma H$ approach is no longer feasible as the individual entropy terms are not well defined. This is the regime of interest in our paper. Recently, Gao et al  \cite{gao2017mixture} proposed a mutual-information estimator based on KNN principle, which can handle such continuous-discrete mixture cases, and the consistency was demonstrated. However it is not clear how it generalizes to even Conditional Mutual Information (CMI) estimation, let alone other generalizations of mutual information. In this paper, we build on that estimator in order to design an estimator for general graph divergence measures and establish its consistency for generic probability spaces.

\subsection{Proposed Estimator}

The proposed estimator is given in Algorithm \ref{algm:gdm} where $\psi(\cdot)$ is the digamma function and $\mathbf{1}_{\{\cdot\}}$ is the indicator function. The process is schematically shown in Fig. \ref{fig:estimator_schematic} (cf. \textbf{supplementary material}). We used the $\ell_\infty$-norm everywhere in our algorithm and proofs.
\begin{algorithm}[htbp]
\caption{Estimating \textbf{Graph Divergence Measure} $\gls{graph_distance_measure}(X, \mathcal{G})$}
\label{algm:gdm}.

\begin{algorithmic}[1]
\Require{\textbf{Parameter:} $k\in\mathbf{Z}^+$, \textbf{Samples:} $x^{(1)}, x^{(2)}, \ldots, x^{(N)}$, \textbf{Bayesian Network:} $\mathcal{G}$ on \textbf{Variables:} ${\cal X}=(X_1,X_2,\cdots,X_d)$}
\Ensure{$\widehat{\gls{graph_distance_measure}}^{(N)}(X,\mathcal{G})$}

 \For{\texttt{$i=1$ to $N$}}
 \State \textbf{Query:}
 \State $\rho_{k,i}$ = $\ell_\infty$-distance to the $k$th nearest neighbor of $x^{(i)}$ in the space $\mathcal{X}$
\State \textbf{Inquire:}
\State  $\tilde{k}_i$ = \# points within the $\rho_{k,i}$-neighborhood of $x^{(i)}$ in the space $\mathcal{X}$
\State $n_{\gls{pa_xl_var}}^{(i)}$ = \# points within the $\rho_{k,i}$-neighborhood of $x^{(i)}$ in the space $\gls{pa_xl_space}$
\State $n_{\gls{pa_n_xl_var}}^{(i)}$ = \# points within the $\rho_{k,i}$-neighborhood of $x^{(i)}$ in the space $\gls{pa_n_xl_space}$
\State \textbf{Compute:}
\State $\zeta_i = \psi(\tilde{k}_i) + \sum_{l=1}^d \left(  \mathbf{1}_{\{  \gls{pa_xl_var} \neq \emptyset \}} \log \left( n_{\gls{pa_xl_var}}^{(i)} +1 \right) - \log \left( n_{\gls{pa_n_xl_var}}^{(i)} +1 \right) \right)$
\EndFor
\State Final Estimator:
\begin{equation}\label{eq:estimator}
\widehat{\gls{graph_distance_measure}}^{(N)}(X,\mathcal{G}) = \frac{1}{N} \sum_{i=1}^N \zeta_i + \left( \sum_{l=1}^d \mathbf{1}_{\{ \gls{pa_xl_var} =\emptyset \}}-1 \right) \log N
\end{equation}
\end{algorithmic}
\end{algorithm}

The estimator intuitively estimates the $\gls{graph_distance_measure}$ by the \textit{resubstitution estimate} $\frac{1}{N} \sum_{i=1}^N \log \hat{f}(x^{(i)})$ in which $\hat{f}(x^{(i)})$ is the estimation of Radon-Nikodym derivative at each sample $x^{(i)}$. If $x^{(i)}$ lies in a region where there is a density, the RN derivative is equal to $g_X(x^{(i)})/\bar{g}_X(x^{(i)})$ in which $g_X(.)$ and $\bar{g}_X(.)$ are density functions corresponding to $\mathbb{P}_X$ and $\gls{graph_dist}_X$ respectively. On the other hand, if  $x^{(i)}$ lies on a point where there is a discrete mass, the RN derivative will be equal to $h_X(x^{(i)})/\bar{h}_X(x^{(i)})$ in which $h_X(.)$ and $\bar{h}_X(.)$ are mass functions corresponding to $\mathbb{P}_X$ and $\gls{graph_dist}_X$ respectively.

The density function $\bar{g}_X(x^{(i)})$ can be written as $\prod_{l=1}^d \left( g_{\gls{pa_n_xl_var}}(\gls{pa_n_xl_val}^{(i)}) / g_{\gls{pa_xl_var}}(\gls{pa_xl_val}^{(i)}) \right)$ for continuous components. Equivalently, the mass function $\bar{h}_X(x^{(i)})$ can be written as $\prod_{l=1}^d \left( h_{\gls{pa_n_xl_var}}(\gls{pa_n_xl_val}^{(i)}) / h_{\gls{pa_xl_var}}(\gls{pa_xl_val}^{(i)}) \right)$. Thus we need to estimate the density functions $g(.)$ and the mass functions $h(.)$ according to the type of $x^{(i)}$. The existing $k$th nearest neighbor algorithms will suffer while estimating the mass functions $h(.)$, since $\rho_{n_S,i}$  (the distance to the $n_S$-th nearest neighbor in subspace $\mathcal{S}$) for such points will be equal to zero for large $N$. Our algorithm, however, is designed in a way that it's capable of approximating both $g(.)$ functions as $\approx \frac{n_S}{N} \frac{1}{(\rho_{n_S,i})^{d_S}}$ and $h(.)$ functions as $\approx \frac{n_S}{N}$ dynamically for any subset $S \subseteq X$. This is achieved by setting $\rho_{n_S,i}$ terms such that all of them cancel out, yielding the estimator as in Eq. (\ref{eq:estimator}).

\section{Proof of Consistency}\label{sec:consistency}

The proof of consistency for our estimator consists of two steps: First we prove that the expected value of the estimator in Eq. (\ref{eq:estimator}) converges to the true value as $N \rightarrow \infty$ , and second we prove that the variance of the estimator converges to zero as $N \rightarrow \infty$.
Let's begin with the definition of $P_X(x,r)$:
\begin{equation} P_X(x,r) = \mathbb{P}_X \big\{ a \in \mathcal{X} : \| a- x \|_{\infty} \leq r \big\} = \mathbb{P}_X \Big\{ B_r(x) \Big\} \label{eq:P_xr}\end{equation}
Thus $P_X(x,r)$ is the probability of a hypercube with the edge length of $2r$ centered at the point $x$. We then state the following assmuptions:

\begin{assumption}\label{assumptions}
We make the following assumptions to prove the consistency of our method:

\begin{enumerate}
\item $k$ is set such that $\lim_{N \rightarrow \infty} k = \infty$ and $\lim_{N \rightarrow \infty} \frac{k \log N}{N}=0$. \label{assumption_kn}
\item The set of discrete points $\{ x : P_X (x, 0) > 0 \}$ is finite. \label{assumption_finitediscrete}
\item $\int_{\mathcal{X}} \big| \log f(x) \big| d\mathbb{P}_X < +\infty$, where $f \equiv d\mathbb{P}_X/d\gls{graph_dist}_X$ is Radon-Nikodym derivative. \label{assumption_logintegrable}
\end{enumerate}

\end{assumption}

The Assumption \ref{assumptions}.\ref{assumption_kn} with \ref{assumptions}.\ref{assumption_finitediscrete} controls the boundary effect between the continuous and the discrete regions; with this assumption we make sure that all the $k$ nearest neighbors of each point belong to the same region almost surely (i.e. all of them are either continuous or discrete).
Assumption \ref{assumptions}.\ref{assumption_logintegrable} is the log-integrability of the Radon-Nikodym derivative. These assumptions are satisfied under mild technical conditions whenever the distribution $\mathbb{P}_X$ over the set $\mathcal{X}$ is (1) finitely discrete; (2) continuous; (3) finitely discrete over some dimensions and continuous over some others; (4) a mixture of the previous cases; (5) has a joint density supported over a lower dimensional manifold. These cases represent almost all the real world data.

As an example of a case not conforming to the aforementioned cases, we can consider singular distributions, among which the \textit{Cantor distribution} is a significant example whose cumulative distribution function is the Cantor function. This distribution has neither a probability density function nor a probability mass function, although its cumulative distribution function is a continuous function. It is thus neither a discrete nor an absolutely continuous probability distribution, nor is it a mixture of these.

The Theorem \ref{thm:mean_convergence} formally states the mean-convergence of the estimator while Theorem \ref{thm:variance_convergence} formally states that convergence of the variance to zero.
\begin{theorem}\label{thm:mean_convergence}
Under the Assumptions \ref{assumptions}, we have $\lim_{N \rightarrow \infty} \mathbb{E} \left[ \widehat{\gls{graph_distance_measure}}^{(N)}(X,\mathcal{G}) \right] = \gls{graph_distance_measure}(X,\mathcal{G})$.
\end{theorem}

\begin{theorem}\label{thm:variance_convergence}
In addition to the Assumptions \ref{assumptions}, assume that we have $( k_N \log N )^2 / N \rightarrow 0 $ as $N$ goes to infinity. Then we have $\lim_{N \rightarrow \infty} \text{Var} \left[ \widehat{\gls{graph_distance_measure}}^{(N)}(X,\mathcal{G}) \right] = 0$.
\end{theorem}

The Theorems \ref{thm:mean_convergence} and \ref{thm:variance_convergence} combined yield the consistency of the estimator \ref{eq:estimator}.

The proof of the Theorem \ref{thm:mean_convergence} starts with writing the Radon-Nikodym derivative explicitly. Then we need to upper-bound the term
$\big| \mathbb{E} \big[ \widehat{\gls{graph_distance_measure}}^{(N)}(X,\mathcal{G}) \big] - \gls{graph_distance_measure}(X,\mathcal{G}) \big|$. To achieve this goal, we segregate the domain of $\mathcal{X}$ into three parts as $\mathcal{X} = \Omega_1 \cup \Omega_2 \cup \Omega_3 $ where  $\Omega_1 = \{x: f(x)=0 \}$,  $\Omega_2 = \{x: f(x)>0, P_X(x,0)>0 \}$ and $\Omega_3 = \{x: f(x)>0, P_X(x,0)=0 \}$.
We will show that $\mathbb{P}_X(\Omega_1)=0$. The sets $\Omega_2$ and $\Omega_3$ correspond to the discrete and continuous regions respectively. Then for each of the two regions, we introduce an upperbound which goes to zero as $N$ grows boundlessly. Thus equivalently we show the mean of the estimate $\zeta_1$ is close to $\log f(x)$ for any $x$.

The proof of the Theorem \ref{thm:variance_convergence} is based on the Efron-Stein inequality, which upperbounds any estimator for any quantity from the observed samples $x^{(1)}, \ldots, x^{(N)}$. For any sample $x^{(i)}$, we then upperbound the term $\left| \zeta_i(X) - \zeta_i(X_{\setminus j}) \right|$ by segregating the samples into various cases, and examining each case separately. $\zeta_i(X)$ is the estimate using all the samples $x^{(1)}, \ldots, x^{(N)}$ and $\zeta_i(X_{\setminus j})$ is the estimate when the $j$th sample is removed. Summing up over all the $i$'s, we obtain an upper-bound which will converge to $0$ as $N$ goes to infinity.

\section{Empirical Results}\label{sec:experiments}

In this section, we evaluate the performance of our proposed estimator in comparison with other estimators via numerical experiments. The estimators evaluated here are our estimator referred to as \textit{GDM}, the plain KSG-based estimators for continuous distributions to which we refer by \textit{KSG}, the \textit{binning} estimators and the noise-induced $\Sigma H$ estimators. A more detailed discussion can be found in Section~\ref{sec:exp_details}.

\textbf{Experiment 1: Markov chain model with continuous-discrete mixture.}
For the first experiment, we simulated an $X$-$Z$-$Y$ Markov chain model in which the random variable $X$ is a uniform random variable $\mathcal{U}(0,1)$ clipped at a threshold $0 <\alpha_1 < 1$ from above. Then $Z = \min \left( X, \alpha_2 \right)$ and $Y = \min \left( Z, \alpha_3 \right)$ in which $ 0 < \alpha_3 <\alpha_2 <\alpha_1$.
We simulated this system for various numbers of samples, setting $\alpha_1=0.9$, $\alpha_2=0.8$ and $\alpha_3=0.7$. For each set of samples we estimated $I(X;Y|Z)$ via different methods. The theory value for $I(X;Y|Z)$ is $0$. The results are shown in Figure \ref{fig:exp1_cmi_vs_n}. We can see that in this regime, only the GDM estimator can correctly converge. The KSG estimator and the $\Sigma H$ estimator show high negative biases and the binning estimator shows a positive bias.

\textbf{Experiment 2: Mixture of AWGN and BSC channels with variable error probability.}
For the second scheme of our experiments, we considered an Additive White Gaussian Noise (AWGN) Channel in parallel with a Binary Symmetric Channel (BSC) where only one of the two can be activated at a time. The random variable $Z = \min ( \alpha, \tilde{Z} )$ where $\tilde{Z} \sim U(0,1)$ controls which channel is activated; i.e. if $Z$ is lower than the threshold $\beta$, activate the AWGN channel, otherwise initiate the BSC channel where $Z$ also determines the error probability at each time point. We set $\alpha=0.3$, $\beta=0.2$, BSC channel input as $X \sim \text{Bern}(0.5)$, and AWGN input and noise deviation as $\sigma_X=1$ and $\sigma_N=0.1$ respectively, and obtained the estimates of $I(X;Y|Z,Z^2,Z^3)$ for various estimators. While the theory value is equal to $I(X;Y|Z)=0.53241$, yet it's conditioned over a low-dimensional manifold in a high-dimensional space. The results are shown in Figure \ref{fig:exp2_cmi_ld_vs_n}. Similar to the previous experiment, the GDM estimator can correctly converge to the true value. The $\Sigma H$ and binning estimators show a negative bias, and the KSG estimator gets totally lost.

\textbf{Experiment 3: Total Correlation for independent mixtures.}
In this experiment, we estimate the total correlation of three independent variables $X$, $Y$ and $Z$. The samples for the variable $X$ are generated in the following fashion:
First toss a fair coin, if heads appears we fix $X$ at $\alpha_X$, otherwise we draw $X$ from a uniform distribution between $0$ and $1$. samples from $Y$ and $Z$ are also generated in the same way independently with parameters $\alpha_Y$ and $\alpha_Z$ respectively. For this setup, $TC(X,Y,Z)=0$. We set $\alpha_X=1$, $\alpha_Y=1/2$ and $\alpha_Z=1/4$, and generated various datasets with different lengths. Then estimated total correlation values are shown in the Figure \ref{fig:exp3_tc_vs_n}.

\textbf{Experiment 4: Total Correlation for independent uniforms with correlated zero-inflation.}
Here we first consider four auxiliary uniform variables $\tilde{X}_1$, $\tilde{X}_2$, $\tilde{X}_3$ and $\tilde{X}_4$ which are taken from $\mathcal{U}(0.5,1.5)$. Then each sample is erased with a Bernoulli probability; i.e. $X_1=\alpha_1 \tilde{X}_1$, $X_2=\alpha_1 \tilde{X}_2$ and $X_3=\alpha_2 \tilde{X}_3$, $X_4=\alpha_2 \tilde{X}_4$   in which $\alpha_1 \sim \text{Bern}(p_1)$ and $\alpha_2 \sim \text{Bern}(p_2)$. As we see, after zero-inflation $X_1$ and $X_2$ become correlated, so do $X_3$ and $X_4$ while still $(X_1,X_2) \rotatebox[origin=c]{90}{$\models$} (X_3,X_4)$.
 In the experiment, we set $p_1=p_2=0.6$. The results of running different algorithms over the data can be seen in Figure \ref{fig:exp4_tc_vs_n}. For the total correlation experiments 3 and 4, similar to that of conditional mutual information in experiments 1 and 2, only the GDM estimator can best estimate the true value. The estimator $\Sigma H$ was removed from the figures due to its high inaccuracy.

\begin{figure}
\centering
	\begin{subfigure}[t]{.45\textwidth}
 	\centering
	\includegraphics[width=\textwidth,trim={.3cm 0 1cm 0},clip]{3node_truncate.pdf}
	\caption{}
	\label{fig:exp1_cmi_vs_n}
	\end{subfigure}
	\begin{subfigure}[t]{.45\textwidth}
	\centering
	\includegraphics[width=\textwidth,trim={.3cm 0 1cm 0},clip]{AWGN_BSC_low_dimension.pdf}
	\caption{}
	\label{fig:exp2_cmi_ld_vs_n}
	\end{subfigure}

	\begin{subfigure}[t]{.45\textwidth}
	\centering
	\includegraphics[width=\textwidth,trim={.3cm 0 1cm 0},clip]{3node_independent_mixture.pdf}
	\caption{}
	\label{fig:exp3_tc_vs_n}
	\end{subfigure}
	\begin{subfigure}[t]{.45\textwidth}
	\centering
	\includegraphics[width=\textwidth,trim={.3cm 0 1cm 0},clip]{4node_zeroinflation.pdf}
	\caption{}
	\label{fig:exp4_tc_vs_n}
	\end{subfigure}

	\begin{subfigure}[t]{.45\textwidth}
	\centering
	\includegraphics[width=\textwidth,trim={.3cm 0 1cm 0},clip]{GRN_auc_vs_N.pdf}
	\caption{}
	\label{fig:grn_auc_vs_n}
	\end{subfigure}
	\begin{subfigure}[t]{.45\textwidth}
	\centering
	\includegraphics[width=\textwidth,trim={.3cm 0 1cm 0},clip]{FeatureSelection.pdf}
	\caption{}
	\label{fig:featureselection_auc_vs_n}
	\end{subfigure}

\caption{The results for the experiments versus the number of samples: \ref{fig:exp1_cmi_vs_n}: The estimated CMI for the X-Z-Y Markov chain. \ref{fig:exp2_cmi_ld_vs_n}: CMI for the AWGN+BSC channels with low-dimensional $Z$ manifold. \ref{fig:exp3_tc_vs_n}: The estimated TC values for three independent variables. \ref{fig:exp4_tc_vs_n}: The estimated TC for zero-inflated variables. \ref{fig:grn_auc_vs_n}: The AUROC values for gene regulatory network inference. The error bars show the standard deviation scaled down by $0.2$. \ref{fig:featureselection_auc_vs_n}: The AUROC values for feature selection accuracy. The error bars show the standard deviations scaled down by $0.2$.}
\end{figure}

\textbf{Experiment 5: Gene Regulatory Networks.}
In this experiment we use different estimators to do Gene Regulatory Network inference based on the conditional Restricted Directed Information (cRDI) \cite{rahimzamani2016network}. We do our test on the simulated neuron cells' development process, based on a model from \cite{qiu2012understanding}. In this model, the time series vector $X$ consists of $13$ random variables each of which corresponding to a single gene in the development process. We simulated the development process for various lengths of time-series in which the noise $N \sim \mathcal{N}(0,.03)$ is added for all the genes, and every single sample is then subject to erasure (i.e. be replaced by 0s) with a probability of $0.5$. Then we applied the cRDI method utilizing various CMI estimators and then calculated the Area-Under-ROC curve (AUROC). The results are shown in Figure \ref{fig:grn_auc_vs_n}. It's seen that the cRDI method implemented with the GDM estimator outperform the other estimators by at least $\%10$ in terms of AUROC. In the tests, cRDI for each $(X_i,X_j)$ is conditioned over the node $k \neq i$ with the highest RDI value to $j$. We notice that the causal signals are highly destroyed due to the zero-inflation, so we won't expect high performance of the causal inference over the data. We did not include the $\Sigma H$ estimator results due to its very low performance.

\textbf{Experiment 6: Feature Selection by Conditional Mutual Information Maximization.}
Feature selection is an important pre-processing step in many learning tasks. The application of information theoretic measures in feature selection is well studied in the literature \cite{vergara2014review}. Among the well-known methods is the conditional mutual information maximization (CMIM) first introduced by Flueret \cite{fleuret2004fast}, a variation of which was later introduced called CMIM-2 \cite{vergara2010cmim2}. Both methods use conditional mutual information as their core measure to select the features. Hence the performance of the estimators can significantly influence the performance of the methods. In our experiment, we generated a vector $X=(X_1, \ldots, X_{15})$ of 15 random variables in which all the random variables are taken from $\mathcal{N}(0,1)$ and then each random variable $X_i$ is clipped from above at $\alpha_i$ which is initially taken randomly from $\mathcal{U}(0.25,0.3)$ and then kept constant during the sample generation. Then $Y$ is generated as $Y=\cos \big( \sum_{i=1}^5 X_i \big)$. Then we did the CMIM-2 algorithm with various CMI estimators to evaluate the performance of the estimators in extracting the relevant features $X_1, \ldots, X_5$. The AUROC values for each algorithm versus the number of samples generated are shown in Figure \ref{fig:featureselection_auc_vs_n}. The feature selection methods implemented with the GDM estimator outperform the other estimators.

\section{Discussion and Future Work}
A general paradigm of graph divergence measures and novel estimators thereof, for general probability spaces are proposed, which estimate several generalizations of mutual information. In future, we would like to derive further efficient estimators for high dimensional data. In the current work, estimators are shown to be consistent with infinite scaling of parameter $k$. In future, we would like to understand the finite $k$ performance of the estimators as well as guarantees on sample complexity and rates of convergence. Another potential direction to follow is to study the variational characterization of the graph divergence measure to design estimators. Improving the computational efficiency of the estimator is another direction of future work. Recent literature including \cite{noshad2018scalable} provide a new methodology to estimate mutual information in a computationally efficient manner and leveraging these ideas for the generalized measures and general proabability distributions can be a promising direction ahead.

\section{Acknowledgement}
This work was partially supported by NSF grants 1651236, 1703403 and NIH grant 5R01HG008164.

\newpage