Nonparametric Estimation of Renyi Divergence and Friends

Introduction

Given samples from two distributions, one fundamental and classical question to ask is: how close are the two distributions? First, one must specify what it means for two distributions to be close, for which a number of divergences have been proposed. Then there is the statistical question: how does one estimate divergence given samples from two distributions. In this paper, we propose and analyze estimators for three common divergences.

Divergence estimation has a number of applications across machine learning and statistics. In statistics, one can use these estimators to construct two-sample and independence tests . In machine learning, it is often convenient to view training data as a set of distributions and use divergences to estimate dissimilarity between examples. This idea has been used in neuroscience, where the neural response pattern of an individual is modeled as a distribution, and divergence is used to compare responses across subjects . It has also enjoyed success in computer vision, where features are computed for each patch of an image and these feature vectors are modeled as independent draws from an underlying distribution .

For these applications and others, it is crucial to accurately estimate divergences given samples drawn independently from each distribution. In the nonparametric setting, a number of authors have proposed various estimators which are provably consistent. However, apart from a few examples, the actual rates of convergence of these estimators and the minimax optimal rates are still unknown.

In this work, we propose three estimators for the $L_{2}^{2}$, Rényi-$\alpha$, and Tsallis-$\alpha$ divergence between two continuous distributions. Our strategy is to correct an initial plug-in estimator by estimates of the higher order terms in the von Mises expansion of the divergence functional. We establish the rates of convergence for these estimators under the assumption that both densities belong to a Hölder class of smoothness $s$. Concretely, we show that the plug-in estimator achieves rate $n^{\frac{-s}{2s+d}}$ while correcting by the first order terms in the expansion results in an $n^{-\min\{\frac{2s}{2s+d},1/2\}}$-estimator and correcting further by the second order terms gives an $n^{-\min\{\frac{3s}{2s+d},1/2\}}$-estimator. These last two estimators achieve the parametric $n^{-1/2}$ rate as long as the smoothness $s$ is larger than $d/2,d/4$, respectively, where $d$ is the dimension. Moreover the first-order estimator, while worse statistically than the second-order estimator, is computationally very elegant. These results contribute to our fairly limited knowledge on this important problems .

We also address the issue of statistical optimality by deriving a minimax lower bound on the convergence rate. Specifically, we show that one cannot estimate these quantities at better than $n^{\frac{-4s}{4s+d}}$-rate when $s\leq d/4$ and $n^{-1/2}$-rate otherwise. This establishes the optimality of our best estimator in the smooth regime and also that $d/4$ is the critical smoothness for this problem.

The remainder of this manuscript is organized as follows. After discussing some related work on divergence estimation and the closely-related entropy estimation in Section 2, we present our estimators and main results in Sections 3 and 4. We provide proof sketches in Section 5. We present some numerical simulations in Section 6 and conclude with some open questions in Section 7. We defer many proof details and several calculations to the appendices.

Technically, these divergences are functionals on distributions, rather than densities, but we will abuse notation and write them as above. As a unification, we consider estimating functionals of the form, $T(p,q)=\int p^{\alpha}(x)q^{\beta}(x)d\mu(x)$ for given $\alpha,\beta$. Various settings of $\alpha,\beta$ yield the main terms in the divergences, and we will verify that estimators for $T(p,q)$ result in good divergence estimators.

The sine qua non of our work is the von Mises expansionSee Chapter 20 of van der Vaart’s book for an introduction to von Mises calculus .. Given a functional $T$ mapping distributions to the reals, the first-order von Mises expansion is:

where $F$ and $G$ are distributions, $R_{2}$ is a remainder term, and $dT(G;F-G)$ is the Gateaux derivative of $T$ at $G$ in the direction of $F-G$:

In our work, $T$ is always of the form $T(F)=\int\phi(f)d\mu$ where $f=dF/d\mu$ is the Radon-Nikodym derivative and $\phi$ is differentiable. In this case, the von Mises expansion reduces to a functional Taylor expansion on the densitiesSee Lemma 8 in the Appendix.:

We generalize these ideas to functionals of two distributions and with higher order expansions analogous to the Taylor expansion. We often write $T(f)$ instead of $T(F)$.

Related Work

Divergence estimation and its applications have received considerable attention over the past several decades. Pardo provides a fairly comprehensive discussion of methods and applications in the context of discrete distributions .

Only recently has attention shifted to the continuous, nonparametric setting, where a number of efforts have established consistent estimators. Many of the approaches are based on nearest-neighbor graphs . For example, Póczos and Schneider use a $k$-nearest-neighbor estimator and show that one does not need a consistent density estimator to consistently estimate Rényi-$\alpha$ and Tsallis-$\alpha$ divergences. A number of other authors have also proposed consistent estimators via the empirical CDF or histograms . Unfortunately, the rates of convergence for all of these methods are still unknown.

Singh and Poczos recently established a rate of convergence for an estimator based on simply plugging kernel density estimates into the divergence functional. Their estimator converges at $n^{\frac{-s}{s+d}}$-rate when $s<d$ and $n^{-1/2}$ otherwise which matches some existing results on estimating entropy functionals . In comparison, we show that corrections of the plug-in estimator lead to faster convergence rates and that the $n^{-1/2}$ rate can be achieved at the much lower smoothness of $s>d/4$. Moreover we establish a minimax lower bound for this problem, which shows that $d/4$ is the critical smoothness index.

Källberg and Seleznjev study an $\epsilon$-nearest neighbor estimator for the $L_{2}^{2}$-divergence that enjoys the same rate of convergence as our projection-based estimator. They prove that the estimator is asymptotically normal in the $s>d/4$ regime, which one can also show for our estimator. In the more general setting of estimating polynomial functionals of the densities, they only show consistency of their estimator, while we also characterize the convegence rate.

A related and flourishing line of work is on estimating entropy functionals. The majority of the methods are graph-based, involving either nearest neighbor graphs or spanning trees over the data . One exception is the KDE-based estimator for mutual information and joint entropy of Liu, Lafferty, and Wasserman . A number of these estimators come with provable convergence rates.

