Breaking the Bandwidth Barrier: Geometrical Adaptive Entropy Estimation

Introduction

Unsupervised representation learning is one of the major themes of modern data science; a common theme among the various approaches is to extract maximally “informative" features via information-theoretic metrics (entropy, mutual information and their variations) – the primary reason for the popularity of information theoretic measures is that they are invariant to one-to-one transformations and that they obey natural axioms such as data processing. Such an approach is evident in many applications, as varied as computational biology , sociology and information retrieval , with the citations representing a mere smattering of recent works. Within mainstream machine learning, a systematic effort at unsupervised clustering and hierarchical information extraction is conducted in recent works of . The basic workhorse in all these methods is the computation of mutual information (pairwise and multivariate) from i.i.d. samples. Indeed, sample-efficient estimation of mutual information emerges as the central scientific question of interest in a variety of applications, and is also of fundamental interest to statistics, machine learning and information theory communities.

While these estimation questions have been studied in the past three decades (and summarized in ), the renewed importance of estimating information theoretic measures in a sample-efficient manner is persuasively argued in a recent work , where the authors note that existing estimators perform poorly in several key scenarios of central interest (especially when the high dimensional random variables are strongly related to each other). The most common estimators (featured in scientific software packages) are nonparametric and involve $k$ nearest neighbor (NN) distances between the samples. The widely used estimator of mutual information is the one by Kraskov and Stögbauer and Grassberger and christened the KSG estimator (nomenclature based on the authors, cf. ) – while this estimator works well in practice (and performs much better than other approaches such as those based on kernel density estimation procedures), it still suffers in high dimensions. The basic issue is that the KSG estimator (and the underlying differential entropy estimator based on nearest neighbor distances by Kozachenko and Leonenko (KL) ) does not take advantage of the fact that the samples could lie in a smaller dimensional subspace (more generally, manifold) despite the high dimensionality of the data itself. Such lower dimensional structures effectively act as boundaries, causing the estimator to suffer from what is known as boundary biases.

Ameliorating this deficiency is the central theme of recent works , each of which aims to improve upon the classical KL (differential) entropy estimator of . A local SVD is used to heuristically improve the density estimate at each sample point in , while a local Gaussian density (with empirical mean and covariance weighted by NN distances) is heuristically used for the same purpose in . Both these approaches, while inspired and intuitive, come with no theoretical guarantees (even consistency) and from a practical perspective involve delicate choice of key hyper parameters. An effort towards a systematic study is initiated in which connects the aforementioned heuristic efforts of to the local log-likelihood density estimation methods from theoretical statistics.

The local density estimation method is a strong generalization of the traditional kernel density estimation methods, but requires a delicate normalization which necessitates the solution of certain integral equations (cf. Equation (9) of ). Indeed, such an elaborate numerical effort is one of the key impediments for the entropy estimator of to be practically valuable. A second key impediment is that theoretical guarantees (such as consistency) can only be provided when the bandwidth is chosen globally (leading to poor sample complexity in practice) and consistency requires the bandwidth $h$ to be chosen such that $nh^{d}\to\infty$ and $h\to 0$, where $n$ is the sample size and $d$ is the dimension of the random variable of interest. More generally, it appears that a systematic application of local log-likelihood methods to estimate functionals of the unknown density from i.i.d. samples is missing in the theoretical statistics literature (despite local log-likelihood methods for regression and density estimation being standard textbook fare ). We resolve each of these deficiencies in this paper by undertaking a comprehensive study of estimating the (differential) entropy and mutual information from i.i.d. samples using sample dependent bandwidth choices (typically fixed $k$-NN distances). This effort allows us to connect disparate threads of ideas from seemingly different arenas: NN methods, local log-likelihood methods, asymptotic order statistics and sample-dependent heuristic, but inspired, methods for mutual information estimation suggested in the work of .

Main Results: We make the following contributions.

Density estimation: Parameterizing the log density by a polynomial of degree $p$, we derive simple closed form expressions for the local log-likelihood maximization problem for the cases of $p\leq 2$ for arbitrary dimensions, with Gaussian kernel choices. This derivation, posed as an exercise in [20, Exercise 5.2], significantly improves the computational efficiency upon similar endeavors in the recent efforts of .

