Fast rate of a nearest neighbor estimate for a global parameter

Author: Yen-Chi Chen (University of Washington)
Date: 07/08/2026

Local versus global parameter

We consider a simple nonparametric regression model where our data consists of IID random variables
$$(X_1,Y_1),\cdots, (X_n,Y_n)$$
from some unknown distributions such that $X_i\in\mathbb{R}^d$.
We further assume the following data generating process:
$$Y_i = m(X_i) + \epsilon_i,$$
where $X_i\perp \epsilon_i$ and $\epsilon_1,\cdots, \epsilon_n$ are IID mean $0$ noises with variance ${\sf Var}(\epsilon_i) = \sigma^2<\infty$.

It is well-knonwn that we cannot estimate the local parameter $m(x)$ at a very fast rate unless it is super smooth. Roughly speaking, if $m(x)$ has bounded $s$-times derivatives, the fastest rate that we can achieve is
$$\hat m(x) - m(x) = O_P\left(n^{-s/(2s+d)}\right).$$
The kNN estimator achieves this rate with a good choice of $k \sim n^{-s/(2s+d)}$.
Note that conventional choice is $s=2$ since we often assume that the function $m(x)$ is twice differentiable.

Now we consider estimating a global parameter
$$\theta = \mathbb{E}(m^2(X)).$$
This parameter is often of interest due to its relation to the variance $\sigma^2 = \mathbb{E}((Y-m(X))^2)$.

How to estimate $\theta$ in a nonparametric setting?

1-NN estimator

It turns out that a simple and accurate estimator is the 1-NN estimator
$$\hat \theta_n = \frac{1}{n} \sum_{i=1}^n Y_i Y_{N(i)},$$
where
$$N(i) = {\sf argmin}_{j\neq i} \|X_j-X_i\|$$
is the index of the $i$-th observation.

The above rate has two components:

• Variance. The $\frac{1}{\sqrt{n}}$ is from the variance of the estimator, which is the conventional rate.
• Bias. The second component $n^{-1/d}$ is the bias. If you are familiar with the kNN theory, you may notice that this is exactly rate of the distance between $\|X_i - X_{N(i)}\|$. In fact, if we use $k$-NN rather than 1-NN, it remains at the same rate as long as $k$ does not scale with respect to the sample size $n$.

Convergenace rate

Here is how we obtain the above rate. First, observed that
$$Y_i = m(X_i) + \epsilon_i,\qquad Y_{N(i)} = m(X_{N(i)}) + \epsilon_{N(i)}$$
so as long as $m$ is Lipschitz continuous, we have
$$m(X_{N(i)}) = m(X_i) + \delta_i,\qquad \delta_i = O_P(\|X_i - X_{N(i)}\|) = O_P(n^{-1/d}).$$
Putting this back into the estimator $\hat \theta_n$, we obtain
$$\begin{align*} \hat \theta_n &= \frac{1}{n}\sum_{i=1}^nY_i Y_{N(i)} \\ &= \frac{1}{n}\sum_{i=1}^n(m(X_i) +\epsilon_i) (m(X_i) + \delta_i + \epsilon_{N(i)}) \\ &= \frac{1}{n}\sum_{i=1}^n \underbrace{m^2(X_i)}_\text{signal} + \underbrace{m(X_i) \delta_i}_\text{bias} + \underbrace{\epsilon_{N(i)} Y_i + \epsilon_i Y_{N(i)}}_\text{noise}. \end{align*}$$

Clearly, the signal component $\frac{1}{n}\sum_{i=1}^n m^2(X_i)$ results in a $\sqrt{n}$-consistent estimator of $\theta$. So we just need to control the remaining quantities.

Bias component

Since it is well-known that $\delta_i = O_P(n^{-1/d})$, we immediately have
$$\frac{1}{n}\sum_{i=1}^n m(X_i) \delta_i =O_P(n^{-1/d})$$
under very mild conditions of $m(x)$ and the underly density function of the covariate $p_X$.

Noise component (part of the variance)

The remaining quantity
$$\bar V_n = \frac{1}{n}\sum_{i=1}^n \underbrace{\epsilon_{N(i)} Y_i +\epsilon_i Y_{N(i)}}_{=V_i}$$
behaves like average of noises and under proper moment conditions of $\epsilon_i$, we have
$${\sf Var} (V_i) = \sigma^2_V<\infty.$$