While it is not clear how to port these ideas to divergence estimation, it is still worth comparing rates. The estimator of Liu et al. converges at rate $n^{\frac{-s}{s+d}}$, achieving the parametric rate when $s>d$. Similarly, Sricharan et al. show that when $s>d$ a $k$-NN style estimator achieves rate $n^{-2/d}$ (in absolute error) ignoring logarithmic factors. In a follow up work, the authors improve this result to $O(n^{-1/2})$ using an ensemble of weak estimators, but they require $s>d$ orders of smoothness . In contrast, our estimators achieve the parametric $n^{-1/2}$ rate at lower smoothness ($s>d/2,d/4$ for the first-order and second-order estimators, respectively) and enjoy a faster rate of convergence uniformly over smoothness.

Interestingly, while many of these methods are plug-in-based, the choice of tuning parameter typically is sub-optimal for density estimation. This contrasts with our technique of correcting optimal density estimators.

We are not aware of any lower bounds for divergence estimation, although analogous results have been established for the entropy estimation problem. Specifically, Birgé and Massart prove a $n^{\frac{-4s}{4s+d}}$-lower bound for estimating integral functionals of a density. Hero et al. give a matching lower bound for estimating Rényi-$\alpha$ entropies.

Finally, our estimators and proof techniques are based on several classical works on estimating integral functionals of a density. The goal here is to estimate $\int\phi(f(x))d\mu(x)$, for some known function $\phi$, given samples from $f$. A series of papers show that $n^{-1/2}$ rate of convergence is attainable if and only if $s>d/4$, which is analogous to our results . Of course, our results pertain to the two-density setting, which encompasses the divergences of interest. We also generalize some of these results to the multi-dimensional setting.

The Estimators

Recall that we are interested in estimating integral functionals of the form $T(p,q)=\int p^{\alpha}(x)q^{\beta}(x)$. As an initial attempt, with estimators $\hat{p}$ and $\hat{q}$ for $p$ and $q$, we can use the plug-in estimator $\widehat{T}_{pl}=T(\hat{p},\hat{q})$. Via the von Mises expansion of $T(p,q)$, the error is of the form:

Classical results on density estimation then suggest that $\widehat{T}_{pl}$ will enjoy a $n^{\frac{-s}{2s+d}}$-rate .

A better convergence rate can be achieved by correcting the plug-in estimator with estimates of the linear term in the von Mises expansion. Informally speaking, the remainder of the first order expansion is $O(\|\hat{p}-p\|_{2}^{2}+\|\hat{q}-q\|_{2}^{2})$ which decays with $n^{\frac{-2s}{2s+d}}$, while the linear terms can be estimated at $n^{-1/2}$-rate. This estimator, which we call $\widehat{T}_{lin}$ enjoys a faster convergence rate than $\widehat{T}_{pl}$.

It is even better to augment the plug-in estimator with both the first and second-order terms of the expansion. Here the remainder decays at rate $n^{-\frac{3s}{2s+d}}$ while the linear and quadratic terms can be estimated at $n^{-1/2}$ and $n^{\frac{-4s}{4s+d}}$ rate respectively. This corrected estimator $\widehat{T}_{quad}$ achieves the parametric rate whenever the smoothness $s>d/4$ which we will show to be minimax optimal.

We now formalize these heuristic developmentsSee Appendices A and B for details omitted in this section.. Below we enumerate the terms in the first and second order von Mises expansions that we will estimate or compute:

These definitions allow us to succinctly write the expansions of $T(p,q)$ about $T(\hat{p},\hat{q})$:

with remainders, $R_{a}=O(\|p-\hat{p}\|_{a}^{a}+\|q-\hat{q}\|_{a}^{a})$.

The terms $\theta_{(\cdot),2}^{(\cdot)}$ are of the form:

again with known $\psi$. To estimate these terms, we have samples $X_{1}^{n}\sim f,Y_{1}^{n}\sim g$. If $\{\phi_{k}\}_{k\in D}$ is an orthonormal basis for $L_{2}(^{d})$ then the estimator for the bilinear term is:

where $M\subset D$ is chosen to tradeoff the bias and the variance. To develop some intuition, if we knew $f$, we would simply use the sample mean $\frac{1}{n}\sum_{j=1}^{n}f(Y_{j})\psi(Y_{j})$. Since $f$ is actually unknown, we replace it with an estimator formed by truncating its Fourier expansion. Specifically, we replace $f$ with $\hat{f}(\cdot)=\sum_{k\in M}\hat{a}_{k}\phi_{k}(\cdot)$ with $\hat{a}_{k}=\frac{1}{n}\sum_{i=1}^{n}\phi_{k}(X_{i})$.

For the quadratic functional, a projection estimator was proposed and analyzed by Laurent :

where $b_{k,k^{\prime}}(\psi)=\int\phi_{k}(x)\phi_{k^{\prime}}(x)\psi(x)dx$. The first term in the estimator is motivated by the same line of reasoning as in the bilinear estimator while the second term significantly reduces the bias without impacting the variance.

Before proceeding to our theoretical analysis, we mention some algorithmic considerations. We estimate $\hat{p},\hat{q}$ with kernel density estimators, which, except for in $\widehat{T}_{pl}$, we only train on half of the sample. This gives us independent samples to estimate the $\hat{\theta}_{(\cdot),(\cdot)}^{(\cdot)}$ terms. Second, in our analysis, we will require that the KDEs are bounded above and below. Under the assumption that $p$ and $q$ are bounded above and below, we will show that clipping the original KDE will not affect the convergence rate.

Another important issue with density estimation over bounded domains, that applies to our setting, is that the standard KDE suffers high bias near the boundary. To correct this bias, we adopt the strategy used by Liu et al. of “mirroring” the data set over the boundaries. We do not dwell too much on this issue, noting that this technique can be shown to suitably correct for boundary bias without substantially increasing the variance. This augmented estimator can be shown to match the rates of convergence in the literature .