Entropy estimation: Using resubstitution of the local density estimate, we derive a simple closed form estimator of the entropy using a sample dependent bandwidth choice (of $k$-NN distance, where $k$ is a fixed small integer independent of the sample size): this estimator outperforms state of the art entropy estimators in a variety of settings. Since the bandwidth is data dependent and vanishes too fast (because $k$ is fixed), the estimator has a bias, which we derive a closed form expression for and show that it is independent of the underlying distribution and hence can be easily corrected: this is our main theoretical contribution, and involves new theorems on asymptotic statistics of nearest neighbors generalizing classical work in probability theory , which might be of independent mathematical interest.

Generalized view: We show that seemingly very different approaches to entropy estimation – recent works of and the classical work of fixed $k$-NN estimator of Kozachenko and Leonenko – can all be cast in the local log-likelihood framework as specific kernel and sample dependent bandwidth choices. This allows for a unified view, which we theoretically justify by showing that resubstitution entropy estimation for any kernel choice using fixed $k$-NN distances as bandwidth involves a bias term that is independent of the underlying distribution (but depends on the specific choice of kernel and parametric density family). Thus our work is a strict mathematical generalization of the classical work of .

Mutual Information estimation: The inspired work of constructs a mutual information estimator that subtly altered (in a sample dependent way) the three KL entropy estimation terms, leading to superior empirical performance. We show that the underlying idea behind this change can be incorporated in our framework as well, leading to a novel mutual information estimator that combines the two ideas and outperforms state of the art estimators in a variety of settings.

In the rest of this paper we describe these main results, the sections organized in roughly the same order as the enumerated list.

Local likelihood density estimation (LLDE)

where maximization is over an exponential polynomial family, locally approximating $f(u)$ near $x$:

For higher degree local likelihood, we provide simple closed form solutions and provide a proof in Section 8.1. Somewhat surprisingly, this result has eluded prior works and which specifically attempted the evaluation for $p=2$. Part of the subtlety in the result is to critically use the fact that the parametric family (eg., the polynomial family in (2)) need not be normalized themselves; the local log-likelihood maximization ensures that the resulting density estimate is correctly normalized so that it integrates to 1.

[20, Exercise 5.2] For a degree $p\in\{1,2\}$, the maximizer of local likelihood (1) admits a closed form solution, when using the Gaussian kernel $K(u)=e^{-\frac{\|u\|^{2}}{2}}$. In case of $p=1$,

In case of $p=2$, for $S_{0}$ and $S_{1}$ defined as above,

where it follows from Cauchy-Schwarz that $\Sigma$ is positive semidefinite.

One of the major drawbacks of the KDE and $k$-NN methods is the increased bias near the boundaries. LLDE provides a principled approach to automatically correct for the boundary bias, which takes effect only for $p\geq 2$ . This explains the performance improvement for $p=2$ in the figure below (left panel), and the gap increases with the correlation as boundary effect becomes more prominent. We use the proposed estimators with $p\in\{0,1,2\}$ to estimate the mutual information between two jointly Gaussian random variables with correlation $r$, from $n=500$ samples, using resubstitution methods explained in the next sections. Each point is averaged over $100$ instances.

In the right panel, we generate i.i.d. samples from a 2-dimensional Gaussian with correlation 0.9, and found local approximation $\widehat{f}(u-x^{*})$ around $x^{*}$ denoted by the blue $*$ in the center. Standard $k$-NN approach fits a uniform distribution over a circle enclosing $k=20$ nearest neighbors (red circle). The green lines are the contours of the degree-2 polynomial approximation with bandwidth $h=\rho_{20,x}$. The figure illustrates that $k$-NN method suffers from boundary effect, where it underestimates the probability by over estimating the volume in (4). However, degree-2 LDDE is able to correctly capture the local structure of the pdf, correcting for boundary biases.

Despite the advantages of the LLDE, it requires the bandwidth to be data independent and vanishingly small (sublinearly in sample size) for consistency almost everywhere – both of these are impediments to practical use since there is no obvious systematic way of choosing these hyperparameters. On the other hand, if we restrict our focus to functionals of the density, then both these issues are resolved: this is the focus of the next section where we show that the bandwidth can be chosen to be based on fixed $k$-NN distances and the resulting universal bias easily corrected.

