Non-Parametric Inference Adaptive to Intrinsic Dimension

TL;DR

This paper introduces a non-parametric method for high-dimensional conditional inference that adapts to the intrinsic low-dimensional structure of data, enabling accurate estimation even when the ambient dimension exceeds the sample size.

Contribution

It develops an adaptive sub-sampled $k$-NN estimator that achieves optimal rates based on intrinsic dimension, with a data-driven method for tuning parameters.

Findings

Estimation error scales as $n^{-1/(d+2)}$ with intrinsic dimension $d$.

Estimator is asymptotically normal regardless of ambient dimension.

Proposed method adapts to unknown intrinsic dimension, improving high-dimensional inference.

Abstract

We consider non-parametric estimation and inference of conditional moment models in high dimensions. We show that even when the dimension $D$ of the conditioning variable is larger than the sample size $n$, estimation and inference is feasible as long as the distribution of the conditioning variable has small intrinsic dimension $d$, as measured by locally low doubling measures. Our estimation is based on a sub-sampled ensemble of the $k$-nearest neighbors ($k$-NN) $Z$-estimator. We show that if the intrinsic dimension of the covariate distribution is equal to $d$, then the finite sample estimation error of our estimator is of order $n^{-1/(d+2)}$ and our estimate is $n^{1/(d+2)}$-asymptotically normal, irrespective of $D$. The sub-sampling size required for achieving these results depends on the unknown intrinsic dimension $d$. We propose an adaptive data-driven approach for choosing this parameter and prove that it achieves the desired rates. We discuss extensions and applications to heterogeneous treatment effect estimation.

Figures

Equations

$$\mathbf{E}[\psi (Z;\theta (x))|X=x]=0,$$

$$\hat{\theta }(x)\text{ solves}:\sum _{i=1}^{n}K(x,X_i,S)\psi (Z_i;\theta )=0,$$

$$ϵ(s)=\mathbf{E}\left[\sup \{\Vert X_i-x{\Vert }_2:i\in S_b,K(x,X_i,S_b)>0,|S_b|=s\}\right].$$

$$\eta (s)=\mathbf{E}\left[\mathbf{E}{\left[K(x,X_i,S_b)|X_i\right]}^2\right].$$

$$\mathrm{Var}(\hat{\theta }(x))=\frac{{\sigma }^2(x)s^2}{n(2s-1)}.$$

$$\theta (x)\text{ solves: }m(x;\theta )=\mathbf{E}[\psi (Z;\theta )\mid X=x]=0.$$

$$\Psi (x;\theta )=\frac{1}{B}\sum _{b=1}^B{\Psi }_{S_b,{\omega }_b}(x;\theta )=\frac{1}{B}\sum _{b=1}^B\sum _{i\in S_b}{\alpha }_{S_b,{\omega }_b}(X_i)\psi (Z_i;\theta ).$$

$${\Psi }_0(x;\theta )={\left(\genfrac{}{}{0ex}{}{n}{s}\right)}^{-1}\sum _{S\in [n]:|S|=s}{\mathbf{E}}_{\omega }\left[\sum _{i\in S}{\alpha }_{S,\omega }(X_i)\psi (Z_i;\theta )\right].$$