Lastly, the estimators all require integration of the term $T(\hat{p},\hat{q})$, which can be computationally burdensome, particularly in high dimension. However, whenever $\alpha+\beta=1$, as in the Rényi-$\alpha$ and Tsallis-$\alpha$ divergences, the constants $C_{1},C_{2}$ are zero, so the first term may be omitted. In this case $\widehat{T}_{lin}$ is remarkably simple; it involves training KDEs and estimating a specific linear functional of them via the sample mean. Although this estimator is not minimax optimal, it enjoys a fairly fast rate of convergence while being computationally practical. Unfortunately, even when $C_{2}=0$, the quadratic estimator still involves integration of the $b_{i,i^{\prime}}$ terms. We therefore advocate for $\widehat{T}_{lin}$ over $\widehat{T}_{quad}$ in practice, as $\widehat{T}_{lin}$ exhibits a better tradeoff between computational and statistical efficiency.

Theoretical Results

For our theoretical analysis, we will assume that the densities $p,q$ belong to $\Sigma(s,L)$, the periodic Hölder class of smoothness $s$, defined as follows:

For any tuple $r=(r_{1},\ldots,r_{d})$ define $D^{r}=\frac{\partial^{r_{1}+\ldots+r_{d}}}{\partial x_{1}^{r_{1}}\ldots\partial x_{d}^{r_{d}}}$. The periodic Hölder class $\Sigma(s,L)$ is the subset of $L_{2}(^{d})$ where for each $f\in\Sigma(s,L)$, the $r$th derivative is periodic for any tuple $r$ with $\sum_{j}r_{j}<s$ and:

for all $x,y$ and for all tuples $r$ with $\sum_{j}r_{j}=\lfloor s\rfloor$ the largest integer strictly smaller than $s$.

We are now ready to state our main assumptions:

$p,q\in\Sigma(s,L)$ for some known smoothness $s$.

The densities are bounded above and below by known parameters $\kappa_{l},\kappa_{u}$. Formally $0<\kappa_{l}\leq p(x),q(x)\leq\kappa_{u}<\infty$ for all $x\in^{d}$.

Set the KDE bandwidth $h\asymp n^{\frac{-1}{2s+d}}$. For any projection-style estimator, set the number of basis elements $m\asymp n^{\frac{2d}{4s+d}}$.

The Hölderian assumption is standard in the nonparametric literature while the periodic assumption subsumes more standard boundary smoothness conditions . It is fairly straightforward to construct kernels meeting Assumption 3 , while the boundedness assumption is common in the literature on estimating integral functionals of a density .

The following theorem characterizes the rate of convergence of our estimators $\widehat{T}_{pl},\widehat{T}_{lin},\widehat{T}_{quad}$:

All expectations are taken with respect to $X_{1}^{n},Y_{1}^{n}$. When $s=d/4$, $\widehat{T}_{quad}$ enjoys $O(n^{-1/2+\epsilon})$ rate of convergence for any $\epsilon>0$The constant is exponential in $\epsilon$ and is infinite for $\epsilon=0$.. $\widehat{T}_{lin}$ and $\widehat{T}_{quad}$ achieve the parametric rate when $s>d/2,d/4$ respectively.

Before commenting on the upper bound and presenting some consequences, we address the question of statistical efficiency. Clearly $\widehat{T}_{pl}$ and $\widehat{T}_{lin}$ are not rate-optimal, since $\widehat{T}_{quad}$ achieves a faster rate of convergence, but is $\widehat{T}_{quad}$ minimax optimal? We make some progress in this direction with a minimax lower bound on the rate of convergence.

Under Assumptions 1 and 2, as long as both $\alpha,\beta\neq 0,1$, then with $\gamma_{\star}=\min\{4s/(4s+d),1/2\}$ and for any $\epsilon>0$:

For a pictorial understanding of the rates of convergence and the lower bound, we plot the exponent $\gamma$ for each of the terms in the von Mises expansion as a function of the smoothness $s$ in Figure 1. The estimator $\widehat{T}_{quad}$ has three terms, with rates $n^{-1/2},n^{\frac{-4s}{4s+d}}$, and $n^{\frac{-3s}{2s+d}}$ respectively which achieves the parametric rate $n^{-1/2}$ when $s>d/4$ and is $n^{\frac{-3s}{2s+d}}$ in the low-smoothness regime. The linear estimator only achieves the parametric rate while $s>d/2$ while $\widehat{T}_{pl}$ only approaches the parametric rate as $s\rightarrow\infty$. Consequently these estimators are statistically inferior to $\widehat{T}_{quad}$. In the last plot we show a lower bound on the rate of convergence from Theorem 2, which is $n^{\frac{-4s}{4s+d}}$ when $s\leq d/4$ and $n^{-1/2}$ when $s>d/4$.

The lower bound rate deviates slightly from the upper bound for $\widehat{T}_{quad}$ in the low-smoothness regime, showing that $\widehat{T}_{quad}$ is also not minimax-optimal uniformly over $s$. This sub-optimality appears even when estimating integral functionals of a single density . In that context, achieving the optimal rate of convergence in the non-smooth regime involves further correction by the third order term in the expansion . It seems as if the same ideas can be adapted to the two-density setting, although we believe computational considerations would render these estimators impractical.

In the smooth regime ($s>d/4$) we see that the parametric $n^{-1/2}$ rate is both necessary and sufficient. This critical smoothness index of $s=d/4$ was also observed in the context of estimating integral functionals of densities .

When $s=d/4$, the quadratic estimator achieves $n^{-1/2+\epsilon}$ rate for any $\epsilon>0$, where the constant is exponential in $\epsilon$, and thus deviates slightly from the lower bound. This phenomenon arises from using the projection-based estimators for the quadratic term. Establishing the rate of convergence for these estimators requires working in a Sobolev space rather than the Hölder class. In translating back to the Hölderian assumption, we lose a small factor in the smoothness, since the Sobolev space only contains the Hölder space if the former is less smooth than the latter.