k𝑘k-LNN Entropy Estimator

We consider resubstitution entropy estimators of the form ${\widehat{H}}(x)=-(1/n)\sum_{i=1}^{n}\log\widehat{f}_{n}(X_{i})$ and propose to use the local likelihood density estimator in (7) and a choice of bandwidth that is local (varying for each point $x$) and adaptive (based on the data). Concretely, we choose, for each sample point $X_{i}$, the bandwidth $h_{X_{i}}$ to be the the distance to its $k$-th nearest neighbor $\rho_{k,i}$. Precisely, we propose the following $k$-Local Nearest Neighbor ($k$-LNN) entropy estimator of degree-$2$:

The truncation is important for computational efficiency, but the analysis works as long as $m=O(n^{{1/(2d)}-\varepsilon})$ for any positive $\varepsilon$ that can be arbitrarily small. For a larger $m$, for example of $\Omega(n)$, those neighbors that are further away have a different asymptotic behavior. We show in Theorem 1 that the asymptotic bias is independent of the underlying distribution and hence can be precomputed and removed, under mild conditions on a twice continuously differentiable pdf $f(x)$ (cf. Lemma 3.1 below).

This proves the $L_{1}$ and $L_{2}$ consistency of the $k$-LNN estimator; we relegate the proof to Section 10 for ease of reading the main part of the paper. The proof assumes Ansatz 1 (also stated in Section 10, which states that a certain exchange of limit holds. As noted in , such an assumption is common in the literature on consistency of $k$-NN estimators, where it has been implicitly assumed in existing analyses of entropy estimators including , without explicitly stating that such assumptions are being made. Our choice of a local adaptive bandwidth $h_{X_{i}}=\rho_{k,i}$ is crucial in ensuring that the asymptotic bias $B_{k,d}$ does not depend on the underlying distribution $f(x)$. This relies on a fundamental connection to the theory of asymptotic order statistics made precise in Lemma 3.1, which also gives the explicit formula for the bias below.

In practice, we propose using a fixed small $k$ such as five. For $k\leq 3$ the estimator has a very large variance, and numerical evaluation of the corresponding bias also converges slowly. For some typical choices of $k$, we provide approximate evaluations below, where $0.0183(\pm 6)$ indicates empirical mean $\mu=183\times 10^{-4}$ with confidence interval $6\times 10^{-4}$. In these numerical evaluations, we truncated the summation at $m=50,000$. Although we prove that $B_{k,d}$ converges in $m$, in practice, one can choose $m$ based on the number of samples and $B_{k,d}$ can be evaluated for that $m$.

Empirical contribution: Numerical experiments suggest that the proposed estimator outperforms state-of-the-art entropy estimators, and the gap increases with correlation. The idea of using $k$-NN distance as bandwidth for entropy estimation was originally proposed by Kozachenko and Leonenko in , and is a special case of the $k$-LNN method we propose with degree and a step kernel. We refer to Section 4 for a formal comparison. Another popular resubstitution entropy estimator is to use KDE in (3) , which is a special case of the $k$-LNN method with degree , and the Gaussian kernel is used in simulations. As comparison, we also study a new estimator based on von Mises expansion (as opposed to simple re-substitution) which has an improved convergence rate in the large sample regime. In Figure 2 (left), we draw $100$ samples i.i.d. from two standard Gaussian random variables with correlation $r$, and plot resulting mean squared error averaged over $100$ instances. The ground truth, in this case is $H(X)=\log(2\pi e)+0.5\log(1-r^{2})$. On the right, we repeat the same simulation for fixed $r=0.99999$ and varying number of samples and $m=7\log_{e}n$.

In Figure 3, we repeat the same simulation for 6 standard Gaussian random variables with ${\rm Cov}(X_{1},X_{2})={\rm Cov}(X_{3},X_{4})={\rm Cov}(X_{5},X_{6})=r$ and ${\rm Cov}(X_{i},X_{j})=0$ for other pairs $(i,j)$. On the left, we draw $100$ i.i.d. samples with various $r$. We plot resulting mean squared error averaged over $100$ instances. The ground truth is $H(X)=3\log(2\pi e)+1.5\log(1-r^{2})$. On the right, we repeat the same simulation for fixed $r=0.99999$ and varying number of samples and $m=7\log_{e}n$.

In Figure 4 (left), we draw $100$ samples i.i.d. from a mixture of two joint Gaussian distributions with zero mean and covariance $\begin{pmatrix}1&r\\ r&1\end{pmatrix}$ and $\begin{pmatrix}1&-r\\ -r&1\end{pmatrix}$, respectively, and plot resulting average estimate over $100$ instances. Here we plot an upper bound of the ground truth $H(X)\leq\log(2)+\log(2\pi e)+0.5\log(1-r^{2})$ for $r\geq 0.9$. On the right, we repeat the same simulation for fixed $r=0.99999$ and varying number of samples and $m=7\log_{e}n$.

Universality of the k𝑘k-LNN approach

for some constant ${\widetilde{B}}_{k,d,p,K}$ that only depends on $k,p,K$ and $d$.

We provide a proof in Section 11. Although in general there is no simple analytical characterization of the asymptotic bias ${\widetilde{B}}_{k,p,K,d}$ it can be readily numerically computed: since ${\widetilde{B}}_{k,p,K,d}$ is independent of the underlying distribution, one can run the estimator over i.i.d. samples from any distribution and numerically approximate the bias for any choice of the parameters. However, when the maximization $\widehat{a}(x)=\arg\max_{a}{\cal L}_{x}(f_{a,x})$ admits a closed form solution, as is the case with proposed $k$-LNN, then $\widetilde{B}_{k,p,K,d}$ can be characterized explicitly in terms of uniform order statistics.

k𝑘k-LNN Mutual information estimator

Given an entropy estimator ${\widehat{H}}_{\rm KL}$, mutual information can be estimated: ${\widehat{I}}_{\rm 3KL}={\widehat{H}}_{\rm KL}(X)+{\widehat{H}}_{\rm KL}(Y)-{\widehat{H}}_{\rm KL}(X,Y)$. In , Kraskov and Stögbauer and Grassberger introduced ${\widehat{I}}_{\rm KSG}(X;Y)$ by coupling the choices of the bandwidths. The joint entropy is estimated in the usual way, but for the marginal entropy, instead of using $k$NN distances from $\{X_{j}\}$, the bandwidth $h_{X_{i}}=\rho_{k,i}(X,Y)$ is chosen, which is the $k$ nearest neighbor distance from $(X_{i},Y_{i})$ for the joint data $\{(X_{j},Y_{j})\}$. Consider ${\widehat{I}}_{\rm 3LNN}(X;Y)={\widehat{H}}_{k{\rm LNN}}(X)+{\widehat{H}}_{k{\rm LNN}}(Y)-{\widehat{H}}_{k{\rm LNN}}(X,Y)$. Inspired by , we introduce the following novel mutual information estimator we denote by ${\widehat{I}}_{\rm LNN-KSG}(X;Y)$. where for the joint $(X,Y)$ we use the LNN entropy estimator we proposed in (9), and for the marginal entropy we use the bandwidth $h_{X_{i}}=\rho_{k,i}(X,Y)$ coupled to the joint estimator. Empirically, we observe ${\widehat{I}}_{\rm KSG}$ outperforms ${\widehat{I}}_{\rm 3KL}$ everywhere, validating the use of correlated bandwidths. However, the performance of ${\widehat{I}}_{{\rm LNN-KSG}}$ is similar to ${\widehat{I}}_{3{\rm LNN}}$–sometimes better and sometimes worse.

In Figure 5 (left), we estimate mutual information under the same setting as in Figure 2 (left). For most regimes of correlation $r$, both 3LNN and LNN-KSG outperforms other state-of-the-art estimators. The gap increases with correlation $r$. On the right, we draw i.i.d. samples from two random variables $X$ and $Y$, where $X$ is uniform over $$and$Y=X+U$, where$U$is uniform over$[0,0.01]$independent of$X$. In the large sample limit, all estimators find the correct mutual information. The plot show how sensitive the estimates are, in the small sample regime. Both LNN and LNN-KSG are significantly more robust compared to other approaches. Mutual information estimators have been recently proposed in based on local likelihood maximization. However, they involve heuristic choices of hyper-parameters or solving elaborate optimization and numerical integrations, which are far from being easy to implement.

In Figure 6, we test the mutual information estimators for $Y=f(X)+U$, where $X$ is uniformly distributed over $$and$U$is uniformly distributed over$[0,\theta]$, independent of$X$, for some noise level$\theta$. Similar simulation were studied in . We draw 2500 i.i.d. sample points for each relationship. The plot show that for small noise level$\theta$, i.e., near-functional related random variables, our proposed estimators${\widehat{I}}_{3LNN}$and${\widehat{I}}_{LNN-KSG}$ perform much better than 3KL and KSG estimators. Also our proposed estimators can handle both linear and nonlinear functional relationships.

In Figure 7, we test our estimators on linear and nonlinear relationships for both low-dimensional ($D=2$) and high-dimensional ($D=5$). Here $X_{i}$’s are uniformly distributed over $$and$U$is uniformly distributed over$[-3^{8}/2,3^{8}/2]$, independently of$X_{i}$’s. Similar simulation were studied in . We can see that our estimators${\widehat{I}}_{3LNN}$and${\widehat{I}}_{LNN-KSG}$converges much faster than$\hat{I}_{3KL}$and${\widehat{I}}_{KSG}$.

Breaking the bandwidth barrier

While $k$-NN distance based bandwidth are routine in practical usage , the main finding of this work is that they also turn out to be the “correct" mathematical choice for the purpose of asymptotically unbiased estimation of an integral functional such as the entropy: $-\int f(x)\log f(x)$; we briefly discuss the ramifications below. Traditionally, when the goal is to estimate $f(x)$, it is well known that the bandwidth should satisfy $h\to 0$ and $nh^{d}\to\infty$, for KDEs to be consistent. As a rule of thumb, $h=1.06\widehat{\sigma}n^{-1/5}$ is suggested when $d=1$ where $\widehat{\sigma}$ is the sample standard deviation [41, Chapter 6.3]. On the other hand, when estimating entropy, as well as other integral functionals, it is known that resubstitution estimators of the form $-(1/n)\sum_{i=1}^{n}\log\widehat{f}(X_{i})$ achieve variances scaling as $O(1/n)$ independent of the bandwidth . This allows for a bandwidth as small as $O(n^{-1/d})$.

The bottleneck in choosing such a small bandwidth is the bias, scaling as $O(h^{2}+(nh^{d})^{-1}+E_{n})$ , where the lower order dependence on $n$, dubbed $E_{n}$, is generally not known. The barrier in choosing a global bandwidth of $h=O(n^{-1/d})$ is the strictly positive bias whose value depends on the unknown distribution and cannot be subtracted off. However, perhaps surprisingly, the proposed local and adaptive choice of the $k$-NN distance admits an asymptotic bias that is independent of the unknown underlying distribution. Manually subtracting off the non-vanishing bias gives an asymptotically unbiased estimator, with a potentially faster convergence as numerically compared below. Figure 8 illustrates how $k$-NN based bandwidth significantly improves upon, say a rule-of-thumb choice of $O(n^{-1/(d+4)})$ explained above and another choice of $O(n^{-1/(d+2)})$. In the left figure, we use the setting from Figure 2 (right) but with correlation $r=0.999$. On the right, we generate $X\sim{\cal N}(0,1)$ and $U$ from uniform $[0,0.01]$ and let $Y=X+U$ and estimate $I(X;Y)$. Following recent advances in , the proposed local estimator has a potential to be extended to, for example, Renyi entropy, but with a multiplicative bias as opposed to additive.

Discussion

The topic of estimation of an integral functional of an unknown density from i.i.d. samples is a classical one in statistics and we tie together a few pertinent topics from the literature in the context of the results of this manuscript.

Only the convergence analysis of the distances, and not the directions, is required for traditional $k$-NN based estimators, such as the entropy estimator of . In the seminal paper, introduced resubstitution entropy estimators of the form ${\widehat{H}}(X)=-(1/n)\sum_{i=1}^{n}\log\widehat{f}_{n}(X_{i})$ with $\widehat{f}_{n}(x)={k}/(n\,C_{d}\,\rho_{k,x}^{d})$ (as defined in (4)). This $k$-NN estimator has a non-vanishing asymptotic bias, which was computed as $B_{k,d}=(\psi(k)-\log(k))$ with the digamma function $\psi(\cdot)$ and was suggested to be manually removed. For $k=1$ this was proved in the original paper of , which later was extended in to general $k$. This mysterious bias term $B_{k,d}=(\psi(k)-\log(k))$ whose original proofs in provided little explanation for, can be alternatively proved with both rigor and intuition by making connections to uniform order statistics. For a special case of $k=1$, with extra assumptions on the support being compact, such an elegant proof is provided in [2, Theorem 7.1] which explicitly applies the convergence of the nearest neighbor distance to uniform order statistics. Namely,

2 Convergence rate of the bias

Establishing the convergence rate of the KL estimator is a challenging problem, and is not quite resolved despite work over the past three decades. The $O(1/n)$ convergence rate of the variance is established in under various assumptions. Establishing the convergence rate of the bias is more challenging. It has been first studied in , where root-$n$ consistency is shown in 1-dimension with bounded support and assuming $f(x)$ is bounded below. is the first to prove a root mean squared error convergence rate of $O(1/\sqrt{n})$ for general densities with unbounded support in 1-dimension and exponentially decaying tail, such as the Gaussian density. These assumptions are relaxed in , where zeroes and fat tails are allowed in $f(x)$. In general $d$-dimensions, prove bounds on the convergence rate of the bias for finite $k=O(1)$, and for $k=\Omega(\log n)$. Establishing the convergence rate for the bias of the proposed local estimator is an interesting open problem – it is interesting to see if the superior empirical performance of the local estimator is captured in the asymptotics of rate of convergence of the bias.

It is intuitive that kernel density estimators can capture the structure in the distribution if the distribution lies on a lower dimensional manifold. This is made precise in , which also shows improved convergence rates for distributions whose support is on low dimensional manifolds. However, the estimator in critically uses the geodesic distances between the sample points on the manifold. Given that the proposed estimators fit distributions locally, a concrete question of interest is whether such an improvement can be achieved without such an explicit knowledge of the geodesic distances, i.e., whether the local estimators automatically adapt to underlying lower dimensional structures.

3 Ensemble estimators

Recent works have proposed ensemble estimators, which use known estimators based on kernel density estimators and $k$-NN methods and construct a new estimate by taking the weighted linear combination of those methods with varying bandwidth or $k$, respectively. With a proper choice of the weights, which can be computed analytically by solving a simple linear program, a boosting of the convergence rate can be achieved. The key property that allows the design of such ensemble estimators is that the leading terms (in terms of the sample size $n$) of the bias have a multiplicative constant that only depends on the unknown distribution. An intuitive explanation for this phenomenon is provided in in the context of $k$-NN methods; it is interesting to explore if such a phenomenon continues in the $k$-LNN scenario studied in this paper. Such a study would potentially lead to ensemble-based estimators in the local setting and also naturally allow a careful understanding of the rate of convergence of the bias term.

Proofs

We first prove the derivation of the LLDE with degree $p=2$ in Equation (7). The gradient of the local likelihood evaluated at the maximizer is zero , which gives a computational tool for finding the maximizer:

where $K(x)=\exp\{-\|x\|^{2}/2\}$ is the Gaussian kernel. Notice that the left-hand side of the equations are $S_{0}/n$, $S_{1}/n$ and $S_{2}/n$, respectively. The RHS can be written in closed forms as:

where $M=h^{-2}I_{d\times d}-2a_{2}$ assuming $h$ sufficiently small such that $M$ is positive definite. We want to derive $\hat{f}(x)=\exp\{a_{0}\}$ from the equations. From (20) we get $M^{-1}a_{1}=S_{1}(h/S_{0})$. Together with (21), we get $M^{-1}+M^{-1}a_{1}a_{1}^{T}M^{-1}=S_{2}(h^{2}/S_{0})$. Hence, $M^{-1}=(S_{2}/S_{0}-(S_{1}/S_{0})(S_{1}/S_{0})^{T})h^{2}=h^{2}\Sigma$. Plug them in (19), we obtain the desired expression.

Analogously, for the derivation of the LLDE with degree $p=1$ in Equation (5), we get

This gives $a_{1}=(1/(hS_{0}))S_{1}$, and $e^{a_{0}}=(S_{0}/(n(2\pi)^{d/2}h^{d}))\,\exp\{-0.5\|S_{1}\|^{2}/S_{0}^{2}\}$.

Proof of Lemma 3.1

Now consider the first term in (25). We consider two cases separately.

Case 2. If $\|Z_{m,i}\|<(\sqrt{n}c_{d}f(x))^{-1/d}$, let $\overline{B}=\{t:(c_{d}nf(x))^{1/d}t\in B\textrm{ and }\|t_{m}\|<(\sqrt{n}c_{d}f(x))^{-1/d}\}$ and $\overline{B_{\theta}}=\{t:(c_{d}nf(x))^{1/d}t\in B_{\theta}\textrm{ and }t_{m}<(\sqrt{n}c_{d}f(x))^{-1/d}\}$. Note that

where the first inequality follows from the fact that $\int_{\theta\in(S^{d-1})^{m}}(\int_{\overline{B_{\theta}}}g(t_{m})dt)d(\sigma^{d-1})^{m}(\theta)=\int_{\overline{B}}g(\|t_{m}\|)dt$. Since $f$ is continuously differentiable, by mean value theorem, there exists $a,b\in B(x,(\sqrt{n}c_{d}f(x))^{-1/d})$ such that

By the assumption, there exists a ball $B(x,\varepsilon)$ such that $\|\nabla f(a)\|=O(1)$ and $f(a)>0$ for all $a\in B(x,\varepsilon)$, so for sufficiently large $n$ such that $(\sqrt{n}c_{d}f(x))^{-1/d}<\varepsilon$, there exists some constant $C$ such that $\sup_{\|t\|\leq(\sqrt{n}c_{d}f(x))^{-1/d}}f(x+t)\leq(1+Cn^{-1/(2d)})\inf_{\|t\|\leq(\sqrt{n}c_{d}f(x))^{-1/d}}f(x+t)$. Therefore, (28) is upper bounded by $(1+Cn^{-1/(2d)})^{m}$. Similarly, (28) is lower bounded by $(1-Cn^{-1/(2d)})^{m}$.

For simplicity, let $\mathcal{E}=\{\|Z_{m,i}\|<(\sqrt{n}c_{d}f(x))^{-1/d}\}$. Then combining the two cases, the first term in (25) is bounded by:

Now consider the second term of (25). We will use Corollary 5.5.5 of to show that this term vanishes for $m=O(\log n)$ and as $n$ grows.

Let $Y_{1},Y_{2},\dots,Y_{n}$ be i.i.d. samples from unknown distribution with pdf $f$. Let $Y_{1:n}\leq Y_{2:n}\leq\dots\leq Y_{n:n}$ be the order statistics. Assume the density $f$ satisfies $|\log f(y)|\leq Ly^{\delta}$ for $0<y<y_{0}$ and $f(y)=0$ for $y<0$, where $L$ and $\delta$ are constants. Then

where $C_{0}>0$ is a constant. $E_{1},\dots,E_{m}$ are i.i.d standard exponential random variables.

where $r_{t}=(t/(c_{d}f(x)))^{1/d}$. Here we have:

If $f$ is twice continuously differentiable, we have:

where $S^{d-1}$ is the $(d-1)$-sphere centered at $x$ with radius $r_{t}$ and $\sigma^{d-1}$ is the spherical measure. By mean value theorem, there exists $a(y)\in B(x,r_{t})$ such that $f(y)-f(x)=(y-x)^{T}\nabla f(x)+(a(y)-x)^{T}H_{f}(a(y))(a(y)-x)$, where $a(y)$ depends on $y$. Therefore,

Since there exists a ball $B(x,\varepsilon)$ such that $\|H_{f}(a)\|=O(1)$ for all $a\in B(x,\varepsilon)$. Therefore, for sufficiently small $t$ such that $r_{t}<\varepsilon$, we have:

Therefore, by combing (30) and (37), we have:

Proof of Theorem 1

We first compute the asymptotic bias. We define new notations to represent the estimate as

Let $H_{i}\equiv h((c_{d}nf(X_{i}))^{1/d}Z_{k,i},S_{0,i},S_{1,i},S_{2,i}))-\log f(X_{i})$. Since the terms $H_{1},H_{2},\dots,H_{n}$ are identically distributed, the expected value of ${\widehat{H}}_{k}^{(n)}$ converges to