$$\alpha (X_{(i)})=\{\begin{array}{ll}\frac{1}{k}{\left(\genfrac{}{}{0ex}{}{n}{s}\right)}^{-1}\left(\genfrac{}{}{0ex}{}{n-i}{s-1}\right)& \text{ if }i\le k\\ \frac{1}{k}{\left(\genfrac{}{}{0ex}{}{n}{s}\right)}^{-1}{\sum }_{j=0}^{k-1}\left(\genfrac{}{}{0ex}{}{i-1}{j}\right)\left(\genfrac{}{}{0ex}{}{n-i}{s-1-j}\right)& \text{ if }i\ge k+1.\end{array}$$

$$\mathbf{E}\left[\sup \left\{\Vert x-X_i{\Vert }_2:K(x,X_i,S)>0\right\}\right]=ϵ(s).$$

$$\sup \left\{\Vert x-X_i{\Vert }_2:K(x,X_i,S)>0\right\}\le ϵ(s,\delta ).$$

$$\Vert \mathrm{Var}\left(\psi (Z;\theta )\mid X=x\right)-\mathrm{Var}\left(\psi (Z;\theta )\mid X=x^{\prime }\right){\Vert }_F\le L_{mm}\Vert x-x^{\prime }{\Vert }_2.$$

$$\Vert \hat{\theta }-\theta (x){\Vert }_2\le \frac{2}{\lambda }\left(L_mϵ(s)+O\left({\psi }_{\max }\sqrt{\frac{ps}{n}\left(\log \log (n/s)+\log (p/\delta )\right)}\right)\right).$$

$$\sqrt{\mathbf{E}\left[\Vert \hat{\theta }-\theta (x){\Vert }_2^2\right]}\le \frac{2}{\lambda }\left(L_mϵ(s)+O\left({\psi }_{\max }\sqrt{\frac{ps}{n}\log \log (pn/s)}\right)\right).$$

$${\sigma }_{n,\beta }^2(x)=\frac{s^2}{n}\mathrm{Var}\left[\mathbf{E}\left[\sum _{i=1}^{s}K(x,X_i,\{Z_j{\}}_{j=1}^s)⟨\beta ,M_0^{-1}\psi (Z_i;\theta (x))⟩\mid Z_1\right]\right].$$

$$\max \left(ϵ(s),ϵ(s{)}^{1/4}{\left(\frac{s}{n}\log \log (n/s)\right)}^{1/2},{\left(\frac{s}{n}\log \log (n/s)\right)}^{5/8}\right)=o({\sigma }_{n,\beta }(x)).$$

$$\frac{⟨\beta ,\hat{\theta }-\theta (x)⟩}{{\sigma }_{n,\beta }(x)}{\to }_d\mathsf{N}(0,1).$$

$${\eta }_s=\mathbf{E}\left[\mathbf{E}{\left[K(x,X_1,\{Z_j{\}}_{j=1}^s)\mid X_1\right]}^2\right]\ge \frac{{\left(1-\delta \right)}^2(1/s{)}^2}{{\inf }_{\rho >0}\left(\mu (B(x,ϵ(s,\rho )))+\rho s/\delta \right)}.$$

$$\Vert x-X_{(k)}{\Vert }_2\le {ϵ}_k(s,\delta )=O{\left(\frac{\log (1/\delta )}{s}\right)}^{1/d}.$$

$$\Vert \hat{\theta }-\theta (x){\Vert }_2\le \frac{2}{\lambda }\left(O\left(s^{-1/d}\right)+O\left({\psi }_{\max }\sqrt{\frac{ps}{n}\left(\log \log (n/s)+\log (p/\delta )\right)}\right)\right),$$

$$\sqrt{\mathbf{E}\left[\Vert \hat{\theta }-\theta (x){\Vert }_2^2\right]}\le \frac{2}{\lambda }\left(O\left(s^{-1/d}\right)+O\left({\psi }_{\max }\sqrt{\frac{sp\log \log (pn/s)}{n}}\right)\right).$$

$${\eta }_k(s)=\mathbf{E}\left[\mathbf{E}{\left[K(x,X_1,\{Z_j{\}}_{j=1}^s)\mid X_1\right]}^2\right]=\frac{1}{(2s-1)k^2}\left(\sum _{t=0}^{2k-2}\frac{a_t}{b_t}\right),$$

$$a_t=\sum _{i=\max \left\{0,t-(k-1)\right\}}^{\min \left\{t,k-1\right\}}\left(\genfrac{}{}{0ex}{}{s-1}{i}\right)\left(\genfrac{}{}{0ex}{}{s-1}{t-i}\right)\text{and}{b}_t=\sum _{i=0}^t\left(\genfrac{}{}{0ex}{}{s-1}{i}\right)\left(\genfrac{}{}{0ex}{}{s-1}{t-i}\right)$$

$${\sigma }_{n,j}^2(x)=\frac{s^2}{n}\frac{{\sigma }_j^2(x)}{k^2(2s-1)}{\zeta }_k+o(s/n),$$

$$\frac{{\hat{\theta }}_j(x)-{\theta }_j(x)}{{\sigma }_{n,j}(x)}\to \mathsf{N}(0,1).$$

$$\frac{{\hat{\theta }}_j-{\theta }_j}{{\tilde{\sigma }}_{n,j}^2(x)}=\frac{{\hat{\theta }}_j-{\theta }_j}{{\sigma }_{n,j}^2(x)}\frac{{\sigma }_{n,j}^2(x)}{{\tilde{\sigma }}_{n,j}^2(x)}{\to }_d\mathsf{N}(0,1).$$

$$\mu (B(x,r))=$$

$$\le$$

$$\mathcal{X}=\left\{(x_1,x_2,\dots ,x_D)\in {\mathbb{R}}^D:\sum _{i=1}^{D}1\left\{x_i\ne 0\right\}\le d\right\}.$$

$$\mathcal{X}={\cup }_{1\le i_1<i_2<\cdots i_d\le D}\left\{(x_1,x_2,\cdots ,x_D)\in {\mathbb{R}}^D:x_j=0,j\notin \left\{i_1,i_2,\dots ,i_d\right\}\right\}.$$

Full text

1 Introduction

Many non-parametric estimation problems in econometrics and causal inference can be formulated as finding a parameter vector $\theta(x)\in\mathbb{R}^{p}$ that is a solution to a set of conditional moment equations:

[TABLE]

when given $n$ i.i.d. samples $(Z_{1},\ldots,Z_{n})$ from the distribution of $Z$, where $\psi:\mathcal{Z}\times\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$ is a known vector valued moment function, $\mathcal{Z}$ is an arbitrary data space, $X\in\mathcal{X}\subset\mathbb{R}^{D}$ is the feature vector that is included $Z$. Examples include non-parametric regression111$Z=(X,Y)$, where $Y\in\mathbb{R}^{p}$ is the dependent variable, and $\psi(Z;\theta(x))=Y-\theta(x)$., quantile regression222$Z=(X,Y)$ and $\psi(Z;\theta(x))=1\{Y\leq\theta(x)\}-\alpha$, for some $\alpha\in[0,1]$ that denotes the target quantile., heterogeneous treatment effect estimation333$Z=(X,T,Y)$, where $T\in\mathbb{R}^{p}$ is a vector of treatments, and $\psi(Z;\theta(x))=(Y-\left\langle\theta(x),T\right\rangle)\,T$., instrumental variable regression444$Z=(X,T,W,Y)$, where $T\in\mathbb{R}$ is a treatment, $W\in\mathbb{R}$ an instrument and $\psi(Z;\theta(x))=(Y-\theta(x)\,T)\,W$., local maximum likelihood estimation555Where the distribution of $Z$ admits a known density $f(z;\theta(x))$ and $\psi(Z;\theta(x))=\nabla_{\theta}\log(f(Z;\theta(x))$. and estimation of structural econometric models (see e.g., Reiss and Wolak (2007) and examples in Chernozhukov et al. (2016, 2018b)). The study of such conditional moment restriction problems has a long history in econometrics (see e.g., Newey (1993); Ai and Chen (2003); Chen and Pouzo (2009); Chernozhukov et al. (2015)). However, the majority of the literature assumes that the conditioning variable $X$ is low dimensional, i.e. $D$ is a constant as the sample size $n$ grows (see e.g., Athey et al. (2019)). High dimensional variants have primarily been analyzed under parametric assumptions on $\theta(x)$, such as sparse linear forms (see e.g., Chernozhukov et al. (2018a)). There are some papers that address the fully non-parametric setup (see e.g., Lafferty and Wasserman (2008); Dasgupta and Freund (2008); Kpotufe (2011); Biau (2012); Scornet et al. (2015)) but those are focused on the estimation problem, and do not address inference (i.e., constructing asymptotically valid confidence intervals).

The goal of this work is to address estimation and inference in conditional moment models with a high-dimensional conditioning variable. As is obvious without any further structural assumptions on the problem, the exponential in dimension rates of approximately $n^{1/D}$ (see e.g., Stone (1982)) cannot be avoided. Thereby estimation is in-feasible even if $D$ grows very slowly with $n$. Our work, follows a long line of work in machine learning (Dasgupta and Freund, 2008; Kpotufe, 2011; Kpotufe and Garg, 2013), which is founded on the observation that in many practical applications, even though the variable $X$ is high-dimensional (e.g. an image), one typically expects that the coordinates of $X$ are highly correlated. The latter intuition is formally captured by assuming that the distribution of $X$ has a small doubling measure around the target point $x$.

We refer to the latter notion of dimension, as the intrinsic dimension of the problem. Such a notion has been studied in the statistical machine learning literature, so as to establish fast estimation rates in high-dimensional kernel regression settings (Dasgupta and Freund, 2008; Kpotufe, 2011; Kpotufe and Garg, 2013; Xue and Kpotufe, 2018; Chen and Shah, 2018; Kim et al., 2018; Jiang, 2017). However, these works solely address the problem of estimation and do not characterize the asymptotic distribution of the estimates, so as to enable inference, hypothesis testing and confidence interval construction. Moreover, they only address the regression setting and not the general conditional moment problem and consequently do not extend to quantile regression, instrumental variable regression or treatment effect estimation.

From the econometrics side, pioneering works of Wager and Athey (2018); Athey et al. (2019) address estimation and inference of conditional moment models with all the aforementioned desiderata that are required for the application of such methodologies to social sciences, albeit in the low dimensional regime. In particular, Wager and Athey (2018) consider regression and heterogeneous treatment effect estimation with a scalar $\theta(x)$ and prove $n^{1/D}$-asymptotic normality of a sub-sampled random forest based estimator and Athey et al. (2019) extend it to the general conditional moment settings.

These results have been extended and improved in multiple directions, such as improved estimation rates through local linear smoothing Friedberg et al. (2018), robustness to nuisance parameter estimation error Oprescu et al. (2018) and improved bias analysis via sub-sampled nearest neighbor estimation Fan et al. (2018). However, they all require low dimensional setting and the rate provided by the theoretical analysis is roughly $n^{-1/D}$, i.e. to get a confidence interval of length $\epsilon$ or an estimation error of $\epsilon$, one would need to collect $O(\epsilon^{-D})$ samples which is prohibitive in most target applications of machine learning based econometrics.

Hence, there is a strong need to provide theoretical results that justify the success of machine learning estimators for doing inference, via their adaptivity to some low dimensional hidden structure in the data. Our work makes a first step in this direction and provides estimation and asymptotic normality results for the general conditional moment problem, where the rate of estimation and the asymptotic variance depend only on the intrinsic dimension, independent of the explicit dimension of the conditioning variable.

Our analysis proceeds in four parts. First, we extend the results by Wager and Athey (2018); Athey et al. (2019) on the asymptotic normality of sub-sampled kernel estimators to the high-dimensional, low intrinsic dimension regime and to vector valued parameters $\theta(x)$. Concretely, when given a sample $S=(Z_{1},\ldots,Z_{n})$, our estimator is based on the approach proposed in Athey et al. (2019) of solving a locally weighted empirical version of the conditional moment restriction

[TABLE]

where $K(x,X_{i},S)$ captures proximity of $X_{i}$ to the target point $x$. The approach dates back to early work in statistics on local maximum likelihood estimation (Fan et al., 1998; Newey, 1994; Stone, 1977; Tibshirani and Hastie, 1987). As in Athey et al. (2019), we consider weights $K(x,X_{i},S)$ that take the form of an average over $B$ base weights: $K(x,X_{i},S)=\frac{1}{B}\sum_{b=1}^{B}K(x,X_{i},S_{b})\,1\{i\in S_{b}\}$, where each $K(x,X_{i},S_{b})$ is calculated based on a randomly drawn sub-sample $S_{b}$ of size $s<n$ from the original sample. We will typically refer to the function $K$ as the kernel. In Wager and Athey (2018); Athey et al. (2019) $K(x,X_{i},S_{b})$ is calculated by building a tree on the sub-sample, while in Fan et al. (2018) it is calculated based on the $1$-NN rule on the sub-sample.

Our main results are general estimation rate and asymptotic normality theorems for the estimator $\hat{\theta}(x)$ (see Theorems 1 and 2), which are stated in terms of two high-level assumptions, specifically an upper bound $\epsilon(s)$ on the rate at which the kernel “shrinks” and a lower bound $\eta(s)$ on the “incrementality” of the kernel. Notably, the explicit dimension of the conditioning variable $D$ does not enter the theorem, so it suffices in what follows to show that $\epsilon(s)$ and $\eta(s)$ depend only on $d$ rather than $D$.

The shrinkage rate $\epsilon(s)$ is defined as the $\ell_{2}$-distance between the target point $x$ and the furthest point on which the kernel places positive weight $X_{i}$, when trained on a data set of $s$ samples, i.e.,

[TABLE]

The shrinkage rate of the kernel controls the bias of the estimate (small $\epsilon(s)$ implies low bias). The sub-sampling size $s$ is a lever to trade off bias and variance; larger $s$ achieves smaller bias, since $\epsilon(s)$ is smaller, but increases the variance, since for any fixed $x$ the weights $K(x,X_{i},S_{b})$ will tend to concentrate on the same data points, rather than averaging over observations. Both estimation and asymptotic normality results require the bias to be controlled through the shrinkage rate.

Incrementality of a kernel describes how much information is revealed about the weight of a sample $i$ solely by knowledge of $X_{i}$, and is captured by the second moment of the conditional expected weight

[TABLE]

The incrementality assumption is used in the asymptotic normality proof to argue that the weights have sufficiently high variance that all data points have some influence on the estimate. From the technical side, we use the Hájek projection to analyze our $U$-statistic estimator. Incrementality ensures that there is sufficiently weak dependence in the weights across a sequence of sub-samples and hence the central limit theorem applies. As discussed, the sub-sampling size $s$ can be used to control the variance of the weights, and so incrementality and shrinkage are related. We make this precise, proving that incrementality can be lower bounded as a function of kernel shrinkage, so that having a sufficiently low shrinkage rate enables both estimation and inference. These general results could be of independent interest beyond the scope of this work.

For the second part of our analysis, we specialize to the case where the base kernel is the $k$-NN kernel, for some constant $k$. We prove that both shrinkage and incrementality depend only on the intrinsic dimension $d$, rather than the explicit dimension $D$. In particular, we show that $\epsilon(s)=O(s^{-1/d})$ and $\eta(s)=\Theta(1/s)$. These lead to our main theorem that the sub-sampled $k$-NN estimate achieves an estimation rate of order $n^{1/(d+2)}$ and is also $n^{1/(d+2)}$-asymptotically normal.

In the third part, we provide a closed form characterization of the asymptotic variance of the sub-sampled $k$-NN estimate, as a function of the conditional variance of the moments, which is defined as $\sigma^{2}(x)=\operatorname{Var}\left(\psi(Z;\theta)\mid X=x\right)$. For example, for the $1$-NN kernel, the asymptotic variance is given by

[TABLE]

This strengthens prior results of Fan et al. (2018) and Wager and Athey (2018), which only proved the existence of an asymptotic variance without providing an explicit form (and thereby relied on bootstrap approaches for the construction of confidence intervals). Our tight characterization enables an easy construction of plugin normal-based intervals that only require a preliminary estimate of $\sigma(x)$. Our Monte Carlo study shows that such intervals provide very good finite sample coverage in a high dimensional regression setup (see Figure 1)666See Appendix C for detailed explanation of our simulations..

Finally in the last part, we discuss an adaptive data-driven approach for picking the sub-sample size $s$ so as to achieve estimation or asymptotic normality with rates that only depend on the unknown intrinsic dimension. This allows us to achieve near-optimal rates while adapting to the unknown intrinsic dimension of data (see Propositions 1 and 2). Figure 2 depicts the performance of our adaptive approach compared to two benchmarks, one constructed based on theory for intrinsic dimension $d$ which may be unknown, and the other one constructed naïvely based on the known but sub-optimal extrinsic dimension $D$. As it can be observed from this figure, setting $s$ based on intrinsic dimension $d$ allows us to build more accurate and smaller confidence intervals, which is crucial for drawing inference in the high-dimensional finite sample regime. Our adaptive approach uses samples to pick $s$ very close to the value suggested by our theory and therefore leads to a compelling finite sample coverage777A preliminary implementation of our code is available via http://github.com/khashayarkhv/np_inference_intrinsic..

Our results shed some light on the importance of using adaptive machine learning based estimators, such as nearest neighbor based estimates, when performing estimation and inference in high-dimensional econometric problems. Such estimators address the curse of dimensionality by adapting to a priori unknown latent structure in the data. Moreover, coupled with the powerful sub-sampling based averaging approach, such estimators can maintain their adaptivity, while also satisfying asymptotic normality and thereby enabling asymptotically valid inference; a property that is crucial for embracing such approaches in econometrics and causal inference.

Structure of the paper.

The rest of the paper is organized as follows. In §2, we provide preliminary definitions, in §2.1 and §2.2 we explain our algorithms, in §2.3 we explain doubling dimension (see Appendix B for examples). In §3 we state our assumptions, in §4 we provide general estimation and inference results for kernels that satisfy shrinkage and incrementality conditions, and in §5 we apply such results to the $k$-NN kernel and prove estimation and inference rates for such kernels that only depend on intrinsic dimension. We defer a discussion on the extension to heterogeneous treatment effect estimation and also technical proofs to Appendices.

2 Preliminaries

Suppose we have a data set $M$ of $n$ observations $Z_{1},Z_{2},\ldots,Z_{n}$ drawn independently from some distribution $\mathcal{D}$ over the observation domain $\mathcal{Z}$. We focus on the case that $Z_{i}=(X_{i},Y_{i})$, where $X_{i}$ is the vector of covariates and $Y_{i}$ is the outcome. In Appendix D, we briefly discuss how our results can be extended to the setting where nuisance parameters and treatments are included in the model.

Suppose that the covariate space $\mathcal{X}\subset\mathbb{R}^{D}$ is contained in a ball with unknown diameter $\Delta_{\mathcal{X}}$. Denote the marginal distribution of $X$ by $\mu$ and the empirical distribution of $X$ on $n$ sample points by $\mu_{n}$. Let $B(x,r)=\left\{z\in\mathbb{R}^{D}:\|x-z\|_{2}<r\right\}$ be the $\ell_{2}$-ball centered at $x$ with radius $r$ and denote the standard basis for $\mathbb{R}^{p}$ by $\left\{e_{1},e_{2},\ldots,e_{p}\right\}$.

Let $\psi:\mathcal{Z}\times\mathbb{R}^{p}\to\mathbb{R}^{p}$ be a score function that maps observation $Z$ and parameter $\theta\in\mathbb{R}^{p}$ to a $p$-dimensional score $\psi(Z;\theta)$. For $x\in\mathcal{X}$ and $\theta\in\mathbb{R}^{p}$ define the expected score as $m(x;\theta)=\mathbf{E}[\psi(Z;\theta)\mid X=x]$. The goal is to estimate the quantity $\theta(x)$ via local moment condition, i.e.

[TABLE]

2.1 Sub-Sampled Kernel Estimation

Base Kernel Learner.

Our learner $\mathcal{L}_{k}$ takes a data set $S$ containing $m$ observations as input and a realization of internal randomness $\omega$, and outputs a kernel weighting function $K_{\omega}:\mathcal{X}\times\mathcal{X}\times\mathcal{Z}^{m}\to[0,1]$. In particular, given any target feature $x$ and the set $S$, the weight of each observation $Z_{i}$ in $S$ with feature vector $X_{i}$ is $K_{\omega}(x,X_{i},S)$. Define the weighted score on a set $S$ with internal randomness $\omega$ as $\Psi_{S}(x;\theta)=\sum_{i\in S}K_{\omega}(x,X_{i},S)\psi(Z_{i};\theta)$. When it is clear from context we will omit $\omega$ from our notation for succinctness and essentially treat $K$ as a random function. For the rest of the paper, we are going to use notations $\alpha_{S,\omega}(X_{i})=K_{\omega}(x,X_{i},S)$ interchangeably.

Averaging over $B$ sub-samples of size $s$.

Suppose that we consider $B$ random and independent draws from all $\binom{n}{s}$ possible subsets of size $s$ and internal randomness variables $\omega$ and look at their average. Index these draws by $b=1,2,\ldots,B$ where $S_{b}$ contains samples in $b$th draw and $\omega_{b}$ is the corresponding draw of internal randomness. We can define the weighted score as

[TABLE]

Estimating $\theta(x)$.

We estimate $\theta(x)$ as a vanishing point of $\Psi(x;\theta)$. Letting $\hat{\theta}$ be this point, then $\Psi(x;\hat{\theta})=\frac{1}{B}\sum_{b=1}^{B}\sum_{i=1}^{n}\alpha_{S_{b},\omega_{b}}(X_{i})\psi(Z_{i};\hat{\theta})=0.$ This procedure is explained in Algorithm 1.

2.2 Sub-Sampled $k$-NN Estimation

We specially focus on the case that the weights are distributed across the $k$-NN of $x$. In other words, given a data set $S$, the weights are given according to $K_{\omega}(x,X_{i},S)=\,1\left\{X_{i}\in H_{k}(x,S)\right\}/k$, where $H_{k}(x,S)$ are $k$-NN of $x$ in the set $S$. The pseudo-code for this can be found in Algorithm 2.

Complete $U$-statistic.

The expression in Equation (5) is an incomplete $U$-statistic. Complete $U$-statistic is obtained if we allow each subset of size $s$ from $n$ samples to be included in the model exactly once. In other words, this is achieved if $B=\binom{n}{s}$, all subsets $S_{1},S_{2},\ldots,S_{B}$ are distinct, and we also take expectation over the internal randomness $\omega$. Denoting this by $\Psi_{0}(x;\theta)$, we have

[TABLE]

Note in the case of $k$-NN estimator we can also represent $\Psi_{0}$ in terms of order statistics, i.e., $\Psi_{0}$ is an $L$-statistics (see e.g., Serfling (2009)). By sorting samples in $\mathbf{X}=\left\{X_{1},X_{2},\ldots,X_{n}\right\}$ based on their distance with $x$ as $\|X_{(1)}-x\|\leq\|X_{(2)}-x\|\leq\cdots\leq\|X_{(n)}-x\|$, we can write $\Psi_{0}(x;\theta)=\sum_{i=1}^{n}\alpha(X_{(i)})\,\psi(Z_{(i)};\theta)$ where the weights are given by

[TABLE]

2.3 Local intrinsic dimension

We are interested in settings that the distribution of $X$ has some low dimensional structure on a ball around the target point $x$. The following notions are adapted from Kpotufe (2011), which we present here for completeness.

Definition 1**.**

The marginal $\mu$ is called doubling measure if there exists a constant $C_{db}>0$ such that for any $x\in\mathcal{X}$ and any $r>0$ we have $\mu(B(x,r))\leq C_{db}\mu(B(x,r/2))$.

An equivalent definition of this notion is that, the measure $\mu$ is doubling measure if there exist $C,d>0$ such that for any $x\in\mathcal{X},r>0$, and $\theta\in(0,1)$ we have $\mu(B(x,r))\leq C\theta^{-d}\mu(B(x,\theta r))$.

One example is given by Lebesgue measure on the Euclidean space $\mathbb{R}^{d}$, where for any $r>0,\theta\in(0,1)$ we have $\operatorname{vol}(B(x,\theta r))=\operatorname{vol}(B(x,r))\theta^{d}$. Building upon this, let $\mathcal{X}\in\mathbb{R}^{D}$ be a subset of $d$-dimensional hyperplane and suppose that for any ball $B(x,r)$ in $\mathcal{X}$ we have $\operatorname{vol}(B(x,r)\cap\mathcal{X})=\Theta(r^{d})$. If $\mu$ is almost uniform, then we also have $\mu(B(x,\theta r))/\mu(B(x,r))=\Theta(\theta^{d})$.

Unfortunately, this global notion of doubling measure is restrictive and most probability measures are globally complex. Rather, once restricted to local neighborhoods, the measure becomes lower dimensional and intrinsically less complex. The following definition captures this intuition better.

Definition 2**.**

Fix $x\in\mathcal{X}$ and $r>0$. The marginal $\mu$ is $(C,d)$-homogeneous on $B(x,r)$ if for any $\theta\in(0,1)$ we have $\mu(B(x,r))\leq C\theta^{-d}\mu(B(x,\theta r))$.

Intuitively, this definition requires the marginal $\mu$ to have a local support that is intrinsically $d$-dimensional. This definition covers low-dimensional manifolds, mixture distributions, $d$-sparse data, and also any combination of these examples. These examples are explained in Appendix B.

3 Assumptions

For non-parametric estimators the bias is connected to the kernel shrinkage, as noted by Athey et al. (2019); Wager and Athey (2018); Oprescu et al. (2018).

Definition 3** (Kernel Shrinkage in Expectation).**

The kernel weighting function output by learner $\mathcal{L}_{k}$ when it is given $s$ i.i.d. observations drawn from distribution $\mathcal{D}$ satisfies

[TABLE]

Definition 4** (Kernel Shrinkage with High Probability).**

The kernel weighting function output by learner $\mathcal{L}_{k}$ when it is given $s$ i.i.d. observations drawn from distribution $\mathcal{D}$ w.p. $1-\delta$ over the draws of the $s$ samples satisfies

[TABLE]

As shown in Wager and Athey (2018), for trees that satisfy some regularity condition, $\epsilon(s)\leq s^{-c/D}$ for a constant $c$. We are interested in shrinkage rates that scale as $s^{-c/d}$, where $d$ is the local intrinsic dimension of $\mu$ on $B(x,r)$. Similar to Oprescu et al. (2018); Athey et al. (2019), we rely on the following assumptions on the moment and score functions.

Assumption 1**.**

The moment $m(x;\theta)$ corresponds to the gradient w.r.t. $\theta$ of a $\lambda$-strongly convex loss $L(x;\theta)$. This also means that the Jacobian $M_{0}=\nabla_{\theta}m(x;\theta(x))$ has minimum eigenvalue at least $\lambda$. 3. 2.

For any fixed parameters $\theta$, $m(x;\theta)$ is a $L_{m}$-Lipschitz function in $x$ for some constant $L_{m}$. 4. 3.

There exists a bound $\psi_{\max}$ such that for any observation $z$ and any $\theta$, $\|\psi(z;\theta)\|_{\infty}\leq\psi_{\max}$. 5. 4.

The bracketing number $N_{[]}({\cal F},\epsilon,L_{2})$ of the function class: ${\cal F}=\{\psi(\cdot;\theta):\theta\in\Theta\}$, satisfies $\log(N_{[]}({\cal F},\epsilon,L_{2}))=O(1/\epsilon)$.

Assumption 2**.**

For any coordinate $j$ of the moment vector $m$, the Hessian $H_{j}(x;\theta)=\nabla^{2}_{\theta\theta}m_{j}(x;\theta)$ has eigenvalues bounded above by a constant $L_{H}$ for all $\theta$. 3. 2.

Maximum eigenvalue of $M_{0}$ is upper bounded by $L_{J}$. 4. 3.

Second moment of $\psi(x;\theta)$ defined as $\operatorname{Var}\left(\psi(Z;\theta)\mid X=x\right)$ is $L_{mm}$-Lipschitz in $x$, i.e.,

[TABLE] 5. 4.

Variogram is Lipschitz: $\sup_{x\in\mathcal{X}}\|\operatorname{Var}(\psi(Z;\theta)-\psi(Z;\theta^{\prime})\mid X=x)\|_{F}\leq L_{\psi}\|\theta-\theta^{\prime}\|_{2}$.

The condition on variogram always holds for a $\psi$ that is Lipschitz in $\theta$. This larger class of functions $\psi$ allows estimation in more general settings such as $\alpha$-quantile regression that involves a $\psi$ which is non-Lipschitz in $\theta$. Similar to Athey and Imbens (2016); Athey et al. (2019), we require kernel $K$ to be honest and symmetric.

Assumption 3**.**

The kernel $K$, built using samples $\left\{Z_{1},Z_{2},\ldots,Z_{s}\right\}$, is honest if the weight of sample $i$ given by $K(x,X_{i},\left\{Z_{j}\right\}_{j=1}^{s})$ is independent of $Y_{j}$ conditional on $X_{j}$ for any $j\in[s]$.

Assumption 4**.**

The kernel $K$, built using samples $\left\{Z_{1},Z_{2},\ldots,Z_{s}\right\}$, is symmetric if for any permutation $\pi:[s]\to[s]$, the distribution of $K(x,X_{i},\left\{Z_{j}\right\}_{j=1}^{s})$ and $K(x,X_{\pi(i)},\left\{Z_{\pi(j)}\right\}_{j=1}^{s})$ are equal. In other words, the kernel weighting distribution remains unchanged under permutations.

For a deterministic kernel $K$, the above condition implies that $K(x,X_{i},\{Z_{j}\}_{j=1}^{s})=K(x,X_{i},\{Z_{\pi(j)}\}_{j=1}^{s})$, for any $i\in[s]$. In the next section, we provide general estimation and inference results for a general kernel based on the its shrinkage and incrementality rates.

4 Guarantees for sub-sampled kernel estimators

Our first result establishes estimation rates, both in expectation and high probability, for kernels based on their shrinkage rates. The proof of this theorem is deferred to Appendix E.

Theorem 1** (Finite Sample Estimation Rate).**

Let Assumptions 1 and 3 hold. Suppose that Algorithm 1 is executed with $B\geq n/s$. If the base kernel $K$ satisfies kernel shrinkage in expectation, with rate $\epsilon(s)$, then w.p. $1-\delta$

[TABLE]

Moreover,

[TABLE]

The next result establishes asymptotic normality of sub-sampled kernel estimators. In particular, it provides coordinate-wise asymptotic normality of our estimate $\hat{\theta}$ around its true underlying value $\theta(x)$. The proof of this theorem is deferred to Appendix F.

Theorem 2** (Asymptotic Normality).**

Let Assumptions 1, 2, 3, and 4 hold. Suppose that Algorithm 1 is executed with $B\geq(n/s)^{5/4}$ and the base kernel $K$ satisfies kernel shrinkage, with rate $\epsilon(s,\delta)$ in probability and $\epsilon(s)$ in expectation. Let $\eta(s)$ be the incrementality of kernel $K$ defined in Equation (4) and $s$ grow at a rate such that $s\rightarrow\infty$, $n\eta(s)\rightarrow\infty$, and $\epsilon(s,\eta(s)^{2})\rightarrow 0$. Consider any fixed coefficient $\beta\in\mathbb{R}^{p}$ with $\|\beta\|\leq 1$ and define the variance as

[TABLE]

Then it holds that $\sigma_{n,\beta}(x)=\Omega\left(s\sqrt{\eta(s)/n}\right)$. Moreover, suppose that

[TABLE]

Then,

[TABLE]

Theorems 1 and 2 generalize existing estimation and asymptotic normality results of Athey et al. (2019); Wager and Athey (2018); Fan et al. (2018) to an arbitrary kernel that satisfies appropriate shrinkage and incrementality rates (see Remark 1 in Appendix F). The following lemma relates these two and provides a lower bound on the incrementality in terms of kernel shrinkage. The proof uses the Paley-Zygmund inequality and is left to Appendix G.

Lemma 1**.**

For any symmetric kernel $K$ (Assumption 4) and for any $\delta\in[0,1]$:

[TABLE]

Thus if $\mu(B(x,\epsilon(s,1/(2s^{2}))))=O(\log(s)/s)$, then picking $\rho=1/(2s^{2})$ and $\delta=1/2$ implies that $\mathbf{E}[\mathbf{E}[K(x,X_{1},\{Z_{j}\}_{j=1}^{s})|X_{1}]^{2}]=\Omega(1/s\log(s)).$

Corollary 1**.**

If $\epsilon(s,\delta)=O((\log(1/\delta)/s)^{1/d})$ and $\mu$ satisfies a two-sided version of the doubling measure property on $B(x,r)$, defined in Definition 2, i.e., $c\theta^{d}\mu(B(x,r))\leq\mu(B(x,\theta r))\leq C\theta^{d}\mu(B(x,r))$ for any $\theta\in(0,1)$. Then, $\mathbf{E}[\mathbf{E}[K(x,X_{1},\{Z_{j}\}_{j=1}^{s})|X_{1}]^{2}]=\Omega(1/(s\log(s)))$.

Even without this extra assumption, we can still characterize the incrementality rate of the $k$-NN estimator, as we observe in the next section.

5 Main theorem: adaptivity of $k$-NN estimator

In this section, we provide estimation guarantees and asymptotic normality of the $k$-NN estimator by using Theorems 1 and 2. We first establish shrinkage and incrementality rates for this kernel.

5.1 Estimation guarantees for the $k$-NN estimator

In this section we provide shrinkage results for the $k$-NN kernel. As observed in Theorem 1, shrinkage rates are sufficient for bounding the estimation error. The shrinkage result that we present in the following would only depend on the local intrinsic dimension of $\mu$ on $B(x,r)$.

Lemma 2** (High probability shrinkage for the $k$-NN kernel).**

Suppose that the measure $\mu$ is $(C,d)$-homogeneous on $B(x,r)$. Then, for any $\delta$ satisfying $2\exp\left(-\mu(B(x,r))s/(8C)\right)\leq\delta\leq\frac{1}{2}\exp(-k/2),$ w.p. at least $1-\delta$ we have

[TABLE]

We can turn this into a shrinkage rate in expectation as follows. In fact, by the very convenient choice of $\delta=s^{-1/d}$ combined with the fact that $\mathcal{X}$ has diameter $\Delta_{\mathcal{X}}$, we can establish $O\left((\log(s)/s)^{1/d}\right)$ rate on expected kernel shrinkage. However, a more careful analysis would help us to remove the $\log(s)$ dependency in the bound and is stated in the following corollary:

Corollary 2** (Expected shrinkage for the $k$-NN kernel).**

Suppose that the conditions of Lemma 2 hold. Let $k$ be a constant and $\epsilon_{k}(s)$ be the expected shrinkage for the $k$-NN kernel. Then, for any $s$ larger than some constant we have $\epsilon_{k}(s)=\mathbf{E}\left[\|x-X_{(k)}\|_{2}\right]=O\left(\frac{1}{s}\right)^{1/d}$.

We are now ready to state our estimation result for the $k$-NN kernel, which is honest and symmetric. Therefore, we can substitute the expected shrinkage rate established in Corollary 2 in Theorem 1 to derive estimation rates for this kernel.

Theorem 3** (Estimation Guarantees for the $k$-NN Kernel).**

Suppose that $\mu$ is $(C,d)$-homogeneous on $B(x,r)$, Assumption 1 holds and that Algorithm 2 is executed with $B\geq n/s$. Then, w.p. $1-\delta$:

[TABLE]

and

[TABLE]

*By picking $s=\Theta\left(n^{d/(d+2)}\right)$ and $B=\Omega\left(n^{2/(d+2)}\right)$ we get $\sqrt{\mathbf{E}\left[\|\hat{\theta}-\theta(x)\|_{2}^{2}\right]}=\tilde{O}\left(n^{-1/(d+2)}\right)$. *

5.2 Asymptotic normality of the $k$-NN estimator

In this section we prove asymptotic normality of $k$-NN estimator. We start by provide bounds on the incrementality of the $k$-NN kernel.

Lemma 3** ($k$-NN Incrementality).**

Let $K$ be the $k$-NN kernel and let $\eta_{k}(s)$ denote the incrementality rate of this kernel. Then, the following holds:

[TABLE]

where sequences $\left\{a_{t}\right\}_{t=0}^{2k-2}$ and $\left\{b_{t}\right\}_{t=0}^{2k-2}$ are defined as

[TABLE]

We can substitute $\eta_{k}(s)$ in Theorem 2 to prove asymptotic normality of the $k$-NN estimator. Before that, we establish the asymptotic variance of this estimator $\sigma_{n,j}(x)$, up to the smaller order terms.

Theorem 4** (Asymptotic Variance of $k$-NN).**

Let $j\in[p]$ be one of coordinates. Suppose that $k$ is constant while $s\rightarrow\infty$. Then, for the $k$-NN kernel

[TABLE]

where $\sigma_{j}^{2}(x)=\operatorname{Var}\left[\left\langle e_{j},M_{0}^{-1}\psi(Z;\theta(x))\right\rangle\mid X=x\right]$ and $\zeta_{k}=k+\sum_{t=k}^{2k-2}2^{-t}\sum_{i=t-k+1}^{k-1}\binom{t}{i}.$

Combining results of Theorem 2, Theorem 4, Corollary 2, and Lemma 3 we have:

Theorem 5** (Asymptotic Normality of $k$-NN Estimator).**

Suppose that $\mu$ is $(C,d)$-homogeneous on $B(x,r)$. Let Assumptions 1, 2 hold and suppose that Algorithm 2 is executed with $B\geq(n/s)^{5/4}$ iterations. Suppose that $s$ grows at a rate such that $s\rightarrow\infty$, $n/s\rightarrow\infty$, and also $s^{-1/d}(n/s)^{1/2}\rightarrow 0$. Let $j\in[p]$ be one of coordinates and $\sigma^{2}_{n,j}(x)$ be defined in Equation (12). Then,

[TABLE]

Finally, if $s=n^{\beta}$ and $B\geq n^{\frac{5}{4}(1-\beta)}$ with $\beta\in(d/(d+2),1)$. Then,

[TABLE]

Plug-in confidence intervals.

Observe that the Theorem 4 implies that if we define $\tilde{\sigma}_{n,j}^{2}(x)=\frac{s^{2}}{n}\frac{\sigma_{j}^{2}(x)}{2s-1}\,\frac{\zeta_{k}}{k^{2}}$ as the leading term in the variance, then $\frac{\sigma_{n,j}^{2}(x)}{\tilde{\sigma}_{n,j}^{2}(x)}\rightarrow_{p}1$. Thus, due to Slutsky’s theorem

[TABLE]

Hence, we have a closed form solution to the variance in our asymptotic normality theorem. If we have an estimate $\hat{\sigma}_{j}^{2}(x)$ of the variance of the conditional moment around $x$, then we can build plug-in confidence intervals based on the normal distribution with variance $\frac{s^{2}}{n}\frac{\hat{\sigma}_{j}^{2}(x)}{2s-1}\,\frac{\zeta_{k}}{k^{2}}$. Note that $\zeta_{k}$ can be calculated easily for desired values of $k$. For instance, we have $\zeta_{1}=1,\zeta_{2}=\frac{5}{2},$ and $\zeta_{3}=\frac{33}{8}$ and for $k=1,2,3$ the asymptotic variance becomes $\frac{s^{2}}{n}\frac{\hat{\sigma}_{j}^{2}(x)}{2s-1},\frac{5}{8}\frac{s^{2}}{n}\frac{\hat{\sigma}_{j}^{2}(x)}{2s-1},$ and $\frac{11}{24}\frac{s^{2}}{n}\frac{\hat{\sigma}_{j}^{2}(x)}{2s-1}$ respectively.

5.3 Adaptive choice for $s$

According to Theorem 3, picking $s=\Theta(n^{d/(d+2)})$ would trade-off between bias and variance terms. Also, according to Theorem 5, picking $s=n^{\beta}$ with $d/(d+2)<\beta<1$ would result in asymptotic normality of the estimator. However, both choices depend on the unknown intrinsic dimension of $\mu$ on the ball $B(x,r)$. Inspired by Kpotufe (2011), we explain a data-driven way for estimating $s$.

Suppose that $\delta>0$ is given. Let $C_{n,p,\delta}=2\log(2pn/\delta)$ and pick $\Delta\geq\Delta_{\mathcal{X}}$. For any $k\leq s\leq n$, let $H(s)$ be the $U$-statistic estimator for $\epsilon(s)$ defined as $H(s)=\sum_{S\in[n]:|S|=s}\max_{X_{i}\in H_{k}(x,S)}\|x-X_{i}\|_{2}/\binom{n}{s}$. Each term in the summation computes the distance of $x$ to its $k$-nearest neighbor on $S$ and $H(s)$ is the average of these numbers over all $\binom{n}{s}$ possible subsets $S$ (see Remark 3 in Appendix H regarding to efficient computation of $H(s)$). Define $G_{\delta}(s)=\Delta\sqrt{C_{n,p,\delta}ps/n}$. Iterate over $s=n,\cdots,k$. Let $s_{2}$ be the smallest $s$ for which we have $H(s)>2G_{\delta}(s)$ and let $s_{1}=s_{2}+1$. Note that $\epsilon_{k}(s)$ is decreasing in $s$ and $G_{\delta}(s)$ is increasing in $s$. Therefore, there exists a unique $1\leq s^{*}\leq n$ such that $\epsilon_{k}(s^{*})\leq G_{\delta}(s^{*})$ and $\epsilon_{k}(s^{*}-1)>G_{\delta}(s^{*}-1)$. We have following results.

Proposition 1** (Adaptive Estimation).**

Let Assumptions of Theorem 3 hold. Suppose that $s_{1}$ is the output of the above process. Let $s_{*}=9s_{1}+1$ and suppose that Algorithm 2 is executed with $s=s_{*}$ and $B\geq n/s_{*}$. Then w.p. at least $1-2\delta$ we have $\|\hat{\theta}-\theta(x)\|_{2}=O(G_{\delta}(s^{*}))=O\left(\left(\frac{n}{p\log(2pn/\delta)}\right)^{-1/(d+2)}\right)$. Further, for $\delta=1/n$ we have $\sqrt{\mathbf{E}\left[\|\hat{\theta}-\theta(x)\|_{2}^{2}\right]}=\tilde{O}\left(n^{-1/(d+2)}\right)$.

Proposition 2** (Adaptive Asymptotic Normality).**

Let Assumptions of Theorem 5 hold. Suppose that $s_{1}$ is the output of the above process when $\delta=1/n$ and $s_{*}=9s_{1}+1$. For any $\zeta\in(0,(\log(n)-\log(s_{1})-\log\log^{2}(n))/\log(n)))$ define $s_{\zeta}=s_{*}n^{\zeta}$. Suppose that Algorithm 2 is executed with $s=s_{\zeta}$ and $B\geq(n/s_{\zeta})^{5/4}$, then for any coordinate $j\in[p]$, we have $\frac{\hat{\theta}_{j}(x)-\theta_{j}(x)}{\sigma_{n,j}(x)}\rightarrow\mathsf{N}(0,1)\,.$

Appendix A Related work

There exists a vast literature on average treatment effect estimation in high-dimensional settings. The key challenge in such settings is the problem of overfitting which is usually handled by adding regularization terms. However, this leads to a shrinked estimate for the average treatment effect and therefore not desirable. The literature has taken various approaches to solve this issue. For instance, Belloni et al. [2014a, b] used a two-step method for estimating average treatment effect where in the first step feature-selection is accomplished via a lasso and then treatment effect is estimated using selected features. Athey et al. [2018] studied approximate residual balancing where a combination of weight balancing and regression adjustment is used for removing undesired bias and for achieving a double robust estimator. Chernozhukov et al. [2018a, a] considered a more general semi-parametric framework and studied debiased/double machine learning methods via first order Neyman orthogonality condition. Mackey et al. [2017] extended this result to higher order moments. Please refer to Athey and Imbens [2017], Mullainathan and Spiess [2017], Belloni et al. [2017] for a review on this literature.

However, in many applications, researchers are interested in estimating heterogeneous treatment effect on various sub-populations. One effective solution is to use one of the methods described in previous paragraph to estimate problem parameters and then project such estimations onto the sub-population of interest. However, these approaches usually perform poorly when there is a model mis-specification, i.e., when the true underlying model does not belong to the parametric search space. Consequently, researchers have studied non-parametric estimators such as $k$-NN estimators, kernel estimators, and random forests. While these non-parametric estimators are very robust to model mis-specification and work well under mild assumptions on the function of interest, they suffer from the curse of dimensionality (see e.g., Bellman [1961], Robins and Ritov [1997], Friedman et al. [2001]). Therefore, for applying these estimators in high-dimensional settings it is necessary to design and study non-parametric estimators that are able to overcome curse of dimensionality when possible.

The seminal work of Wager and Athey [2018] utilized random forests originally introduced by Breiman [2001] and adapted them nicely for estimating heterogeneous treatment effect. In particular, the authors demonstrated how the recursive partitioning idea, explained in Athey and Imbens [2016] for estimating heterogeneity in causal settings, can be further analyzed to establish asymptotic properties of such estimators. The main premise of random forests is that they are able to adaptively select nearest neighbors and that is very desirable in high-dimensional settings where discarding uninformative features is necessary for combating the curse of dimensionality. In a follow-up work, they extended these results and introduced Generalized Random Forests for more general setting of solving generalized method of moment (GMM) equations Athey et al. [2019]. There has been some interesting developments of such ideas to other settings. Fan et al. [2018] introduced Distributional Nearest Neighbor (DNN) where they used $1$-NN estimators together with sub-sampling and explained that by precisely combining two of these estimators for different sub-sampling sizes, the first order bias term can be efficiently removed. Friedberg et al. [2018] paired this idea with a local linear regression adjustment and introduced Local Linear Forests in order to improve forest estimations for smooth functions. Oprescu et al. [2018] incorporated the double machine learning methods of Chernozhukov et al. [2018a] into GMM framework of Athey et al. [2019] and studied Orthogonal Random Forests in partially linear regression models with high-dimensional controls. Although forest kernels studied in Wager and Athey [2018] and Athey et al. [2019] seem to work well in high-dimensional applications, to the best of our knowledge, there still does not exists a theoretical result supporting it. In fact, all existing theoretical results suffer from the curse of dimensionality as they depend on the dimension of problem $D$.

The literature on machine learning and non-parametric statistics has recently studied how these worst-case performances can be avoided when the intrinsic dimension of problem is smaller than $D$. Please refer to Cutler [1993] for different notions of intrinsic dimension in metric spaces. Dasgupta and Freund [2008] studied random projection trees and showed that the structure of these trees do not depend on the actual dimension $D$, but rather on the intrinsic dimension $d$. They used the notion of Assouad Dimension, introduced by Assouad [1983], and proved that using random directions for splitting, the number of levels required for halving the diameter of a leaf scales as $O(d\log d)$. The follow-up work Verma et al. [2009] generalized these results for some other notions of dimension. Kpotufe and Dasgupta [2012] extended this idea to the regression setting and proved integrated risk bounds for random projection trees that were only dependent on intrinsic dimension. Kpotufe [2011], Kpotufe and Garg [2013] studied this in the context of $k$-NN and kernel estimations and established uniform point-wise risk bounds only depending on the local intrinsic dimension.

Our work is deeply rooted in the literature on intrinsic dimension explained above, literature on $k$-NN estimators (see e.g, Mack [1981], Samworth et al. [2012], Györfi et al. [2006], Biau and Devroye [2015], Berrett et al. [2019], Fan et al. [2018]), and generalized method of moments (see e.g., [Tibshirani and Hastie, 1987, Staniswalis, 1989, Fan et al., 1998, Hansen, 1982, Stone, 1977, Lewbel, 2007, Mackey et al., 2017]). We adapt the framework of Athey et al. [2019] and Oprescu et al. [2018] and solve a generalized moment problem using a DNN estimator, originally introduced and studied by Fan et al. [2018]. We establish consistency and inference properties of this estimator and prove that these properties only depend on the local intrinsic dimension of problem. In particular, we prove that the finite sample estimation error of order $n^{-1/(d+2)}$ together with $n^{1/(d+2)}$-asymptotically normality result of DNN estimator for solving the generalized moment problem regardless of how big the actual dimension $D$ is.

Our result differs from existing literature on intrinsic dimension (e.g., Kpotufe [2011], Kpotufe and Garg [2013]) since in addition to estimation guarantees for the regression setting, we also allow valid inference in solving conditional moment equations. Our asymptotic normality result is different from existing results for $k$-NN (see e.g., Mack [1981]), generalized method of moments (see e.g., Lewbel [2007]). This paper complements the work of Fan et al. [2018] and extends it to the generalized method of moment setting. Furthermore, we relax the common assumption on the existence of density for covariates and prove that DNN estimators are adaptive to intrinsic dimension.

We also provide the exact expression for the asymptotic variance of DNN estimator built using a $k$-NN kernel, which enables plug-in construction of confidence intervals, rather than the bootstrap method of [Efron, 1982] which was used by [Wager and Athey, 2018, Athey et al., 2019, Fan et al., 2018]. While establishing consistency and asymptotic normality of our estimator, we also provide more general bounds on kernel shrinkage rate and also incrementality which can be useful for establishing asymptotic properties in other applications. One such application is given in high-dimensional settings where the exact nearest neighbor search is computationally expensive and Approximate Nearest Neighbor (ANN) search is often replaced in order to reduce this cost. Our flexible result allows us to use the state-of-the-art ANN algorithms (see e.g., Andoni et al. [2017, 2018]) while maintaining consistency and asymptotic normality.

Appendix B Examples of spaces with small intrinsic dimension

In this section we provide examples of metric spaces that have small local intrinsic dimension. Our first example covers the setting where the distribution of data lies on a low-dimensional manifold (see e.g., Roweis and Saul [2000], Tenenbaum et al. [2000], Belkin and Niyogi [2003]). For instance, this happens for image inputs. Even though images are often high-dimensional (e.g., $4096$ in the case of $64$ by $64$ images), all these images belong intrinsically to a $3$-dimensional manifold.

Example 1** (Low dimensional manifold (adapted from Kpotufe [2011])).**

Consider a $d$-dimensional submanifold $\mathcal{X}\subset\mathbb{R}^{D}$ and let $\mu$ have lower and upper bounded density on $\mathcal{X}$. The local intrinsic dimension of $\mu$ on $B(x,r)$ is $d$, provided that $r$ is chosen small enough and some conditions on curvature hold. In fact, Bishop-Gromov theorem (see e.g., Carmo [1992]) implies that under such conditions, the volume of ball $B(x,r)\cap\mathcal{X}$ is $\Theta(r^{d})$. This together with the lower and upper bound on the density implies that $\mu(B(x,r)\cap\mathcal{X})/\mu(B(x,\theta r)\cap\mathcal{X})=\Theta(\theta^{d})$, i.e. $\mu$ is $(C,d)$-homogeneous on $B(x,r)$ for some $C>0$.

Another example which happens in many applications, is sparse data. For example, in the bag of words representation of text documents, we usually have a vocabulary consisting of $D$ words. Although $D$ is usually large, each text document contains only a small number of these words. In this application, we expect our data (and measure) to have smaller intrinsic dimension. Before stating this example, let us discuss a more general example about mixture distributions.

Example 2** (Mixture distributions (adapted from Kpotufe [2011])).**

Consider any mixture distribution $\mu=\sum_{i}\pi_{i}\mu_{i}$, with each $\mu_{i}$ defined on $\mathcal{X}$ with potentially different supports. Consider a point $x$ and note that if $x\not\in\operatorname{supp}(\mu_{i})$, then there exists a ball $B(x,r_{i})$ such that $\mu_{i}(B(x,r_{i}))=0$. This is true since the support of any probability measure is always closed, meaning that its complement is an open set. Now suppose that $r$ is chosen small enough such that for any $i$ satisfying $x\in\operatorname{supp}(\mu_{i})$, $\mu_{i}$ is $(C_{i},d_{i})$-homogeneous on $B(x,r)$, while for any $i$ satisfying $x\not\in\operatorname{supp}(\mu_{i})$ we have $\mu_{i}(B(x,r))=0$. Then,

[TABLE]

where $C=\max_{i:\mu_{i}(B(x,r))>0}C_{i}$ and $d=\max_{i:\mu_{i}(B(x,r))>0}d_{i}$ and we used the fact that if $\mu_{i}(B(x,r))=0$ then $\mu_{i}(B(x,\theta r))=0$. Therefore, $\mu$ is $(C,d)$-homogeneous on $B(x,r)$.

This result applies to the case of $d$-sparse data and is explained in the following example.

Example 3** ($d$-sparse data).**

Suppose that $\mathcal{X}\subset\mathbb{R}^{D}$ is defined as

[TABLE]

Let $\mu$ be a probability measure on $\mathcal{X}$. In this case, we can write $\mathcal{X}$ as the union of $k=\binom{D}{d},d$-dimensonal hyperplanes in $\mathbb{R}^{D}$. In fact,

[TABLE]

Letting $\mu_{i_{1},i_{2},\ldots,i_{d}}$ be the probability measure restricted to the hyperplane defined by $x_{j}=0,j\not\in\left\{i_{1},i_{2},\ldots,i_{d}\right\}$, we can express $\mu=\sum_{1\leq i_{1}<i_{2}<\cdots i_{d}\leq D}\pi_{i_{1},i_{2},\ldots,i_{d}}\mu_{i_{1},i_{2},\ldots,i_{d}}$. Therefore, the result of Example 2 implies that for any $x\in\mathcal{X}$, for $r$ that is small enough $\mu$ is $(C,d)$-homogeneous on $B(x,r)$.

Our final example is about the product measure. This allows us to prove that any concatenation of spaces with small intrinsic dimension has a small intrinsic dimension as well.

Example 4** (Concatenation under the product measure).**

Suppose that $\mu_{i}$ is a probability measure on $\mathcal{X}_{i}\subset\mathbb{R}^{D_{1}}$. Define $\mathcal{X}=\left\{(z_{1},z_{2})\mid z_{1}\in\mathcal{X}_{1},z_{2}\in\mathcal{X}_{2}\right\}$ and let $\mu=\mu_{1}\times\mu_{2}$ be the product measure on $\mathcal{X}$, i.e., $\mu(E_{1}\times E_{2})=\mu_{1}(E_{1})\times\mu_{2}(E_{2})$ for $E_{i}$ that is $\mu_{i}$-measurable, $i=1,2$. Suppose that $\mu_{i}$ is $(C_{i},d_{i})$-homogeneous on $B(x_{i},r_{i})$ and let $x=(x_{1},x_{2})$. Then, $\mu$ is $(C,d)$-homogeneous on $B(x,r)$, where $d=d_{1}+d_{2},r=\min\left\{r_{1},r_{2}\right\}$ and $C=(C_{1}\,C_{2}\,r^{-(d_{1}+d_{2})}\,2^{(d_{1}+d_{2})/2})/(r_{1}^{-d_{1}}\,r_{2}^{-d_{2}})$. To establish this, let $r=\min\left\{r_{1},r_{2}\right\}$ and note that for any $\theta\in(0,1)$ we have

[TABLE]

where we used two simple inequalities that $\|(z_{1},z_{2})-(x_{1},x_{2})\|_{2}\leq r$ implies $\|z_{i}-x_{i}\|_{2}\leq r,~{}i=1,2,$ and further $\|z_{i}-x_{i}\|_{2}\leq r/\sqrt{2},~{}i=1,2,$ implies $\|(z_{1},z_{2})-(x_{1},x_{2})\|_{2}\leq r$.

Appendix C Simulation Setting

Here we explain the settings for simulations shown in Figures 1 and 2.

C.1 Single test point

The data for single test point simulation, shown in Figure 1, has been generated as follows. Here $p=1$, $D=20$ and $d=2$. All the points are generated using $X_{i}=AX_{i}^{\text{low}}$, where $A\in\mathbb{R}^{D\times d}$ and entries of $A$ are independently sampled from $U[-1,1]$. Components of each $X_{i}^{\text{low}}$ are also generated independently from $U[-1,1]$. We generate a fix test point $x_{\text{test}}=Ax_{\text{test}}^{\text{low}}$ and keep the matrix $A$ throughout all Monte-Carlo iterations fixed. In each Monte-Carlo iteration, we generate $n=20000$ training points as mentioned before. The values of $Y_{i}$ are generated according to $Y_{i}=f(X_{i})+\varepsilon_{i}$, where $f(X)=\frac{1}{1+\exp(-3X[0])}$, and $\varepsilon_{i}\sim\mathsf{N}(0,\sigma_{e}^{2})$ with $\sigma_{e}=1$. We are interested in estimate and inference for $f(x_{\text{test}})$ which is equivalent to solving for $\mathbf{E}[\psi(Z;\theta(x))\mid X=x]=0$ with $\psi(Z;\theta(x))=Y-\theta(x)$ at $x=x_{\text{test}}$. We run DNN (Algorithm 2) for $k=1,2$ and $5$ with parameter $s=s_{\zeta}$ chosen using Proposition 2 with $\zeta=0.1$ over $1000$ Monte-Carlo iterations and report the histogram and quantile-quantile plot of estimates compared to theoretical asymptotic normal distribution of estimates stemming from our characterization. In our simulations, we considered the complete $U$-statistic case, i.e., $B=\binom{n}{s}$.

C.2 Multiple test point

The data for the multiple test point simulation, shown in Figure 2, has been generated very similarly to the single test point setting. The only difference is that instead of generating a single test point we generate $100$ test points. These test points together with matrix $A$ are kept fixed throughout all $1000$ Monte-Carlo iterations. We compare the performance of DNN (Algorithm 2) with parameter $s=s_{\zeta}$ chosen using Proposition 2 with $\zeta=0.1$ with two benchmarks that set $s_{d}=n^{1.05d/(d+2)}$ and $s_{D}=n^{1.05D/(D+2)}$. This process has been repeated for $k=1,2$ and $5$ and the coverage over a single run for all test points, the empirical coverage over $1000$ runs, and chosen $s_{\zeta}$ versus $s_{d}$ are depicted.

Appendix D Nuisance parameters and heterogeneous treatment effects

Using the techniques of Oprescu et al. [2018], our work also easily extends to the case where the moments depend on, potentially infinite dimensional, nuisance components $h_{0}$, that also need to be estimated, i.e.,

[TABLE]

If the moment $m$ is orthogonal with respect to $h$ and assuming that $h_{0}$ can be estimated on a separate sample with a conditional MSE rate of

[TABLE]

then using the techniques of Oprescu et al. [2018], we can argue that both our finite sample estimation rate and our asymptotic normality rate, remain unchanged, as the estimation error only impacts lower order terms. This extension allows us to capture settings like heterogeneous treatment effects, where the treatment model also needs to be estimated when using the orthogonal moment as

[TABLE]

where $y$ is the outcome of interest, $t$ is a treatment, $x,w$ are confounding variables, $q_{0}(x,w)=\mathbf{E}[Y|X=x,W=w]$ and $p_{0}(x,w)=E[T|X=x,W=w]$. The latter two nuisance functions can be estimated via separate non-parametric regressions. In particular, if we assume that these functions are sparse linear in $w$, i.e.:

[TABLE]

Then we can achieve a conditional mean-squared-error rate of the required order by using the kernel lasso estimator of Oprescu et al. [2018], where the kernel is the sub-sampled $k$-NN kernel, assuming the sparsity does not grow fast with $n$.

Appendix E Proof of Theorem 1

Lemma 4**.**

For any $\theta\in\Theta$:

[TABLE]

Proof.

By strong convexity of the loss $L(x;\theta)$ and the fact that $m(x;\theta(x))=0$, we have:

[TABLE]

By convexity of the loss $L(x;\theta)$ we have:

[TABLE]

Combining the latter two inequalities we get:

[TABLE]

Note that if $\|\theta-\theta(x)\|_{2}=0$, then the result is obvious. Otherwise, dividing over by $\|\theta-\theta(x)\|_{2}$ completes the proof of the lemma. ∎

Lemma 5**.**

Let $\Lambda(x;\theta)=m(x;\theta)-\Psi(x;\theta)$. Then the estimate $\hat{\theta}$ satisfies:

[TABLE]

Proof.

Observe that $\hat{\theta}$, by definition, satisfies $\Psi(x;\hat{\theta})=0$. Thus:

[TABLE]

∎

Lemma 6**.**

Suppose that the kernel is built with sub-sampling at rate $s$, in an honest manner (Assumption 3) and with at least $B\geq n/s$ sub-samples. If the base kernel satisfies kernel shrinkage in expectation, with rate $\epsilon(s)$, then w.p. $1-\delta$:

[TABLE]

Proof.

Define

[TABLE]

where we remind that $\Psi_{0}$ denotes the complete $U$-statistic:

[TABLE]

Here the expectation is taken with respect to the random draws of $n$ samples. Then, the following result which is due to Oprescu et al. [2018] holds.

Lemma 7** (Adapted from Oprescu et al. [2018]).**

For any $\theta$ and target $x$

[TABLE]

In other words, Lemma 7 states that, in the expression for $\mu_{0}$ we can simply replace $\psi(Z_{i};\theta)$ with its expectation which is $m(X_{i};\theta)$. We can then express $\Lambda(x;\theta)$ as sum of kernel error, sampling error, and sub-sampling error, by adding and subtracting appropriate terms, as follows:

[TABLE]

The parameters should be chosen to trade-off these error terms nicely. We will now bound each of these three terms separately and then combine them to get the final bound.

Bounding the Kernel error.

By Lipschitzness of $m$ with respect to $x$ and triangle inequality, we have:

[TABLE]

where the second to last inequality follows from the fact that $\sum_{i}|\alpha_{S_{b}}(X_{i})|=1$.

Bounding the Sampling error.

For bounding the sampling error we rely on Lemma 15 and in particular Corollary 4. Observe that for each $j\in\{1,\ldots,p\}$, $\Psi_{0j}(x;\theta)$ is a complete $U$-statistic for each $\theta$. Thus the sampling error defines a $U$-process over the class of symmetric functions $\text{conv}(\mathcal{F}_{j})=\{f_{j}(\cdot;\theta):\theta\in\Theta\}$, with $f_{j}(Z_{1},\ldots,Z_{s};\theta)=\mathbf{E}_{\omega}\left[\sum_{i=1}^{s}\alpha_{Z_{1:s},\omega}(X_{i})\psi_{j}(Z_{i};\theta)\right]$. Observe that since $f_{j}\in\text{conv}(\mathcal{F}_{j})$ is a convex combination of functions in $\mathcal{F}_{j}=\{\psi_{j}(\cdot;\theta):\theta\in\Theta\}$, the bracketing number of functions in $\text{conv}(\mathcal{F}_{j})$ is upper bounded by the bracketing number of $\mathcal{F}_{j}$, which by our assumption, satisfies $\log(N_{[]}({\cal F}_{j},\epsilon,L_{2}))=O(1/\epsilon)$. Moreover, by our assumptions on the upper bound $\psi_{\max}$ of $\psi_{j}(z;\theta)$, we have that $\sup_{f_{j}\in\text{conv}(\mathcal{F}_{j})}\|f_{j}\|_{2},\sup_{f_{j}\in\text{conv}(\mathcal{F}_{j})}\|f_{j}\|_{\infty}\leq\psi_{\max}$. Thus all conditions of Corollary 4 are satisfied, with $\eta=G=\psi_{\max}$ and we get that w.p. $1-\delta/2p$:

[TABLE]

By a union bound over $j$, we get that w.p. $1-\delta/2$:

[TABLE]

Bounding the Sub-sampling error.

Sub-sampling error decays as $B$ is increased. Note that for a fixed set of samples $\left\{Z_{1},Z_{2},\ldots,Z_{n}\right\}$, for a set $S_{b}$ randomly chosen among all $\binom{n}{s}$ subsets of size $s$ from the $n$ samples, we have:

[TABLE]

Therefore, $\Psi(x;\theta)$ can be thought as the sum of $B$ i.i.d. random variables each with expectation equal to $\Psi_{0}(x;\theta)$, where expectation is taken over $B$ draws of sub-samples, each with size $s$. Thus one can invoke standard results on empirical processes for function classes as a function of the bracketing entropy. For simplicity, we can simply invoke Corollary 4 in the appendix for the case of a trivial $U$-process, with $s=1$ and $n=B$ to get that w.p. $1-\delta/2$:

[TABLE]

Thus for $B\geq n/s$, the sub-sampling error is of lower order than the sampling error and can be asymptotically ignored. Putting together the upper bounds on sampling, sub-sampling and kernel error finishes the proof of the Lemma. ∎

The probabilistic statement of the proof follows by combining the inequalities in the above three lemmas. The in expectation statement follows by simply integrating the exponential tail bound of the probabilistic statement.

Appendix F Proof of Theorem 2

We will show asymptotic normality of $\hat{\alpha}=\left\langle\beta,\hat{\theta}\right\rangle$ for some arbitrary direction $\beta\in\mathbb{R}^{p}$, with $\|\beta\|_{2}\leq R$. Consider the complete multi-dimensional $U$-statistic:

[TABLE]

Let

[TABLE]

where $\mu_{0}(x;\theta)=\mathbf{E}\left[\Psi_{0}(x;\theta)\right]$ (as in the proof of Theorem 1) and

[TABLE]

Finally, let

[TABLE]

For shorthand notation let $\alpha_{0}=\left\langle\beta,\theta(x)\right\rangle$, $\psi_{\beta}(Z;\theta)=\left\langle\beta,M_{0}^{-1}(\psi(Z;\theta)-m(X;\theta))\right\rangle$ and

[TABLE]

be a single dimensional complete $U$-statistic. Thus we can re-write:

[TABLE]

We then have the following lemma which its proof is provided in Appendix J:

Lemma 8**.**

Under the conditions of Theorem 2:

[TABLE]

for $\sigma_{n}^{2}(x)=\frac{s^{2}}{n}\operatorname{Var}\left[\mathbf{E}\left[\sum_{i=1}^{s}K(x,X_{i},\{X_{j}\}_{j=1}^{s})\psi_{\beta}(Z_{i};\theta)\mid X_{1}\right]\right]=\Omega(\frac{s^{2}}{n}\eta(s))$.

Invoking Lemma 8 and using our assumptions on the kernel, we conclude that:

[TABLE]

For some sequence $\sigma_{n}^{2}$ which decays at least as slow as $s^{2}\eta(s)/n$. Hence, since

[TABLE]

if we show that $\frac{\hat{\alpha}-\tilde{\alpha}}{\sigma_{n}(x)}\rightarrow_{p}0$, then by Slutsky’s theorem we also have that:

[TABLE]

as desired. Thus, it suffices to show that:

[TABLE]

Observe that since $\|\beta\|_{2}\leq R$, we have $\|\hat{\alpha}-\tilde{\alpha}\|_{2}\leq R\|\hat{\theta}-\tilde{\theta}\|_{2}$. Thus it suffices to show that:

[TABLE]

Lemma 9**.**

Under the conditions of Theorem 2, for $\sigma_{n}^{2}(x)=\Omega\left(\frac{s^{2}}{n}\eta(s)\right)$:

[TABLE]

Proof.

Performing a second-order Taylor expansion of $m_{j}(x;\theta)$ around $\theta(x)$ and observing that $m_{j}(x;\theta(x))=0$, we have that for some $\bar{\theta}_{j}\in\Theta$:

[TABLE]

Letting $\rho=(\rho_{1},\ldots,\rho_{p})$, writing the latter set of equalities for each $j$ in matrix form, multiplying both sides by $M_{0}^{-1}$ and re-arranging, we get that:

[TABLE]

Thus by the definition of $\tilde{\theta}$ we have:

[TABLE]

By the bounds on the eigenvalues of $H_{j}(x;\theta)$ and $M_{0}^{-1}$, we have that:

[TABLE]

Thus we have:

[TABLE]

By our estimation error Theorem 1, we have that the expected value of the second term on the right hand side is of order $O\left(\epsilon(s)^{2},\frac{s}{n}\log\log(n/s)\right)$. Thus by the assumptions of the theorem, both are $o(\sigma_{n})$. Hence, the second term is $o_{p}(\sigma_{n})$.

We now argue about the convergence rate of the first term on the right hand side. Similar to the proof of Theorem 1, since $\Psi(x;\hat{\theta})=0$ we have:

[TABLE]

We can further add and subtract $\mu_{0}$ from $m(x;\hat{\theta})$.

[TABLE]

Combining we have:

[TABLE]

Now similar to proof of Theorem 1 we bound different terms separately and combine the results.

Kernel Error.

Term $C$ is a kernel error and hence is upper bounded by $\epsilon(s)$ in expectation. Since, by assumption $s$ is chosen such that $\epsilon(s)=o(\sigma_{n}(x))$, we ge that $\|C\|_{2}/\sigma_{n}(x)\rightarrow_{p}0$.

Sub-sampling Error.

Term $E$ is a sub-sampling error, which can be made arbitrarily small if the number of drawn sub-samples is large enough and hence $\|E\|_{2}/\sigma_{n}(x)\rightarrow_{p}0$. In fact, similar to the part about bounding sub-sampling error in Lemma 6 we have that that:

[TABLE]

Therefore, $\Psi(x;\theta)$ can be thought as the sum of $B$ independent random variables each with expectation equal to $\Psi_{0}(x;\theta)$. Now we can invoke Corollary 4 in the appendix for the trivial U-process, with $s=1,n=B$ to get that w.p. $1-\delta_{1}$:

[TABLE]

Hence, for $B\geq(n/s)^{5/4}$, due to our assumption that $\left(s/n\log\log(n/s)\right)^{5/8}=o(\sigma_{n}(x))$ we get $\|E\|_{2}/\sigma_{n}(x)\rightarrow_{p}0$.

Sampling Error.

Thus it suffices that show that $\|F\|_{2}/\sigma_{n}(x)\rightarrow_{p}0$, in order to conclude that $\frac{\|m(x;\hat{\theta})+\Psi_{0}(x;\theta(x))\|_{2}}{\sigma_{n}(x)}\rightarrow_{p}0$. Term $F$ can be re-written as:

[TABLE]

Observe that each coordinate $j$ of $F$, is a stochastic equicontinuity term for $U$-processes over the class of symmetric functions $\text{conv}(\mathcal{F}_{j})=\{f_{j}(\cdot;\theta):\theta\in\Theta\}$, with $f_{j}(Z_{1},\ldots,Z_{s};\theta)=\mathbf{E}_{\omega}\left[\sum_{i=1}^{s}\alpha_{Z_{1:s},\omega}(X_{i})(\psi_{j}(Z_{i};\theta(x))-\psi_{j}(Z_{i};\theta))\right]$. Observe that since $f_{j}\in\text{conv}(\mathcal{F}_{j})$ is a convex combination of functions in $\mathcal{F}_{j}=\{\psi_{j}(\cdot;\theta(x))-\psi_{j}(\cdot;\theta):\theta\in\Theta\}$, the bracketing number of functions in $\text{conv}(\mathcal{F}_{j})$ is upper bounded by the bracketing number of $\mathcal{F}_{j}$, which in turn is upper bounded by the bracketing number of the function class $\{\psi_{j}(\cdot;\theta):\theta\in\Theta\}$, which by our assumption, satisfies $\log(N_{[]}({\cal F}_{j},\epsilon,L_{2}))=O(1/\epsilon)$. Moreover, under the variogram assumption and the lipschitz moment assumption we have that if $\|\theta-\theta(x)\|\leq r\leq 1$, then:

[TABLE]

Moreover, $\|f_{j}\|_{\infty}\leq 2\psi_{\max}$. Thus we can apply Corollary 4, with $\eta=\sqrt{L_{\psi}r+L_{J}^{2}r^{2}}=O(\sqrt{r})$ and $G=2\psi_{\max}$ to get that if $\|\hat{\theta}-\theta(x)\|\leq r$, then w.p. $1-\delta/p$:

[TABLE]

Using a union bound this implies that w.p. $1-\delta$ we have

[TABLE]

By our MSE theorem and also Markov’s inequality, w.p. $1-\delta^{\prime}$: $\|\hat{\theta}-\theta(x)\|\leq\nu(s)/\delta^{\prime}$, where:

[TABLE]

Thus using a union bound w.p. $1-\delta-\delta^{\prime}$, we have:

[TABLE]

To improve readability from here we ignore all the constants in our analysis, while we keep all terms (even $\log$ or $\log\log$ terms) that depend on $s$ and $n$. Note that we can even ignore $\delta$ and $\delta^{\prime}$, because they can go to zero at very slow rate such that terms $\log(1/\delta)$ or even $\delta^{\prime 1/4}$ appearing in the analysis grow slower than $\log\log$ terms. Now, by the definition of $\nu(s)$ and $\kappa(r,s,n,\delta^{\prime})$, as well as invoking the inequality $(a+b)^{1/4}\leq a^{1/4}+b^{1/4}$ for $a,b>0$ we have:

[TABLE]

Hence, using our Assumption on the rates in the statement of Theorem 2 we get that both of the terms above are $o(\sigma_{n}(x))$. Therefore, $\|F\|_{2}/\sigma_{n}(x)\rightarrow_{p}0$. Thus, combining all of the above, we get that:

[TABLE]

as desired. ∎

Remark 1**.**

Our notion of incrementality is slightly different from that of Wager and Athey [2018], as there the incrementality is defined as $\operatorname{Var}\left[\mathbf{E}\left[K(x,X_{1},\{Z_{j}\}_{j=1}^{s})\mid X_{1}\right]\right]$. However, using the tower law of expectation

[TABLE]

For a symmetric kernel the term $\mathbf{E}\left[K(x,X_{1},\{Z_{j}\}_{j=1}^{s})\right]^{2}$ is equal to $1/s^{2}$ and is asymptotically negligible compared to $\operatorname{Var}\left[\mathbf{E}\left[K(x,X_{1},\{Z_{j}\}_{j=1}^{s})\mid X_{1}\right]\right]$, which usually decays at a slower rate.

Appendix G Lower Bound on Incrementality as Function of Kernel Shrinkage

We give a generic lower bound on the quantity $\mathbf{E}[\mathbf{E}[K(x,X_{1},\{Z_{j}\}_{j=1}^{s})|X_{1}]^{2}]$ that depends only on the Kernel shrinkage. The bound essentially implies that if we know that the probability that the distribution of $x$’s assigns to a ball of radius $\epsilon(s,1/2s)$ around the target $x$ is of order $1/s$, i.e. we should expect at most a constant number of samples to fall in the kernel shrinkage ball, then the main condition on incrementality of the kernel, required for asymptotic normality, holds. In some sense, this property states that the kernel shrinkage behavior is tight in the following sense. Suppose that the kernel was assigning positive weight to at most a constant number of $k$ samples. Then kernel shrinkage property states that with high probability we expect to see at least $k$ samples in a ball of radius $\epsilon(s,\delta)$ around $x$. The above assumption says that we should also not expect to see too many samples in that radius, i.e. we should also expect to see at most a constant number $K>k$ of samples in that radius. Typically, the latter should hold, if the characterization of $\epsilon(s,\delta)$ is tight, in the sense that if we expected to see too many samples in the radius, then most probably we could have improved our analysis on Kernel shrinkage and given a better bound that shrinks faster.

G.1 Proof of Lemma 1

By the Paley-Zygmund inequality, for any random variable $Z\geq 0$ and for any $\delta\in[0,1]$:

[TABLE]

Let $W_{1}=K(x,X_{1},\{Z_{j}\}_{j=1}^{s})$. Then, applying the latter to the random variable $Z=\mathbf{E}[W_{1}|X_{1}]$ and observing that by symmetry $\mathbf{E}[Z]=\mathbf{E}[W_{1}]=1/s$, yields:

[TABLE]

Moreover, observe that by the definition of $\epsilon(s,\rho)$ for some $\rho>0$:

[TABLE]

This means that at most a mass $\rho\,s/\delta$ of the support of $X_{1}$ in the region $\|X_{1}-x\|\geq\epsilon(s,\rho)$ can have $\Pr[W_{1}>0|X_{1}]\geq\delta/s$. Otherwise the overall probability that $W_{1}>0$ in the region of $\|X_{1}-x\|\geq\epsilon(s,\rho)$ would be more than $\rho$. Thus we have that except for a region of mass $\rho\,s/\delta$, for each $X_{1}$ in the region $\|X_{1}-x\|\geq\epsilon(s,\rho)$: $\mathbf{E}[W_{1}|X_{1}]\leq\delta/s$. Combining the above we get:

[TABLE]

Thus:

[TABLE]

Since $\rho$ was arbitrarily chosen, the latter upper bound holds for any $\rho$, which yields the result.

G.2 Proof of Corollary 1

Thus applying Lemma 1 with $\delta=1/2$ yields:

[TABLE]

Observe that:

[TABLE]

Hence:

[TABLE]

where the last follows by choosing $\rho=1/s^{2}$. Combining all the above yields:

[TABLE]

Appendix H Proofs of Section 5

H.1 Proof of Lemma 2

For proving this result, we rely on Bernstein’s inequality which is stated below:

Proposition 3** (Bernstein’s Inequality).**

Suppose that random variables $Z_{1},Z_{2},\ldots,Z_{n}$ are i.i.d., belong to $[-c,c]$ and $\mathbf{E}[Z_{i}]=\mu$. Let $\bar{Z}_{n}=\frac{1}{n}\sum_{i=1}^{n}Z_{i}$ and $\sigma^{2}=\operatorname{Var}(Z_{i})$. Then, for any $\theta>0$,

[TABLE]

This also implies that w.p. at least $1-\delta$ the following holds:

[TABLE]

Let $A$ be any $\mu$-measurable set. An immediate application of Bernstein’s inequality to random variables $Z_{i}=\,1\left\{X_{i}\in A\right\}$, implies that w.p. $1-\delta$ over the choice of covariates $\left(X_{i}\right)_{i=1}^{s}$, we have:

[TABLE]

In above, we used the fact that $\operatorname{Var}(Z_{i})=\mu(A)(1-\mu(A))\leq\mu(A)$. This result has the following corollary.

Corollary 3**.**

Define $U=2\log(2/\delta)/s$ and let $A$ be an arbitrary $\mu$-measurable set. Then, w.p. $1-\delta$ over the choice of training samples, $\mu(A)\geq 4U$ implies $\mu_{s}(A)\geq U$.

Proof.

Define $U=2\log(2/\delta)/s$. Then, Bernstein’s inequality in Proposition 3 implies that w.p. $1-\delta$ we have

[TABLE]

Assume that $\mu(A)\geq 4U$, we want to prove that $\mu_{s}(A)\geq U$. Suppose, the contrary, i.e., $\mu_{s}(A)<U$. Then, by dividing the above equation by $\mu(A)$ we get

[TABLE]

Note that since $\mu_{s}(A)<U<\mu(A)$, by letting $z=U/\mu(A)\leq 1/4$ the above implies that

[TABLE]

which as $z>0$ only holds for

[TABLE]

This contradicts with $z\leq 1/4$, implying the result. ∎

Now we are ready to finish the proof of Lemma 2. First, note that using the definition of $(C,d)$-homogeneous measure. Note that for any $\theta\in(0,1)$ we have $\mu(B(x,\theta r))\geq(1/C)\theta^{d}\mu(B(x,r))$. Replace $\theta r=\epsilon$ in above. It implies that for any $\epsilon\in(0,r)$

[TABLE]

Pick $\epsilon_{k}(s,\delta)$ according to

[TABLE]

Note that for having $\epsilon_{k}(s,\delta)\in(0,r)$ we need

[TABLE]

Therefore, replacing this choice of $\epsilon_{k}(s,\delta)$ in Equation (35) implies that $\mu(B(x,\epsilon_{k}(s,\delta)))\geq\frac{8\log(2/\delta)}{s}$. Now we can use the result of Corollary 3 for the choice $A=B(x,\epsilon_{k}(s,\delta))$. It implies that w.p. $1-\delta$ over the choice of $s$ training samples, we have

[TABLE]

Note that whenever $\delta\leq\exp(-k/2)/2$ we have

[TABLE]

Therefore, w.p. $1-\delta$ we have

[TABLE]

H.2 Proof of Corollary 2

Lemma 2 shows that for any $t=\epsilon_{k}(s,\delta)=r\left(\frac{8C\log(2/\delta)}{\mu(B(x,r))s}\right)^{1/d}$, such that $t\leq r$ and $t\geq r\left(\frac{4k\,C}{\mu(B(x,r))s}\right)^{1/d}$, we have that:

[TABLE]

Let $\rho=\frac{1}{r}\left(\frac{\mu(B(x,r))}{8C}\right)^{1/d}$, which is a constant. Solving for $\delta$ in terms of $t$ we get:

[TABLE]

for any $t\in\left[\frac{(s/k)^{-1/d}}{\rho},r\right]$. Thus, noting that $X_{i}$’s and target $x$ both belong to $\mathcal{X}$ that has diameter $\Delta_{\mathcal{X}}$, we can upper bound the expected value of $[\|x-X_{(k)}\|_{2}$ as:

[TABLE]

Note that for $s$ larger than some constant, we have $\exp\left\{-\rho^{d}\,r^{d}\,s\right\}\leq s^{-1/d}$. Thus the first and last terms in the latter summation are of order $\left(\frac{1}{s}\right)^{1/d}$. We now show that the same holds for the middle term, which would complete the proof. By setting $u=\rho^{d}\,s\,t^{d}$ and doing a change of variables in the integral we get:

[TABLE]

where $\Gamma$ is the Gamma function. Since by the properties of the Gamma function $z\Gamma(z)=\Gamma(z+1)$, the latter evaluates to: $\frac{s^{-1/d}}{\rho}\Gamma((d+1)/d)$. Since $(d+1)/d\in[1,2]$, we have that $\Gamma((d+1)/d)\leq 2$. Thus the middle term is upper bounded by $\frac{2s^{-1/d}}{\rho}$, which is also of order $\left(\frac{1}{s}\right)^{1/d}$.

H.3 Proof of Lemma 3

Before proving this lemma we state and prove and auxiliary lemma which comes in handy in our proof.

Lemma 10**.**

Let $P_{1}$ denote the mass that the density of the distribution of $X_{i}$ puts on the ball around $x$ with radius $\|x-X_{1}\|_{2}$, which is a random variable as it depends on $X_{1}$. Then, for any $s\geq k$ the following holds:

[TABLE]

Proof.

The proof is an easy consequence of symmetry. Let S_{1}=\,1\{\text{sample 1 is among k nearest neighbors}\}, then we can write

[TABLE]

which simply computes the probability that there are at most $k-1$ other points in the ball with radius $\|x-X_{1}\|$. Now, by using the tower law

[TABLE]

which holds because of the symmetry. In other words, the probability that sample $1$ is among the $k$-NN is equal to $k/s$. Hence, the conclusion follows. ∎

We can finish the proof of Lemma 3. Define S_{1}=\,1\{\text{sample 1 is among k nearest neighbors}\}, then we can write

[TABLE]

Recall that if $P_{1}$ denotes the mass that the density of the distribution of $X_{i}$ puts on the ball around $x$ with radius $\|x-X_{1}\|_{2}$, which is a random variable depending on $X_{1}$. Therefore,

[TABLE]

Now we can write

[TABLE]

Now using Lemma 10 (where $s$ is replaced by $2s-1$) we know that for any value of $0\leq r\leq 2s-2$ we have

[TABLE]

This implies that for any value of $r$ we have $\mathbf{E}\left[b_{r}(1-P_{1})^{2s-2-r}P_{1}^{r}\right]=1/(2s-1)$. The reason is simple. Note that the above is obvious for $r=0$ using Equation (36). For other values of $r\geq 1$, we can write Equation (36) for values $r$ and $r-1$. Taking their difference implies the result. Note that this further implies that $\mathbf{E}\left[(1-P_{1})^{2s-2-r}P_{1}^{r}\right]=1/(b_{r}\,(2s-1))$, as $b_{r}$ is a constant. Therefore, by plugging this back into the expression of $\mathbf{E}[\mathbf{E}[S_{1}\mid X_{1}]^{2}]$ we have

[TABLE]

which implies the desired result.

Remark 2**.**

Note that $b_{t}=\binom{2s-2}{t}$ since we can view $b_{t}$ as follows: how many different subsets of size $t$ can we create from a set of $2s-2$ elements if we pick a number $i=\{0,\ldots,t\}$ and then choose $i$ elements from the first half of these elements and $t-i$ elements from the second half. Observe that this process creates all possible sets of size $t$ from among the $2s-2$ elements, which is equal to $\binom{2s-2}{t}$.

Furthermore, $a_{t}=b_{t}$ for $0\leq t\leq k-1$ and for any $k\leq t\leq 2k-2$, after some little algebra, we have

[TABLE]

This implies that the summation appeared in Lemma 11 satisfies

[TABLE]

H.4 Proof of Theorem 4

Note that according to Lemma 17, the asymptotic variance $\sigma_{n,j}^{2}(x)=\frac{s^{2}}{n}\operatorname{Var}\left[\Phi_{1}(Z_{1})\right]$, where $\Phi_{1}(Z_{1})=\frac{1}{k}\mathbf{E}[\sum_{i\in H_{k}(x,s)}\left\langle e_{j},M_{0}^{-1}\psi(Z_{i};\theta(x))\right\rangle\mid Z_{1}]$. Therefore, once we establish an expression for $\operatorname{Var}\left[\Phi_{1}(Z_{1})\right]$ we can finish the proof of this theorem. The following lemma provides such result.

Lemma 11**.**

Suppose that the kernel $K$ is the $k$-NN kernel and let $\sigma_{j}^{2}(x)=\operatorname{Var}\left(\left\langle e_{j},M_{0}^{-1}\psi(z;\theta(x))\right\rangle\mid X=x\right)$. Moreover, suppose that $\epsilon_{k}(s,1/s^{2})\rightarrow 0$ for any constant $k$. Then:

[TABLE]

where the second equality above holds due to Lemma 3 and sequences $a_{t}$ and $b_{t}$, for $0\leq t\leq 2k-2$, are defined in Lemma 3.

Proof.

In this proof for simplicity we let $Y_{i}=\left\langle e_{j},M_{0}^{-1}\psi(Z_{i};\theta(x))\right\rangle$ and $\mu(X_{i})=\mathbf{E}[Y_{i}]=\left\langle e_{j},M_{0}^{-1}m(X_{i};\theta(x))\right\rangle$. Let $Z_{(i)}$ denote the random variable of the $i$-th closest sample to $x$. For the case of $k$-NN we have that:

[TABLE]

Let S_{1}=\mathbf{1}\{\text{sample 1 is among k nearest neighbors}\}. Then we have:

[TABLE]

Let $\tilde{Y}_{(i)}$ denote the label of the $i$-th closest point to $x$, excluding sample $1$. Then:

[TABLE]

Observe that $\tilde{Y}_{(i)}$ are all independent of $Z_{1}$. Hence:

[TABLE]

Therefore the variance of $\Phi(Z_{1})$ is equal to the variance of the first term on the right hand side. Hence:

[TABLE]

Where we used the fact that:

[TABLE]

Moreover, observe that under the event that $S_{1}=1$, we know that the difference between the closest $k$ values and the closest $k$ values excluding $1$ is equal to the difference between the $Y_{1}$ and $Y_{(k+1)}$. Hence:

[TABLE]

where the last equation holds from the fact that for any $j\neq 1$, conditional on $X_{j}$, the random variable $Y_{j}$ is independent of $Z_{1}$ and is equal to $\mu(X_{j})$ in expectation. Under the event $S_{1}=1$, we know that the $(k+1)$-th closest point is different from sample $1$. We now argue that up to lower order terms, we can replace $\mu(X_{(k+1)})$ with $\mu(X_{1})$ in the last equality:

[TABLE]

Observe that:

[TABLE]

Moreover, by Jensen’s inequality, Lipschitzness of the first moments and kernel shrinkage:

[TABLE]

Hence, for $\delta=1/s^{2}$, the latter is $o(1/s)$. Similarly:

[TABLE]

which for $\delta=1/s^{2}$ is also of order $o(1/s)$. Combining all the above we thus have:

[TABLE]

We now work with the first term on the right hand side. By the tower law of expectations:

[TABLE]

By Lipschitzness of the second moments, we know that the second part is upper bounded as:

[TABLE]

For $\delta=1/s^{2}$ it is of $o(1/s)$. Thus:

[TABLE]

Note that Lemma 3 provides an expression for $\mathbf{E}\left[\mathbf{E}[S_{1}\mid X_{1}]^{2}\right]$ which finishes the proof. ∎

For finishing proof of Theorem 4 we need to prove that $\sum_{t=0}^{2k-2}\frac{a_{t}}{b_{t}}$ is equal to $\zeta_{k}$ plus lower order terms. This is proved in the following lemma.

Lemma 12**.**

Suppose that $s\rightarrow\infty$ and $k$ is fixed. Then

[TABLE]

Proof.

Note that for any $0\leq t\leq k-1$ we have $a_{t}=b_{t}$ according to Remark 2. For any $k\leq t\leq 2k-2$ we have

[TABLE]

where we used the fact that $t$ and $i$ are both bounded above by $2k-2$ which is a constant. Hence,

[TABLE]

as desired. ∎

H.5 Proof of Theorem 5

The goal is to apply Theorem 2. Note that $k$-NN kernel is both honest and symmetric. According to Lemma 2, we have that $\epsilon_{k}(s,\delta)=O\left((\log(1/\delta)/s)^{1/d}\right)$ for $\exp(-Cs)\leq\delta\leq D$, where $C$ and $D$ are constants. Corollary 2 also implies that $\epsilon_{k}(s)=O((1/s)^{1/d})$. Furthermore, according to Lemma 3, the incrementality $\eta_{k}(s)$ is $\Theta(1/s)$. Therefore, as $s$ goes to $\infty$ we have $\epsilon_{k}(s,\eta_{k}(s))=O\left((\log(s)/s)^{1/d}\right)\rightarrow 0$. Moreover, as $\eta_{k}(s)=\Theta(1/s)$, we also get that $n\eta_{k}(s)=O(n/s)\rightarrow\infty$. We only need to ensure that Equation (9) is satisfied. Note that $\sigma_{n,j}(x)=\Theta(\sqrt{s/n})$. Therefore, by dividing terms in Equation (9) it suffices that

[TABLE]

Note that due to our Assumption $n/s\rightarrow\infty$, the last term obviously goes to zero. Also, because of the assumption made in the statement of theorem, the first term also goes to zero. We claim that if the first term goes to zero, the same also holds for the second term. Note that we can write

[TABLE]

and since $n/s\rightarrow\infty$, our claim follows. Therefore, all the conditions of Theorem 2 are satisfied and the result follows.

The second part of result is implied by the first part since if $s=n^{\beta}$ and $\beta\in(d/(d+2),1)$ then $s^{-1/d}\sqrt{\frac{n}{s}}\rightarrow 0$.

H.6 Proof of Proposition 1

For proving this lemma, we need two following auxiliary results. Before that we state the Hoeffding’s inequality for $U$-statistics Hoeffding [1994].

Proposition 4**.**

Suppose that $\mathbf{X}=\left(X_{1},X_{2},\ldots,X_{n}\right)$ are i.i.d. and $q$ is a function that has range $[0,1]$. Define $U_{s}=\binom{n}{s}^{-1}\sum_{i_{1}<i_{2}<\ldots<i_{s}}q(X_{i_{1}},X_{i_{2}},\ldots,X_{i_{s}})$. Then, for any $\epsilon>0$

[TABLE]

Furthermore, for any $\delta>0$, w.p. $1-\delta$ we have

[TABLE]

Lemma 13**.**

Consider the function $H(s)$ defined in Section 5.3 and $G_{\delta}(s)=\Delta\sqrt{2ps/n\log(2np/\delta)}$. Then, w.p. $1-\delta$, for all values of $k\leq s\leq n$ we have

[TABLE]

Proof.

Note that $H(s)$ is the complete $U$-statistic estimator for $\epsilon_{k}(s)$. For each subset $S$ of size $s$ from $[n]$ we have

[TABLE]

Further, $\|x-x^{\prime}\|_{2}\leq\Delta_{\mathcal{X}}\leq\Delta$ holds for any $x^{\prime}\in\mathcal{X}$. Therefore, using Hoeffding’s inequality for $U$-statistics stated in Proposition 4, for any fixed $s$, w.p. $1-\delta$ we have

[TABLE]

Note that $\lfloor{z\rfloor}\geq z/2$ for $z\geq 1$ and therefore the above translates to

[TABLE]

Taking a union bound over $s=k,k+1,\ldots,n$, replacing $\delta=\delta/n$, and using $p\geq 1$, implies the result. ∎

Lemma 14**.**

Consider the selection process mentioned in Section 5.3 and let $s_{1}$ be the output of this process. Then, w.p. $1-\delta$ we have

[TABLE]

Proof.

Note that using Lemma 13, w.p. $1-\delta$, for all values of $s$ we have $|H(s)-\epsilon_{k}(s)|\leq G_{\delta}(s)$. Now consider three different cases:

• •

$s_{1}\geq s_{2}\geq s^{*}:$ Note that based on the choice of $s_{1},s_{2}$, we have $H(s_{2})>2G_{\delta}(s_{2})$. However, $H(s_{2})\leq\epsilon_{k}(s_{2})+G_{\delta}(s_{2})$. Hence, $\epsilon_{k}(s_{2})>G_{\delta}(s_{2})$ which contradicts with the assumption that $s_{2}\geq s^{*}$. Note that this is true since $\epsilon_{k}(s)-G_{\delta}(s)$ is non-positive for $s\geq s^{*}$.

• •

$s_{1}=s^{*},s_{2}=s^{*}-1:$ Obviously $s_{1}\leq s^{*}$.

• •

$s_{2}\leq s_{1}\leq s^{*}-1:$ Note that we have

[TABLE]

Hence, $G_{\delta}(s^{*}-1)<\epsilon_{k}(s^{*}-1)\leq\epsilon_{k}(s_{1})\leq 3G_{\delta}(s_{1})$. This means that $G_{\delta}(s^{*}-1)/G_{\delta}(s_{1})\leq 3$ which implies $\sqrt{(s^{*}-1)/s_{1}}\leq 3$. Therefore, $s_{1}\geq(s^{*}-1)/9$.

This completes the proof. ∎

Now we are ready to finalize the proof of Theorem 1. Note that using the result of Lemma 14, w.p. $1-\delta$, we have

[TABLE]

This basically means that if $s_{*}=9s_{1}+1$, then $s_{*}$ belongs to $[s^{*},10s^{*}]$. Hence, we have $\epsilon_{k}(s_{*})\leq\epsilon_{k}(s^{*})\leq G_{\delta}(s^{*})$ and $G_{\delta}(s_{*})\leq G_{\delta}(10s^{*})=\sqrt{10}G_{\delta}(s^{*})$. Now using Theorem 1, for $B\geq n/s_{*}$ w.p. $1-\delta$ we have

[TABLE]

Note that $G_{\delta}(s_{*})=\Delta\sqrt{\frac{2ps_{*}}{n}\log(2pn/\delta)}$. Therefore,

[TABLE]

Replacing this in above equation together with a union bound implies that w.p. at least $1-2\delta$ we have

[TABLE]

which finishes the first part of the proof. For the second part, note that according to Corollary 2, for the $k$-NN kernel $\epsilon(s)\leq Cs^{-1/d}$, for a constant $C$. Note that at $s=s^{*}-1$ we have

[TABLE]

for a constant $C$. The above implies that

[TABLE]

Hence,

[TABLE]

Remark 3**.**

Note that although computation of $H(s)$ may look complex as it involves calculation of distance to $k$-nearest neighbor of $x$ on all $\binom{n}{s}$ subsets, there is a closed form representation for $H(s)$ according to its representation based on $L$-statistic. In fact, by sorting samples $(X_{1},X_{2},\ldots,X_{n})$ based on their distance to $x$, i.e, $\|x-X_{(1)}\|_{2}\leq\|x-X_{(2)}\|_{2}\leq\ldots\leq\|x-X_{(n)}\|_{2}$, we have

[TABLE]

Therefore, after sorting training samples, we can compute values of $H(s)$ very efficient and fast.

H.7 Proof of Proposition 2

Note that according to Lemma 14, w.p. $1-1/n$, the output of process, $s_{1}$ satisfies

[TABLE]

where $s^{*}$ is the point for which we have $\epsilon_{k}(s^{*})=G_{1/n}(s^{*})$. This basically means that $s_{*}=9s_{1}+1\geq s^{*}$. Note that for the $k$-NN kernel we have $\eta_{k}(s)=\Theta(1/s)$. As $s_{\zeta}\geq n^{\zeta}$, this also implies that $\epsilon_{k}(s_{\zeta},\eta_{k}(s_{\zeta}))=O((\log(s_{\zeta})/s_{\zeta})^{1/d})\rightarrow 0$. Also, according to the inequality $\zeta<\frac{\log(n)-\log(s_{*})-\log\log^{2}(n)}{\log(n)}$ we have $1-\zeta>(\log(s_{*})+\log\log^{2}(n))/\log(n)$ and therefore

[TABLE]

and hence $n\eta_{k}(s_{\zeta})\rightarrow 0$. Finally, note that $\sigma_{n,j}(x)=\Theta(\sqrt{s/n})$ and according to Theorem 2 it suffices that

[TABLE]

Note that for any $\zeta>0,s_{\zeta}\geq s^{*}$ and therefore $\epsilon_{k}(s_{\zeta})\leq\epsilon_{k}(s^{*})=G_{1/n}(s^{*})$. For the first term,

[TABLE]

Now note that $s_{\zeta}=s_{*}n^{\zeta}\geq s^{*}n^{\zeta}$ and hence $\sqrt{s^{*}/s_{\zeta}}\log(n)=O(n^{-\zeta/2}\log(n))\rightarrow 0$. For the second term, note that again $s_{\zeta}\geq s^{*}$ and therefore $\epsilon_{k}(s_{\zeta})\leq\epsilon_{k}(s^{*})=G_{1/n}(s^{*})\leq G_{1/n}(s_{\zeta})$. Now note that since $s_{\zeta}/n\leq 1/\log^{2}(n)$ hence

[TABLE]

Finally, for the last term we have $s_{\zeta}/n\leq 1/\log^{2}(n)$ and hence

[TABLE]

This basically means w.p. $1-1/n$, $s_{\zeta}$ belongs to the interval for which the asymptotic normality result in Theorem 2 holds. As $n\rightarrow\infty$, the conclusion follows.

Appendix I Stochastic Equicontinuity of $U$-statistics via Bracketing

We define here some standard terminology on bracketing numbers in empirical process theory. Consider an arbitrary function space $\mathcal{F}$ of functions from a data space $\mathcal{Z}$ to $\mathbb{R}$, equipped with some norm $\|\cdot\|$. A bracket $[a,b]\subseteq\mathcal{F}$, where $a,b:\mathcal{Z}\rightarrow\mathbb{R}$ consists of all functions $f\in\mathcal{F}$, such that $a\leq f\leq b$. An $\epsilon$-bracket is a bracket $[a,b]$ such that $\|a-b\|\leq\epsilon$. The bracketing number $N_{[]}(\epsilon,\mathcal{F},\|\cdot\|)$ is the minimum number of $\epsilon$-brackets needed to cover $\mathcal{F}$. The functions $[a,b]$ used in the definition of the brackets need not belong to $\mathcal{F}$ but satisfy the same norm constraints as functions in $\mathcal{F}$. Finally, for an arbitrary measure $P$ on $\mathcal{Z}$, let

[TABLE]

Lemma 15** (Stochastic Equicontinuity for $U$-statistics via Bracketing).**

Consider a function space ${\cal F}$ of symmetric functions from some data space $\mathcal{Z}^{s}$ to $\mathbb{R}$ and consider the $U$-statistic of order $s$, with kernel $f$ over $n$ samples:

[TABLE]

Suppose $\sup_{f\in\mathcal{F}}\|f\|_{P,2}\leq\eta$, $\sup_{f\in\mathcal{F}}\|f\|_{P,\infty}\leq G$ and let $\kappa=n/s$. Then for $\kappa\geq\frac{G^{2}}{\log N_{[]}(1/2,\mathcal{F},\|\cdot\|_{P,2})}$, w.p. $1-\delta$:

[TABLE]

Proof.

Let $\kappa=n/s$. Moreover, wlog we will assume that $\mathcal{F}$ contains the zero function, as we can always augment $\mathcal{F}$ with the zero function without changing the order of its bracketing number. For $q=1,\ldots,M$, let $\mathcal{F}_{q}=\cup_{i=1}^{N_{q}}\mathcal{F}_{qi}$ be a partition of $\mathcal{F}$ into brackets of diameter at most $\epsilon_{q}=2\eta/2^{q}$, with $\mathcal{F}_{0}$ containing a single partition of all the functions. Moreover, we assume that $\mathcal{F}_{q}$ are nested partitions. We can achieve the latter as follows: i) consider a minimal bracketing cover of $\mathcal{F}$ of diameter $\epsilon_{q}$, ii) assign each $f\in\mathcal{F}$ to one of the brackets that it is contained arbitrarily and define the partition $\bar{\mathcal{F}}_{q}$ of size $\bar{N}_{q}=N_{[]}(\epsilon_{q},\mathcal{F},\|\cdot\|_{P,2})$, by taking $\mathcal{F}_{qi}$ to be the functions assigned to bracket $i$, iii) let $\mathcal{F}_{q}$ be the common refinement of all partitions $\bar{\mathcal{F}}_{0},\ldots,\bar{\mathcal{F}}_{q}$. The latter will have size at most $N_{q}\leq\prod_{q=0}^{M}\bar{N}_{q}$. Moreover, assign a representative function $f_{qi}$ to each partition $\mathcal{F}_{qi}$, with the representative for the single partition at level $q=0$ is the zero function.