However, this is NOT enough to argue that ${\sf Var}(\bar V_n) = O_P(n^{-1/2})$ because $V_1,\cdots, V_n$ are dependent due to the 1NN component.
Formally, we have
$${\sf Var}(\bar V_n) = \frac{1}{n^2}\sum_{i=1}^n{\sf Var}(V_i) + \frac{1}{n^2}\sum_{i\neq j=1}^n {\sf Cov}(V_i, V_j).$$
The first variance is controlled as we have mentioned. The tricky part is the second component the covariance.

Using the idea of kissing number and the fact that $V_i,V_j$ will be correlated only if $\{X_i, X_{N(i)}\} \cap \{X_j, X_{N(j)}\}\neq \emptyset$, we realize that the set $\{j: {\sf Cov}(V_i,V_j)\neq 0\}$ has at most $C_d \leq 2K_d$ elements. So most pairs ${\sf Cov}(V_i,V_j) = 0$, and we conclude that
$$\frac{1}{n^2}\sum_{i\neq j=1}^n {\sf Cov}(V_i, V_j) \leq \frac{ 2 n K_d \times \sigma^2_V}{n^2} = O(n^{-1}),$$
where $\sigma^2_V<\infty$ is just some constant. Therefore, we conclude that
$${\sf Var}(\bar V_n) = O(n^{-1})$$
and hence
$$\hat \theta_n - \theta = O_P\left(\frac{1}{\sqrt{n}} + n^{-1/d}\right).$$

Estimating higher-order moments

The same estimator can be applied to cases where we are interested in the average of an $L$-th local moment. Specifically, consider $\theta_L = \mathbb{E}(m^L(X))$.
We can use a $(L-1)$-NN estimator:
$$\hat \theta_L = \frac{1}{n}\sum_{i=1}^n Y_i \prod_{\ell=1}^{L-1} Y_{N_{L-1}(i)},$$
where $N_{L-1}(i)$ is the indices of the $(L-1)$-NN of $X_i$.

Using the same derivation, the convergence rate of $\hat \theta_L$ remains the same:
$$\hat \theta_L - \theta_L = O_P\left(\frac{1}{\sqrt{n}} + n^{-1/d}\right).$$

Another example of a slightly faster rate: $\eta = \mathbb{E}(e^{m(X)})$

Now we consider the problem of estimating $\eta = \mathbb{E}(e^{m(X)})$. Since the exponential is a inifinite-order polynomial, the above fixed $k$ approach is not applicable so we have to consider $k=k_n\rightarrow \infty$ as $n\rightarrow\infty$.

Consider a simple plug-in estimator:
$$\hat \eta_k = \frac{1}{n}\sum_{i=1}^n e^{\hat m_k(X_i)}.$$

How does the rate improvement occur?

To see how the rate improve occurs, consider the orcale estimator
$$\tilde \eta = \frac{1}{n}\sum_{i=1}^n e^{m(X_i)}.$$
Clearly, $\tilde \eta - \eta = O_P\left(\frac{1}{\sqrt{n}}\right)$.
But this estimator and the kNN estimator are similar via a Taylor expansion:
$$e^{\hat m_k(X_i)} - e^{m(X_i)} \approx e^{m(X_i)}(\hat m_k(X_i) - m(X_i)) + \frac{1}{2} e^{m(X_i)}(\hat m_k(X_i) - m(X_i))^2.$$

Therefore,
$$\begin{align*} \hat \eta - \tilde \eta \approx \underbrace{\frac{1}{n}\sum_{i=1}^n e^{m(X_i)}(\hat m_k(X_i) - m(X_i))}_\text{(linear)} + \underbrace{\frac{1}{n}\sum_{i=1}^n\frac{1}{2} e^{m(X_i)}(\hat m_k(X_i) - m(X_i))^2}_\text{(MSE)}. \end{align*}$$

• (linear): key improvement of the rate (a.k.a. variance collapse). This part is the key improvement of the convergence rate. Roughly speaking, the estimator
  $$\hat m_k(X_i) = \frac{1}{k}\sum_{j \in N_k(i)} Y_j = \frac{1}{k}\sum_{j\in N_k(i)} m(X_j) +\epsilon_j.$$

  • So the $i$-th term in (linear) has a bias at rate
    $$B_i = \frac{1}{k}\sum_{j \in N_k(i)} e^{m(X_i)}|m(X_j) - m(X_i)| = O_P\left(\left(\frac{k}{n}\right)^{s/d}\right)$$
    and the variance $\frac{1}{k}\sum_{j\in N_k(i)}e^{m(X_i)} \epsilon_j$.
  • The average does not change the bias rate and the variance part has an interesting phenomenon:
    $$\frac{1}{n}\sum_{i=1}^n \sum_{j\in N_k(i)}e^{m(X_i)} \epsilon_j = \frac{1}{n}\sum_{i=1}^n W_i \epsilon_j,$$
    where $W_i$ is related to tthe number of times $X_i$ is within $k$-NN of other observation $X_j$ and the corresponding $e^{m(X_j)}$. Again, due to the geometry (similar to the kissing number), each point can only appear up to certain times in others' kNN, so this weight will be upper bounded.
  • Thus, the average $\frac{1}{n}\sum_{i=1}^n W_i \epsilon_j$ leads to a rate $O_P\left(\frac{1}{\sqrt{n}}\right)$. Putting these together, we conclude that
    $$\text{(linear)} = O_P\left(\frac{1}{\sqrt{n}}+\left(\frac{k}{n}\right)^{s/d}\right).$$