Under this ansatz, perhaps surprisingly, we will show that the expectation inside converges to $-\log f(X_{1})$ plus some bias that is independent of the underlying distribution. Precisely, for almost every $x$ and given $X_{1}=x$,

as $n\to\infty$ where $B_{k,d}$ is a constant that only depends on $k$ and $d$, defined in (44). This implies that

Together with (40), this finishes the proof of the desired claim.

We are now left to prove the convergence of (42). We first give a formal definition of the bias $B_{k,d}$ by replacing the sample defined quantities by a similar quantities defined from order-statistics, and use Lemma 3.1 to prove the convergence. Recall that our order-statistics is defined by two sequences of $m$ i.i.d. random variables: i.i.d. standard exponential random variables $E_{1},\dots,E_{m}$ and i.i.d. random variables $\xi_{1},\dots,\xi_{m}$ uniformly distributed over $S^{d-1}$. We define

and we show that both terms converge to zero for any $m=\Theta(\log n)$. Given that $h$ is continuous and bounded, this implies that

where the last inequality follows from Lemma 3.1. By the assumption that $f$ has open support and $\|\nabla f\|$ and $\|H_{f}\|$ is bounded almost everywhere, this convergence holds for almost every $x$.

Assume $m\to\infty$ as $n\to\infty$ and $k\geq 3$ , then