The lower bound on estimating integral functionals in Theorem 2 almost immediately implies a lower bound for Tsallis-$\alpha$ divergences. For Rényi-$\alpha$, some care must be taken in the translation, but we are able to prove the same lower bound as long as $D_{\alpha}(p,q)$ is bounded. The idea behind these extensions is to translate an estimator $\hat{D}$ for the divergence into an estimator $\hat{T}$ for $T(p,q)$. We then argue that if $\hat{D}$ enjoyed a fast rate of convergence, so would $\hat{T}$, which leads to a contradiction of the theorem. Unfortunately, Theorem 2 does not imply a lower bound for $L_{2}^{2}$ divergence, since we are unable to handle the $\alpha=\beta=1$ case, which is exactly the cross term in the $L_{2}^{2}$-divergence.

Our proof requires that both $\alpha,\beta$ are both not or $1$, which is not entirely surprising. If $\alpha=\beta=0$, $T(p,q)$ is identically zero, so one should not be able to prove a lower bound. Similarly $\alpha=0,\beta=1$ or vice versa, $T(p,q)=1$ for any $p,q$, so we have efficient, trivial estimators.

The only non-trivial case is $\alpha=\beta=1$ and we conjecture that the $n^{-\gamma_{\star}}$ rate is minimax optimal there, although our proof does not apply. Our proof strategy involves fixing $q$ and perturbing $p$, or vice versa. In this approach, one can view the optimal estimator as having knowledge of $q$, so if $\alpha=1$, the sample average is a $n^{-1/2}$-consistent estimator, which prevents us from achieving the $n^{-\gamma_{\star}}$ rate. We believe this is an artifact of our proof, and by perturbing both $p$ and $q$ simultaneously, we conjecture that one can prove a minimax lower bound of $n^{-\gamma_{\star}}$ when $\alpha=\beta=1$.

We now show how an estimate of $T(p,q)$ can be used to estimate the divergences mentioned above. Plugging $\widehat{T}_{quad}$ into the definition of Rényi-$\alpha$ and Tsallis-$\alpha$ divergences, we immediately have the following corollary:

Under Assumptions 1- 4, as long as $D_{\alpha}(p,q)\geq c>0$ for some constant $c$, the estimators:

As we mentioned before, when $\beta=1-\alpha$, for both the linear and quadratic estimators, one can omit the term $T(\hat{p},\hat{q})$ as the constants $C_{1},C_{2}=0$. However, $\widehat{T}_{quad}$ is still somewhat impractical due to the numeric integration in the quadratic terms. On the other hand, the linear estimator $\widehat{T}_{lin}$ is computationally very simple, although its convergence rate is $O(n^{-1/2}+n^{\frac{-2s}{2s+d}})$.

For the $L_{2}^{2}$ divergence, instead of applying Theorem 4 directly, it is better to directly use the quadratic and bilinear estimators for the terms in the factorization. Specifically, let $\theta_{p}=\int p^{2}$ and define $\hat{\theta}_{p}$ by Equation 2 with $\psi(x)=1$. Define $\theta_{q},\hat{\theta}_{q}$ analogously and finally define $\theta_{p,q}=2\int pq$ with $\hat{\theta}_{p,q}$ given by Equation 1 where $\psi(x)=2$. As a corollary of Theorem 6 below, we have:

Under Assumptions 1- 4, the estimator $\hat{L}=\hat{\theta}_{p}+\hat{\theta}_{q}-\hat{\theta}_{p,q}$ for $L_{2}^{2}(p,q)$ satisfies:

Notice that for both quadratic terms, the $b_{i,i^{\prime}}$ terms in Equation 2 are $\mathbf{1}[i=i^{\prime}]$ since $\psi(x)=1$ and since $\{\phi_{k}\}$ is an orthonormal collection. Thus the estimator $\hat{L}$ is computationally attractive, as numeric integration is unnecessary. In addition, we do not need KDEs, removing the need for bandwidth selection, although we still must select the basis functions used in the projection.

Proof Sketches

The rates of convergence for $\widehat{T}_{pl},\widehat{T}_{lin}$, and $\widehat{T}_{quad}$ come from analyzing the kernel density estimators and the estimators for $\hat{\theta}_{(\cdot),(\cdot)}^{(\cdot)}$. Recall that we must use truncated KDEs $\hat{p},\hat{q}$ with boundary correction, so standard analysis does not immediately apply. However, we do have the following theorem establishing that truncation does not affect the rate, which generalizes previous results to high dimension .

Let $f$ be a density satisfying Assumptions 1- 4 and suppose we have $X_{1}^{n}\sim f$. The truncated KDE $\hat{f}_{n}$ satisfies:

Let $f,g$ be densities in $\Sigma(s,L)$ and let $\psi$ be a known bounded function. Let $\phi_{k}$ be the Fourier basis and $M$ the set of basis elements with frequency not exceeding $m_{0}^{1/d}$, where $m_{0}\asymp n^{\frac{2d}{4s^{\prime}+d}}$ for some $s^{\prime}<s$. If $\theta=\int\psi(x)f(x)g(x)$ and $\hat{\theta}$ is given by Equation 1 or if $\theta=\int\psi(x)f^{2}(x)$ and $\hat{\theta}$ is given by Equation 2, then:

Theorem 4 follows from these results, the von Mises expansion, and the triangle inequality.

2 Lower Bound

The first part of the lower bound is an application of Le Cam’s method and generalizes a proof of Birge and Massart . We begin by reducing the estimation problem to a simple-vs.-simple hypothesis testing problem. We will use the squared Hellinger distance, defined as:

Let $T$ be a functional defined on some subset of a parameter space $\Theta\times\Theta$ which contains $(p,q)$ and $(g_{\lambda},q)\forall\lambda$ in some index set $\Lambda$. Define $\bar{G}^{n}=\frac{1}{|\Lambda|}\sum_{\lambda\in\Lambda}G^{n}_{\lambda}$ where $G_{\lambda}$ has density $g_{\lambda}$. If:

where $c_{\gamma}=\frac{1}{2}[1-\sqrt{\gamma(1-\gamma/4)}]$.