• (MSE): the conventional rate. The second term is like a sample version of the mean square error (MSE). So it is at the conventional rate of the MSE
  $$\text{(MSE)} = O_P\left(\frac{1}{k}+\left(\frac{k}{n}\right)^{2s/d}\right)$$

Therefore, the dominating error is the bias from (linear) and the variance from (MSE), which corresponds to the rate
$$\hat \eta_k - \eta = \hat \eta_k - \tilde \eta + \tilde \eta - \eta = O_P\left(\frac{1}{\sqrt{n}}+\frac{1}{k}+\left(\frac{k}{n}\right)^{s/d}\right)$$
and under the optimal choice of $k\sim n^{s/(s+d)}$, we obtain the rate
$$\hat \eta_k - \eta = O_P\left(n^{-s/(s+d)}\right).$$

An example of no improvements: $\rho = \mathbb{E}(I(m(X)>0.5))$

When the parameter of interest is not a polynomial functional, this kNN trick fails and we have to retreat to the conventional local estimator, although we may still have a fast rate with some other conditions.

Consider the quantity
$$\rho = \mathbb{E}(I(m(X) > 0.5)).$$
This is a common quantity in binary classification, where it measures the predictive probability of class $1$ if we convert the regression model into a classifier.

To see this, if you pick any fixed $k$, say $k=5$, the kNN estimator $\hat m_k$ has a non-vanishing variance. So for any $x$ with $m(x) \approx 0.5$, this non-vanishing variance could push it up or down the critical threshold $0.5$, leading to a non-vanishing error regardless of the sample size. Therefore, we have to set $k$ to slowly increase with respect to the sample size to achieve consistency. Moreover, the Taylor expansion method in the previous section is not applicable since the indicator function is not differentiable.

Improving rate by Tsybakov's margin condition

Note that in this scenario, we may improve the convergence rate via the Tsybakov's margin condition:

Under Tsybakov's margin condition, the convergence rate of an $k$-NN estimator of $\rho$ is
$$\hat \rho_k - \rho = O_P\left(n^{-\frac{s \alpha}{2s+d}}\right)$$
when we choose $k \sim n^{s/(2s+d)}$.

Rate comparison for estimating a global parameter

Now we see three regimes of the convergence rate. To put everything in the same condition, we consider $m$ to be only Lipschitz smooth. A sufficient condition is $m$ to have a uniformly bounded derivative. So we can view this as the case of $s=1$.
Recall that the convergence rate for estimating a local parameter $m(x)$ is
$$\hat m_k(x) - m(x) = O_P\left(n^{-1/(2+d)}\right)$$
when we choose $k\sim n^{1/(2+d)}$.

In this case, we can see three convergence rates:

• Ultra fast regime: fixed $k$. For estimating global parameter that is the expectation of a polynomial of $m(X)$, we have rate
  $$\hat \theta_n - \theta = O_P\left(n^{-1/2} + n^{-1/d}\right).$$
  In this case, we achieve the parametric rate when $d\leq 2$.
• Fast regime: quickly diverging $k$. For estimating a global parameter that is the expecatation of a smooth transformation of $m(X)$, we have a rate
  $$\hat \eta_k - \theta = O_P\left(n^{-1/2} + n^{-1/(1+d)}\right)$$
  when we choose $k \sim n^{1/(1+d)}$. We can achieve parametric rate when $d=1$.
• Regular regime. If the global parameter involves a non-smooth transformation of $m(X)$, we may not be able to obtain any rate improvement relative to the local parameter:
  $$\hat \rho_k - \rho =O_P\left(n^{-1/2} + n^{-1/(2+d)}\right)$$
  when we choose $k\sim n^{1/(2+d)}$. We cannot acheive parametric rate in this scenario.