Combine (48) and (50) in (46), this implies the desired claim.

We next prove the upper bound on the variance, following the technique from [2, Section 7.3]. For the usage of Efron-Stein inequality, we need a second set of i.i.d. samples $\{X^{\prime}_{1},X^{\prime}_{2},\dots,X^{\prime}_{n}\}$. For simplicity, denote ${\widehat{H}}={\widehat{H}}_{kLNN}^{(n)}(X)$ be the kLNN estimate base on original sample $\{X_{1},\dots,X_{n}\}$ and ${\widehat{H}}^{(i)}$ be the kLNN estimate based on $\{X_{1},\dots,X_{i-1},X^{\prime}_{i},X_{i+1},\dots X_{n}\}$. Then Efron-Stein theorem states that

Similarly, we can write ${\widehat{H}}^{(j)}=\frac{1}{n}\sum_{i=1}^{n}H_{i}^{(j)}$ for any $j\in\{1,\dots,n\}$. Therefore, the difference of ${\widehat{H}}$ and ${\widehat{H}}^{(j)}$ can be bounded by:

Take expectation over $X_{1}$, we obtain:

where the last inequality comes from the assumption that $\int f(x)(\log f(x))^{2}dx<+\infty$. Combining with (51) and (54), we have

for $t\in[1,+\infty)$. Given the CDF of $R_{j}$, each term in (66) is upper bounded by:

Therefore, in order to establish an upper bound for (66), we need an upper bound for $F_{R_{j}}(t)$. Here we will consider two cases depending on $t$. If $t>(j/2k)^{1/d}$, we just use the trivial upper bound $F_{R_{j}}(t)<1$. If $1\leq t\leq(j/2k)^{1/d}$, since $t^{d}\geq 1$, we have:

here we use the fact that $t\leq(j/2k)^{1/d}$ so $j-t^{d}k>j/2$. Therefore, $F_{R_{j}}(t)\leq 4t^{2d}k/j^{2}$ for $t>(j/2k)^{1/d}$. Combine the two cases and plug into (61), we obtain:

where $C_{d}=\int_{t=1}^{\infty}t^{2d+3}e^{-t^{2}}dt$ is a constant only depend on $d$. Therefore, we can see that

Proof of Theorem 2

The proposed estimator is a solution to a maximization problem $\widehat{a}=\arg\max_{a}{\cal L}_{X_{i}}(f_{a,X_{i}})$. From we know that the maximizer is a fixed point of a series of non-linear equations of the form

Acknowledgement

This work is supported by NSF SaTC award CNS-1527754, NSF CISE award CCF-1553452, NSF CISE award CCF-1617745. We thank the anonymous reviewers for their constructive feedback.

References