To construct the $g_{\lambda}$ functions, we partition the space $^{d}$ into $m$ cubes $R_{j}$ and construct functions $u_{j}$ that are compactly supported on $R_{j}$. We then set $g_{\lambda}=p+K\sum_{j=1}^{m}\lambda_{j}u_{j}$ for $\lambda\in\Lambda=\{-1,1\}^{m}$. By appropriately selecting the functions $u_{j}$, we can ensure that:

Ensuring smoothness requires $K=O(m^{-s/d})$ at which point, making the Hellinger distance $O(1)$ requires $m=\Omega(n^{\frac{2d}{4s+d}})$. With these choices we can apply Lemma 7 and arrive at the lower bound since $K^{2}=m^{-2s/d}=n^{\frac{-4s}{4s+d}}$.

As for the second part of the theorem, the $n^{-1/2}$ lower bound, we use a (to our knowledge) novel proof technique which we believe may be applicable in other settings. The first ingredient of our proof is a lower bound showing that one cannot estimate a wide class of quadratic functionals at better than $n^{-1/2}$ rate. We provide a proof of this result based on Le Cam’s method in the appendix although related results appear in the literature . Then starting with the premise that there exists an estimator $\hat{T}$ for $T(p,q)$ with rate $n^{-1/2-\epsilon}$, we construct an estimator for a particular quadratic functional with $n^{-1/2-\epsilon}$ convergence rate, and thus arrive at a contradiction. A somewhat surprising facet of this proof technique is that the proof has the flavor of an upper bound proof; in particular, we apply Theorem 5 in this argument.

The proof works as follows: Suppose there exists a $\hat{T}_{n}$ such that $|\hat{T}_{n}-T(p,q)|\leq c_{1}n^{-1/2-\epsilon}$ for all $n$. If we are given $2n$ samples, we can use the first half to train KDEs $\hat{p}_{n},\hat{q}_{n}$, and the second half to compute $\hat{T}_{n}$. Armed with these quantities, we can build an estimator for the first and second order terms in the von Mises expansion, which, once $\hat{p}_{n},\hat{q}_{n}$ are fixed, is simply a quadratic functional of the densities. The precise estimator is $\hat{T}_{n}-C_{2}T(\hat{p}_{n},\hat{q}_{n})$. The triangle inequality along with Theorem 5 shows that this estimator converges at rate $n^{-1/2-\epsilon}+n^{\frac{-3s}{2s+d}}$ which is $o(n^{-1/2})$ as soon as $s>d/4$. This contradicts the minimax lower bound for estimating quadratic functionals of Hölder smooth densities. We refer the interested reader to the appendix for details of the proof.

Experiments

We conducted some simulations to examine the empirical rates of convergence of our estimators. We plotted the error as a function of the number of samples $n$ on a $\log$-$\log$ scale in Figure 2 for each estimator and over a number of problem settings. Since our theoretical results are asymptotic in nature, we are not concerned with some discrepancy between the empirical and theoretical rates.

In the top row of Figure 2, we plot the performance of $\widehat{T}_{pl}$ and $\widehat{T}_{lin}$ across four different problem settings: $d=1,s=1$; $d=1,s=2$; $d=2,s=2$; and $d=2,s=4$. The lines fit to the plug-in estimator’s error rate have slopes $-0.25,-0.5,-0.1,-0.2$ from left to right while the lines for the linear estimator have slopes $-0.7,-0.75,-0.65,-0.6$. Qualitatively we see that the $\widehat{T}_{lin}$ is consistently better than $\widehat{T}_{pl}$. We also see that increasing the smoothness $s$ appears to improve the rate of convergence of both estimators.

In the first plot on the bottom row, we record the error rate for $\widehat{T}_{quad}$ with $d=1$ and $s=1.0,2.0$. The fitted lines have slopes $-0.82,-0.93$ respectively, which demonstrate that $\widehat{T}_{quad}$ is indeed a better estimator than $\widehat{T}_{lin}$, at least statistically speaking. Recall that we studied $\widehat{T}_{quad}$ primarily for its theoretical properties and to establish the critical smoothness index of $s>d/4$ for this problem. Computing this estimator is quite demanding, so we did not evaluate it for larger sample size and in higher dimension.

Finally in the last three plots we show the rate of convergence for our divergence estimators, that is $\widehat{T}_{lin}$ plugged into the equations for $D_{\alpha}$ or $T_{\alpha}$ and the quadratic-based estimator for $L_{2}^{2}$. Qualitatively, it is clear that the estimators converge fairly quickly and moreover we can verify that increasing the smoothness $s$ does have some effect on the rate of convergence.

Discussion

In this paper, we address the problem of divergence estimation with corrections of the plug-in estimator. We prove that our estimators enjoy parametric rates of convergence as long as the densities are sufficiently smooth. Moreover, through information theoretic techniques, we show that our best estimator $\widehat{T}_{quad}$ is nearly minimax optimal.

Can we construct divergence estimators that are computationally and statistically efficient? Recall that $\widehat{T}_{quad}$ involves numeric integration and is computationally impractical, yet $\widehat{T}_{lin}$, while statistically inferior, is surprisingly simple when applied to the divergences we consider. At this point we advocate for the use of $\widehat{T}_{lin}$, in spite of its sub-optimality.

What other properties do these estimators enjoy? Can we construct confidence intervals and statistical tests from them? In particular, can we use our estimators to test for independence between two random variables?

Do our techniques yield estimators for other divergences, such as $f$-divergence and the Kullback-Leibler divergence?

Lastly, can one prove a lower bound for the case where $\alpha=\beta=1$, i.e. the $L_{2}$ inner product?

We hope to address these questions in future work.

Acknowledgements

This research is supported by DOE grant DESC0011114, NSF Grants DMS-0806009, IIS1247658, and IIS1250350, and Air Force Grant FA95500910373. AK is supported in part by a NSF Graduate Research Fellowship.

References

Appendix A The von Mises Expansion

Before diving into the auxiliary results of Section 5, let us first derive some properties of the von Mises expansion. It is a simple calculation to verify that the Gateaux derivative is simply the functional derivative of $\phi$ in the event that $T(F)=\int\phi(f)$.