Chaining definitions.

Consider the following random variables, where the dependence on the random input $Z$ is hidden:

[TABLE]

for some sequence of numbers $\alpha_{0},\ldots,\alpha_{M}$, to be chosen later. By noting that $A_{q-1}f=A_{q}f+B_{q}f$ and continuously expanding terms by adding and subtracting finer approximations to $f$, we can write the telescoping sum:

[TABLE]

For simplicity let $\mathbb{P}_{s,n}f=\Psi(f,Z_{1:n})$, $\mathbb{P}f=E[f(Z_{1:s})]$ and $\mathbb{G}_{s,n}$ denote the $U$-process:

[TABLE]

Our goal is to bound $\|\mathbb{P}_{s,n}f\|_{\mathcal{F}}=\sup_{f\in\mathcal{F}}|\mathbb{P}_{s,n}f|$, with high probability. Observe that since $\mathcal{F}_{0}$ contains only the zero function, then $\mathbb{G}_{s,n}f_{0}=0$. Moreover, the operator $\mathbb{G}_{s,n}$ is linear. Thus:

[TABLE]

Moreover, by triangle inequality:

[TABLE]

We will bound each term in each summand separately.

Edge cases.

The final term we will simply bound it by $2\alpha_{M}$, since $|(f-\pi_{M}f)A_{M}f|\leq\alpha_{M}$, almost surely. Moreover, the summand in the first term for $q=0$, we bound as follows. Observe that $B_{0}f=1\{\sup_{f}|f|>\alpha_{0}\}$. But we know that $\sup_{f}|f|\leq G$, hence: $B_{0}f\leq 1\{G>\alpha_{0}\}$.