Let $T(F)=\int\phi(f)d\mu$ where $f=dF/d\mu$ is the Radon-Nikodym derivative, $\phi$ is differentiable and let $G$ be some other distribution with density $g=dG/d\mu$. Then:

We now demonstrate that the remainder for the $t$th order von Mises expansion is $O(\|p-\hat{p}\|_{t+1}^{t+1}+\|q-\hat{q}\|_{t+1}^{t+1})$ under the assumption that $p,\hat{p},q,\hat{q}$ are all bounded above and below.

Let $T(p,q)=\int p^{\alpha}q^{\beta}d\mu$ and uppose that $p,\hat{p},q,\hat{q}$ are all bounded from above and below. Then $R_{t}$, the remainder of the $t$th order von Mises expansion of $T(p,q)$ around $T(\hat{p},\hat{q})$ satisfies:

The $t$th order term in the von Mises expansion is:

where $\prod_{i=0}^{0}a_{i}=1$. If we are to take a $t-1$st order expansion, the remainder is of the same form as the $t$th term, except that the terms $\hat{p}^{\alpha-a}(x),\hat{q}^{\beta-(t-a)}(x)$ are replaced by functions $\xi_{1}^{\alpha-a}(x),\xi_{2}^{\beta-(t-a)}(x)$ for some functions $\xi_{1},\xi_{2}$ that are bounded between $p,\hat{p}$ and $q,\hat{q}$ respectively. In our setting, $p,q\in[\kappa_{l},\kappa_{u}]$ and $\hat{p},\hat{q}\in[\kappa_{l}-\epsilon,\kappa_{u}+\epsilon]$ so $\xi_{1},\xi_{2}$ are bounded functions. With this bound, we can simplify the remainder term $R_{t-1}$ to:

Looking at the integral pointwise, either $|p(x)-\hat{p}(x)|\leq|q(x)-\hat{q}(x)|$ in which case the expression is upper bounded by $|q(x)-\hat{q}(x)|^{t}$ or the opposite is true in which case it is bounded by $|p(x)-\hat{p}(x)|^{t}$. Either way, we can upper bound the integral by the sum. This gives:

In many cases, the constant can be worked out:

If $\alpha=\beta=1$, then $R_{1}=\alpha\beta$ while $R_{2},\ldots,=0$.

If $\alpha,\beta>0,\alpha+\beta=1$ as in the Rényi Divergence, $R_{2}=1$ while $R_{3}=\frac{5}{6}\kappa_{\epsilon}^{-2}\alpha\beta$ where $\kappa_{\epsilon}=\min\{\kappa_{l}-\epsilon,(\kappa_{u}+\epsilon)^{-1}\}$.

The second order expansion is computed similarly. The three second order terms are:

Adding these together along with the linear terms, expanding and regrouping terms we get:

Appendix B Full Specification of the Estimators

Here we write out the complete expressions for the estimators $\widehat{T}_{pl},\widehat{T}_{lin},\widehat{T}_{quad}$. Recall that we have samples $X_{1}^{n}\sim p,Y_{1}^{n}\sim q$ and our goal is to estimate $T(p,q)=\int p^{\alpha}q^{\beta}$. Define:

where $DS$ is used to denote that we are data splitting, and $K_{h}$ is a kernel with bandwidth $h$ meeting Assumption 3. The estimator $\widehat{T}_{pl}$ is formed by simply plugging in $\hat{p},\hat{q}$ into the function $T$. Formally:

The estimator $\widehat{T}_{lin}$ is formed by a first order correction but we must used the data split KDEs to ensure independence between the multiple terms in the estimator.

For the quadratic term we perform an additional correction:

Recall that $\{\phi_{k}\}_{k\in D}$ is an orthonormal basis for $L_{2}(^{d})$, and $M$ is an appropriately chosen subset of $D$. The first line of the estimator is simply the plugin term, while the second lines makes up the two linear terms. The third through sixth lines are the two quadratic terms, one involving the data from $p$ and the other involving the data from $q$. Finally the last line is the bilinear term.

Appendix C Detailed Proofs of Upper Bound

Let us now prove the the auxiliary results stated in Section 5

$\int K(x)dx=1$ and $\int\prod_{i}x_{i}^{p_{i}}K(x)dx=0$ for all tuples $p=(p_{1},\ldots,p_{d})$ with $\sum p_{i}\leq\lfloor s\rfloor$.

Note that we can use the Legendre polynomials to construct kernels meeting these properties .

If we expand the expectation and drop the terms that vanish we get all terms of the form:

where the last expression comes from expanding the integral, performing a substitution and bounding $f(x)\leq\kappa_{u}$. So we can bound by $C(q,\kappa_{u},K)h^{-(q-1)d}$. Plugging this into the expression above, we get:

The second inequality holds for $n$ sufficiently large. The third inequality holds whenever $nh^{d}\geq 1$ which will be true for $n$ sufficiently large, given our setting of $h$. To summarize, all of the terms can be upper bounded by $c(n^{p}/h^{pd})$ and there are a constant-in-$p$ number of terms. Plugging this into Equation 15 we get

The constant here has exponential dependence on $p$ but we are only concerned with cases where $p$ is a small constant (at most $4$).

As for the bias (note that $x,u,t$ are all $d$-dimensional vectors here):

Let us define $m=\lfloor s\rfloor$. Taking the $(m-1)$st order von Mises expansion of $f(x+uh)$ about $f(x)$ we get terms of the form:

which are all zero by our assumption on $K$. The remainder term, gives us:

which we will denote $C(m,K,d)Lh^{s}$. Here the function $\xi$ is between $D^{r}f(x)$ and $D^{r}f(x-uh)$ and to reach the last expression, we use the fact that $|D^{r}f(x)-D^{r}f(x-uh)|\leq L\|uh\|^{s-r}$, i.e. the Hölderian assumption on $f$. In applying the Hölderian assumption, there is another term of the form $D^{r}f(x)\int\prod_{i}u_{i}^{r_{i}}K(u)du$ which is zero by the assumption on $K$. Equipped with this bound, we can bound the bias:

This last inequality is an application of Bernstein’s inequality noting that $|\eta_{i}(x)|\leq\frac{2}{h^{d}}\|K\|_{\infty}$ and $\operatorname*{{\rm Var}}(\eta_{i}(x))\leq h^{-d}\kappa_{u}\|K\|_{2}^{2}$ since:

since the second term goes to zero exponentially quickly in $n$. This proves the theorem.

C.2 Convergence Rate for Estimating Linear Functionals

It is trivial to derive the convergence rate for estimating linear functionals:

C.3 Proof of Theorem 6

For the quadratic terms, we use a result of Laurent :

Let $X_{1}^{n}$ be i.i.d random variables with common density $f$ that belongs to some Hilbert Space $L^{2}(d\mu)$. Let $\{\phi_{i}\}_{i\in D}$ be an orthonormal basis of $L^{2}(d\mu)$. Assume that $f$ is uniformly bounded and belongs to the ellipsoid $\mathcal{E}=\{\sum_{i\in D}a_{i}\phi_{i}:\sum_{i\in D}|a_{i}^{2}/c_{i}^{2}|\leq 1\}$. Let $\psi$ be bounded function and define $\theta=\int\psi(x)f(x)\mu(dx)$ and $\hat{\theta}$ as in Equation 2 where the set $M=M_{n}\subset D$ has size $m$. Then whenever $n\geq n_{0}$ (some absolute constant), we have:

For the bi-linear term $\theta_{2,2}^{p,q}$ we have the following theorem:

Let $X_{1}^{n}$ be i.i.d random variables with common density $f$ and $Y_{1}^{n}$ be i.i.d. with common density $g$. Let $f,g$ belong to some Hilbert space $L^{2}(d\mu)$ and let $\{\phi_{i}\}_{i\in D}$ be an orthonormal basis for $L^{2}(d\mu)$. Assume that $f,g$ are uniformly bounded and both belong to the ellipsoid $\mathcal{E}=\{\sum_{i\in D}a_{i}\phi_{i}:\sum_{i\in D}|a_{i}^{2}/c_{i}^{2}|\leq 1\}$. Let $\theta=\int\psi(x)f(x)g(x)\mu(dx)$ and $\hat{\theta}$ be defined by Equation 1 where the set $M=M_{n}\subset D$ has size $m$. Then whenever $n\geq n_{0}$ (some absolute constant), we have:

where $\alpha_{i}=\int\phi_{i}(x)f(x)$ and $\mathcal{P}_{M}f$ is the projection of $f$ onto the subspace defined by $M$. Define $\beta_{i}=\int\phi_{i}(x)g(x)$. If $f,g$ live in the ellipsoid $\mathcal{E}=\{\sum a_{i}\phi_{i}|\sum|a_{i}|^{2}/|c_{i}|^{2}\leq L\}$ then:

The term inside the parenthesis can be bounded as:

so the bias is $\textrm{Bias}^{2}(\hat{\theta})\leq\|\psi\|_{\infty}^{2}L^{2}\sup_{i\notin M}|c_{i}|^{4}$.

As for the variance, let us define $Q(x)$ to be the $m$-dimensional vector of functions $\phi_{i}(x)-\alpha_{i}$ and $R(x)$ to be the $m$-dimensional vector of functions $\phi_{i}(x)\psi(x)-\int\psi\phi_{i}g$. Further define $A,B$ to be the $m$-dimensional vectors with $i$th components $\alpha_{i}=\int\phi_{i}f$ and $\beta_{i}=\int\psi\phi_{i}g$ respectively. Then our estimator can alternatively be written as:

For $T_{2}$ again by independence we have:

Here the last inequality follows from the fact that $\beta_{i}=\int\psi\phi_{i}g$ is the $i$th fourier coefficient of $\psi g$ so $\sum_{i}\beta_{i}\phi_{i}$ is the projection onto $M$. Of course this quantity is bounded by:

Essentially the same argument reveals that $T_{3}$ is bounded in the same way.

In Lemma 14 we show that the class $\mathcal{W}(s^{\prime},L^{\prime})$ contains $\Sigma(s,L)$ as long as $s^{\prime}<s$ and with appropriate choice of $L^{\prime}$. For now let us work in $\mathcal{W}(s^{\prime},L^{\prime})$.

Thinking of $M_{n}$ as an integer lattice with side lengths $m_{0}=m^{1/d}$ we see that $|M_{n}|=m$. Moreover $\sup_{i\notin M_{n}}|c_{i}|^{4}=L^{2}(2/m)^{4s^{\prime}/d}$. For the quadratic terms, this results in the bound:

and plugging in our definition of $m$ followed by some algebraic simplifications, we get

For the bilinear terms, plugging into Theorem 11, we get

Appendix D Proofs of Corollaries 4 and 3

The proof of Corollary 4 is immediate given the decomposition $\|p-q\|_{2}^{2}=\int p^{2}+\int q^{2}-2\int pq$ and the Theorem 6.

For Corollary 3, if we use our estimator $\hat{T}$ for $T(p,q)=\int p^{\alpha}q^{1-\alpha}$ we can plug $\hat{T}$ into the definition of Rényi divergence to obtain an estimator $\hat{D}_{\alpha}$. The rate of convergence is:

where $\gamma$ is the rate of convergence of our estimator. This is $O(n^{-\gamma})$ as long as $T(p,q)\geq c>0$.

Appendix E Detailed Proofs for Lower Bound

To prove the main part of the theorem, the $\Omega(n^{\frac{-4s}{4s+d}})$ rate, we use Le Cam’s method. We decompose the proof into three parts. In the first part, we adapt Le Cam’s method to our setting. In the second part, we show how the properties established on the functions $u_{j}$, $j\in[p]$ allow us to apply the technique and establish the theorem. In the third part, we prove the existence of such functions $u_{j}$. We conclude this section with a proof of the $\Omega(n^{-1/2})$ when $s>d/4$.

Recall that in our proof we partition $^{d}$ into $m$ cubes $R_{1},\ldots,R_{m}$ of side length $m^{-1/d}$. On each bin we require a function $u_{j}$ such that:

where the last inequality needs to hold for all tuples $r$ with $\sum_{j}r_{j}\leq s+1$. Using these functions $u_{j}$, we construct the alternatives $g_{\lambda}=p+K\sum_{\lambda\in\Lambda}\lambda_{j}u_{j}\mathbf{1}_{R_{j}}$ for all $\lambda\in\Lambda=\{-1,1\}^{m}$. The third property above ensures that $g_{\lambda}$ is a valid density.

Properties 2, 4, and 5 ensure that $T(p,q)-T(g_{\lambda},q)$ is sufficiently large. Indeed, by the von Mises expansion:

Here $\xi$ is the function in the Taylor’s remainder theorem, bounded between $p$ and $g_{\lambda}$, both of which are bounded above and below. $g_{\lambda}$ is bounded above and below by property 5 since $\|D_{0}u_{j}\|_{\infty}=\|u_{j}\|_{\infty}\leq 1$ which means that $g_{\lambda}\in[1-K,1+K]$. $K$ will be decreasing with $n$, so this quantity will certainly be bounded for $n$ large enough. Property 2 allows us to arrive at the last line since each $u_{j}$ is orthogonal to the derivative of $T$, so the first term in the expansion is zero. Finally property 4 allows us to lower bound $\|u_{j}\|_{2}^{2}$.

Property 2 is also critical in ensuring that $h^{2}(P^{n}\times Q^{n},\bar{G}^{n}\times Q^{n})$ is small through the following Theorem of Birge and Massart .

Consider a set of densities $p$ and $p_{\lambda}=p[1+\sum_{j}\lambda_{j}v_{j}(x)]$ for $\lambda\in\Lambda=\{-1,1\}^{m}$. Suppose that (i) $\|v_{j}\|_{\infty}\leq 1$ (ii) $\|1_{R_{j}^{C}}v_{j}\|_{1}=0$, (iii) $\int v_{j}p=0$ and (iv) $\int v_{j}^{2}p=\alpha_{j}>0$ all hold with:

Define $\bar{P}_{\Lambda}^{n}=\frac{1}{|\Lambda|}\sum_{\lambda\in\Lambda}P_{\lambda}^{n}$. Then:

where $C<1/3$ is continuous and non-decreasing with respect to each argument and $C(0,0,0)=1/16$.

In bounding the Hellinger distance $h^{2}(P^{n}\times Q^{n},\bar{G}^{n}\times Q^{n})$ we first use the property that hellinger distance decomposes across product measures:

If we define $v_{j}(x)=Ku_{j}(x)/p(x)$ then we have $g_{\lambda}=p[1+\sum_{j}\lambda_{j}v_{j}]$ as needed by Theorem 12. We immediately satisfy requirements 1, 2, and 3 and we have $\int v_{j}^{2}p=K^{2}\int u_{j}^{2}/p\leq K^{2}\kappa_{l}/m=\alpha_{j}$. Thus in applying the theorem we have:

Property 1 and 5 ensure that $g_{\lambda}\in\Sigma(s,L)$ via the following argument. Defining $u_{\lambda}=K\sum_{j}\lambda_{j}u_{j}$, we will first show that $u_{\lambda}$ is holder smooth and $g_{\lambda}$ will be holder by a final application of the triangle inequality. For $u_{\lambda}$, fix $r$ with $\sum_{j}r_{j}=s$ and fix $x,y$. Let $x_{1}$ be the boundary point of $R_{j}$, the bin containing $x$ along the line between $x$ and $y$ and let $y_{1}$ be the analogous boundary point for $y$.

The first line is an application of the triangle inequality. In the second line we use that $u_{\lambda}$ is zero and has all derivatives equal to zero on the boundaries of the cubes $R_{j}$. This follows from the fact that $u_{j}$ is not supported in the band around the border of $R_{j}$. The third line is an application of the fundamental theorem of calculus, $\gamma(x,x_{1})$ is the path between $x$ and $x_{1}$. The fourth line follows from Hölder’s inequality, we replace each derivative with its supremum and are left with just the path integral, which simplifies to the length of the path, i.e. $\|x-x_{1}\|_{2}$. In the fifth line we use the assumption $\|D^{r}u_{j}\|_{\infty}\leq m^{r/d}$ for any derivative operator with $\sum_{j}r_{j}\leq s+1$. To arrive at the sixth line, notice that since $x,x_{1}$ are in the same box $R_{j}$, we have $\|x-x_{1}\|_{2}\leq\sqrt{d}m^{-1/d}$ (there are $m$ boxes and each one has length $m^{-1/d}$ on each side). The last line is true since $x_{1},y_{1}$ are on the line segment between $x,y$.

In other words, $g_{\lambda}$ is holder smooth as long as $Km^{s/d}\sqrt{d}\asymp L$, imposing the requirement that $K=O(m^{-s/d})$. So if we pick $m=n^{\frac{2d}{4s+d}}$ and $K=m^{-s/d}=n^{\frac{-2s}{4s+d}}$ we get that $g_{\lambda}\in\Sigma(s,L)$ as long as there is some wiggle room around $p$. We also get that the Hellinger distance is bounded by $O(n^{2}n^{\frac{-8s}{4s+d}}n^{\frac{-2d}{4s+d}})=O(1)$ and the distance in our metric is $n^{\frac{-4s}{4s+d}}$ as we desired. We can apply Theorem 7 and arrive at the result.

To wrap up, we need to show that we can in fact find the functions $u_{j}$. We can do this by mapping $R_{j}$ to $^{d}$ and using an orthonormal system $\{\phi_{j}\}_{j=1}^{q}$ for $L^{2}(^{d})$ with $q\geq 3$. Suppose that $\phi_{j}$ satisfy (i) $\phi_{1}=1$, $\phi_{j}(x)=0$ for $x\notin[\epsilon,1-\epsilon]^{d}$ and (iii) $\|D^{r}\phi_{j}\|_{\infty}\leq K<\infty$ for all $j$. Certainly we can find such an orthonormal system.