[TABLE]

Hence, if we assume that $\alpha_{0}$ is large enough such that $\alpha_{0}>G$, then the latter term is zero. By the setting of $\alpha_{0}$ that we will describe at the end, the latter would be satisfied if $\kappa\geq\frac{G^{2}}{\log N_{[]}(1/2,\mathcal{F},\|\cdot\|_{P,2})}$.

$B_{q}$ terms.

For the terms in the first summand we have by triangle inequality:

[TABLE]

Moreover, observe that:

[TABLE]

where we used the fact that the partitions in $\mathcal{F}_{q}$, have diameter at most $\epsilon_{q}$, with respect to the $\|\cdot\|_{P,2}$ norm. Now observe that because the partitions $\mathcal{F}_{q}$ are nested, $\Delta_{q}f\leq\Delta_{q-1}f$. Therefore, $\Delta_{q}fB_{q}f\leq\Delta_{q-1}fB_{q}f\leq\alpha_{q-1}$, almost surely. Moreover, $\|\Delta_{q}fB_{q}f\|_{P,2}\leq\|\Delta_{q}f\|_{P,2}\leq\epsilon_{q}$. By Bernstein’s inequality for $U$ statistics (see e.g. Peel et al. [2010]) for any fixed $f$, w.p. $1-\delta$:

[TABLE]

Taking a union bound over the $N_{q}$ members of the partition, and combining with the bound on $\mathbb{P}\Delta_{q}fB_{q}f$, we have w.p. $1-\delta$:

[TABLE]

$A_{q}$ terms.

For the terms in the second summand, we have that since the partitions are nested, $|(\pi_{q}f-\pi_{q-1}f)A_{q-1}f|\leq\Delta_{q-1}fA_{q-1}f\leq\alpha_{q-1}$. Moreover, $\|(\pi_{q}f-\pi_{q-1}f)A_{q-1}f\|_{P,2}\leq\|\Delta_{q-1}f\|_{P,2}\leq\epsilon_{q-1}\leq 2\epsilon_{q}$. Thus, by similar application of Bernstein’s inequality for $U$-statistics, we have for a fixed $f$, w.p. $1-\delta$:

[TABLE]

As $f$ ranges there are at most $N_{q-1}\,N_{q}\leq N_{q}^{2}$ different functions $(\pi_{q}f-\pi_{q-1}f)A_{q-1}f$. Thus taking a union bound, we have that w.p. $1-\delta$:

[TABLE]

Taking also a union bound over the $2M$ summands and combining all the above inequalities, we have that w.p. $1-\delta$:

[TABLE]

Choosing $\alpha_{q}=\epsilon_{q}\sqrt{\kappa}/\sqrt{\log(2N_{q+1}M/\delta)}$ for $q<M$ and $\alpha_{M}=\epsilon_{M}$, we have for some constant $C$:

[TABLE]

Moreover, since $\log(N_{q})\leq\sum_{t=0}^{q}\log(N_{[]}(\epsilon_{q},\mathcal{F},\|\cdot\|_{P,2})$, we have:

[TABLE]

Combining all the above yields the result. ∎

Corollary 4**.**

Consider a function space ${\cal F}$ of symmetric functions. Suppose that $\sup_{f\in\mathcal{F}}\|f\|_{P,2}\leq\eta$ and $\log(N_{[]}(\epsilon,\mathcal{F},\|\cdot\|_{P,2})=O(1/\epsilon)$. Then for $\kappa\geq O(G^{2})$, w.p. $1-\delta$:

[TABLE]

Proof.

Applying Lemma 15, we get for every $\rho>0$, the desired quantity is upper bounded by:

[TABLE]

Choosing $\rho=\sqrt{\eta}/\sqrt{\kappa}$, yields the desired bound. ∎

Appendix J Proof of Lemma 8

We will argue asymptotic normality of the $U$-statistic defined as:

[TABLE]

under the assumption that for any subset of indices $S_{b}$ of size $s$: $\mathbf{E}\left[\mathbf{E}[\alpha_{S_{b},\omega_{b}}(X_{1})|X_{1}]^{2}\right]=\eta(s)$ and that the kernel satisfies shrinkage in probability with rate $\epsilon(s,\delta)$ such that $\epsilon(s,\eta(s)^{2})\rightarrow 0$ and $n\eta(s)\rightarrow\infty$. For simplicity of notation we let:

[TABLE]

and we then denote:

[TABLE]

Observe that we can then re-write our $U$-statistic as:

[TABLE]

Moreover, observe that by the definition of $Y_{i}$, $\mathbf{E}[Y_{i}\mid X_{i}]=0$ and also

[TABLE]

Invoking Lemma 17, it suffices to show that: $\operatorname{Var}\left[\Phi_{1}(Z_{1})\right]=\Omega(\eta(s))$, where $\Phi_{1}(z_{1})=\mathbf{E}[\Phi(z_{1},Z_{2},\ldots,Z_{s})]$. The following lemma shows that under our conditions on the kernel, the latter property holds.

Lemma 16**.**

Suppose that the kernel $K$ is symmetric (Assumption 4), has been built in an honest manner (Assumption 3) and satisfies:

[TABLE]

Then, the following holds

[TABLE]

Proof.

Note we can write

[TABLE]

Here, $B$ is zero mean conditional on $X_{1}$ and also $A$ and $B$ are uncorrelated, i.e., $\mathbf{E}[AB]=\mathbf{E}[A]\mathbf{E}[B]=0$. Therefore:

[TABLE]

For simplicity of notation let $W_{i}=K(x,X_{1},\{Z_{j}\}_{j=1}^{s})$ denote the random variable which corresponds to the weight of sample $i$. Note that thanks to the honesty of kernel defined in Assumption 3, $W_{i}$ is independent of $Y_{1}$ conditional on $X_{1}$, for $i\geq 2$. Hence all the corresponding terms in the summation are zero. Therefore, the expression inside the variance above simplifies to

[TABLE]

Moreover, by honesty $W_{1}$ is independent of $Y_{1}$ conditional on $X_{1}$. Thus, the above further simplifies to:

[TABLE]

Using $\operatorname{Var}(G)=\mathbf{E}[G^{2}]-\mathbf{E}[G]^{2}$, this can be further rewritten as

[TABLE]

Note that $Y_{1}-\mathbf{E}[Y_{1}\mid X_{1}]$ is uniformly upper bounded by some $\psi_{\max}$. Furthermore, by the symmetry of the kernel we have $\mathbf{E}\left[\mathbf{E}[W_{1}\mid X_{1}]\right]=\mathbf{E}[W_{1}]=1/s$.888Since $\mathbf{E}[W_{i}]$ are all equal to the same value $\kappa$ and $\sum_{i}\mathbf{E}[W_{i}]=1$, we get $\kappa=1/s$. Thus the second term in the latter is of order $1/s^{2}$. Hence:

[TABLE]

Focusing at the first term and letting $\sigma^{2}(x)=\operatorname{Var}(Y|X=x)$, we have:

[TABLE]

The goal is to prove that the second term is $o(1/s)$. For ease of notation let $V_{1}=\mathbf{E}\left[W_{1}\mid X_{1}\right]$. Then we can bound the second term as:

[TABLE]

where we used the fact that $V_{1}\leq 1$, the assumption that $\sigma^{2}(\cdot)$ is $L_{mm}$-Lipschitz, the tower rule and the definition of $\eta(s)$. Furthermore,

[TABLE]

which by definition is at most $\delta$. By putting $\delta=\eta(s)^{2}$ we obtain

[TABLE]

where we invoked our assumption that $\epsilon(s,\eta(s)^{2})\rightarrow 0$. Thus we have obtained that:

[TABLE]

which is exactly the form of the lower bound claimed in the statement of the lemma. This concludes the proof. ∎

J.1 Hájek Projection Lemma for Infinite Order $U$-statistics

The following is a small adaptation of Theorem 2 of Fan et al. [2018], which we present here for completeness.

Lemma 17** (Fan et al. [2018]).**

Consider a $U$-statistic defined via a symmetric kernel $\Phi$:

[TABLE]

where $Z_{i}$ are i.i.d. random vectors and $s$ can be a function of $n$. Let $\Phi_{1}(z_{1})=\mathbf{E}[\Phi(z_{1},Z_{2},\ldots,Z_{s})]$ and $\eta_{1}(s)=\operatorname{Var}_{z_{1}}\left[\Phi_{1}(z_{1})\right]$. Suppose that $\operatorname{Var}{\Phi}$ is bounded, $n\,\eta_{1}(s)\rightarrow\infty$. Then:

[TABLE]

where $\sigma_{n}^{2}=\frac{s^{2}}{n}\eta_{1}(s)$.

Proof.

The proof follows identical steps as the one in Fan et al. [2018]. We argue about the asymptotic normality of a $U$-statistic:

[TABLE]

Consider the following projection functions:

[TABLE]

where $\mathbf{E}\left[\Phi\right]=\mathbf{E}\left[\Phi(Z_{1},\ldots,Z_{s})\right]$. Then we define the canonical terms of Hoeffding’s $U$-statistic decomposition as:

[TABLE]

Subsequently the kernel of the $U$-statistic can be re-written as a function of the canonical terms:

[TABLE]

Moreover, observe that all the canonical terms in the latter expression are un-correlated. Hence, we have:

[TABLE]

We can now re-write the $U$ statistic also as a function of canonical terms:

[TABLE]

Now we define the Hájek projection to be the leading term in the latter decomposition:

[TABLE]

The variance of the Hajek projection is:

[TABLE]

The Hájek projection is the sum of independent and identically distributed terms and hence by the Lindeberg-Levy Central Limit Theorem (see e.g., Billingsley [2008], Borovkov [2013]):

[TABLE]

We now argue that the remainder term: $\frac{U-\mathbf{E}\left[U\right]-\hat{U}}{\sigma_{n}}$ vanishes to zero in probability. The latter then implies that $\frac{U-\mathbf{E}\left[U\right]}{\sigma_{n}}\rightarrow_{d}\mathsf{N}(0,1)$ as desired. We will show the sufficient condition of convergence in mean square: $\frac{\mathbf{E}\left[\left(U-\mathbf{E}\left[U\right]-\hat{U}\right)^{2}\right]}{\sigma_{n}^{2}}\rightarrow 0$. From an inequality due to Wager and Athey [2018]:

[TABLE]

Since $\operatorname{Var}\left[\Phi(Z_{1},\ldots,Z_{n})\right]$ is bounded by a constant $V^{*}$ and $n\,\eta_{1}(s)\rightarrow\infty$, by our assumption, we have:

[TABLE]

∎