The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
30,101 characters
Debiased Nonparametric Regression for Statistical Inference and Distributionally Robustness
\maketitle
\begin{abstract}
This study proposes a debiasing method for smooth nonparametric estimators. While machine learning techniques such as random forests and neural networks have demonstrated strong predictive performance, their theoretical properties remain relatively underexplored. In particular, many modern algorithms lack guarantees of pointwise and uniform risk convergence, as well as asymptotic normality. These properties are essential for statistical inference and robust estimation and have been well-established for classical methods such as Nadaraya-Watson regression. To ensure these properties for various nonparametric regression estimators, we introduce a model-free debiasing method. By incorporating a correction term that estimates the conditional expected residual of the original estimator, or equivalently, its estimation error, into the initial nonparametric regression estimator, we obtain a debiased estimator that satisfies pointwise and uniform risk convergence, along with asymptotic normality, under mild smoothness conditions. These properties facilitate statistical inference and enhance robustness to covariate shift, making the method broadly applicable to a wide range of nonparametric regression problems.
\end{abstract}
\clearpage
\section{Introduction}
This study investigates the problem of nonparametric regression. While modern machine learning methods enable the estimation of complex regression functions, they often lack the theoretical guarantees that traditional estimators possess, such as pointwise and uniform mean squared error (MSE) convergence, as well as asymptotic normality. Asymptotic normality is crucial for statistical inference, whereas uniform risk convergence plays an important role in prediction under covariate shift. In this study, we propose a debiasing method for regression estimators and show that, under mild smoothness conditions, they achieve asymptotic normality and uniform convergence.
Here, we formulate the problem. Let $X \in {\mathcal{X}}$ be a $d$-dimensional covariate and $Y \in {\mathcal{Y}}$ be a target variable, where ${\mathcal{X}} \subseteq {\mathbb{R}}^d$ and ${\mathcal{Y}} \subset {\mathbb{R}}$ represent the covariate and target spaces, respectively. Let $P$ denote the joint distribution of $(X, Y)$, and define the regression function under $P$ as
\[
f_P(X) = {\mathbb{E}}_P\big[Y \mid X\big],
\]
where ${\mathbb{E}}_P\big[\cdot\big]$ represents the expectation operator with respect to $P$. We assume access to observations $\big\{(X_i, Y_i)\big\}_{i=1}^n$, where $(X_i, Y_i)$ are independent and identically distributed samples from a true distribution $P_0$. The true regression function is denoted by $f_0 = f_{P_0}$. Our objective is to estimate $f_0$ from the observations $\big\{(X_i, Y_i)\big\}_{i=1}^n$, aiming for pointwise and uniform risk convergence, as well as asymptotic normality.
\paragraph{Notation.} For a real-valued vector $z$, let $\|z\|$ denote its Euclidean norm. For a measurable function $h\colon {\mathcal{X}} \to {\mathbb{R}}$, let $\|h\|_\infty \coloneqq \sup_{x \in {\mathcal{X}}} |h(x)|$ denote the sup-norm, and $\|h\|_2 \coloneqq \sqrt{{\mathbb{E}}\left[h(X)^2\right]}$ denote the $L^2$ norm.
\subsection{Content of this study}
In Section~\ref{sec:debiased_estimator}, we define our debiased estimator for nonparametric regression. Our proposed estimator consists of the following three steps: (1) estimate the regression function $f_0$ using a nonparametric regression method; (2) estimate the expected conditional residual, or equivalently, the estimation error of the first-stage estimator, using local polynomial regression; and (3) add the second-stage estimator to the first-stage estimator.
In Section~\ref{sec:local_polynomial}, we present an example of our proposed debiased estimator, applying local polynomial regression for conditional expected residual estimation.
In Section~\ref{sec:conv_analysis}, we establish that our debiased nonparametric regression estimator achieves pointwise and uniform MSE convergence if either (a) the first-stage nonparametric regression estimator and the true regression function, or (b) their difference, belongs to the H\"older class. Using pointwise MSE convergence, we also prove asymptotic normality. Furthermore, our estimator possesses the doubly robust property: if either the first-stage nonparametric regression estimator or the second-stage conditional expected residual estimator is consistent, then the resulting debiased estimator is also consistent.
The key takeaway is that our debiased estimator ensures desirable theoretical properties, typically associated with classical methods, for a wide range of modern first-stage estimators. As long as certain smoothness conditions are met, our method effectively corrects bias in a broad class of initial estimators, equipping the resulting estimator with strong theoretical guarantees.
\subsection{Related work}
Doubly robust and debiased estimation has attracted significant attention across various fields, including statistics \citep{tsiatis2007semiparametric}, economics \citep{ChernozhukovVictor2018Dmlf}, epidemiology \citep{BangRobins2005}, and machine learning \citep{Kallus2020double}. With the advancement of machine learning regression models, debiasing methods based on doubly robust estimators have been extensively studied \citep{ChernozhukovVictor2018Dmlf}, particularly in conjunction with sample splitting techniques \citep{klaassen1987,vanderVaart2002,ZhengWenjing2011CTME}. While much of the existing literature focuses on parametric or semiparametric models, this study extends these approaches to nonparametric regression. The application of debiased machine learning to nonparametric regression is discussed in \citet{kennedy2023semiparametricdoublyrobusttargeted}, though their focus differs from ours.
Nonparametric regression in machine learning presents several challenges:
\begin{itemize}
\item Pointwise and uniform MSE convergence is not established for many estimators, such as neural networks and random forests.
\item Asymptotic normality remains unproven for numerous machine learning estimators.
\end{itemize}
Pointwise and uniform MSE convergence and asymptotic normality are critical in statistical analysis, motivating significant efforts to establish these properties. For instance, \citet{Wager2018} and \citet{Mourtada2020} demonstrate asymptotic normality and optimality, respectively, for their modified random forests, but not for the original random forest of \citet{Breiman2001}.
Pointwise and uniform convergence is particularly relevant for robustness under distributional shifts, such as covariate shift \citep{Shimodaira2000}. These convergence properties are rarely established for modern nonparametric estimators, whereas classical methods like local linear and series regression possess them. This gap arises in part due to the data-adaptive nature of modern methods and the frequent use of empirical process arguments, which often focus on population MSE, where expectations are taken over the covariate distribution rather than on pointwise or sup-norm MSE. \citet{SchmidtHieber2024} addresses this concern by demonstrating restricted uniform optimality for neural network regression in the one-dimensional covariate setting, but this approach does not generalize easily.
Our method shares some motivations with debiasing approaches in high-dimensional regression. For example, the Lasso estimator, while widely used, introduces a bias that decreases with sample size at a slower rate than $\sqrt{n}$, thereby hindering asymptotic normality \citep{tibshirani96regression,Buhlmann2011}. To mitigate this, \citet{vandeGeer2014} proposes incorporating bias-correction terms, leading to improved asymptotic properties such as normality and efficiency \citep{Jankova2018}.
This study generalizes our previous work on doubly robust methods for nonparametric regression discontinuity design \citep{kato2024doublyrobustregressiondiscontinuity}. Concurrently with our work, \citet{chernozhukov2024conditionalinfluencefunctions} and \citet{vanderlaan2025automaticdebiasedmachinelearning} also explore double machine learning in the context of nonparametric regression, though their emphases differ substantially from ours. For example, \citet{chernozhukov2024conditionalinfluencefunctions} develop conditional influence functions, which are related to the framework of \citet{Ichimura2022}.
\section{The debiased estimator}
\label{sec:debiased_estimator}
In this section, we define our debiased estimator in a general way. First, we randomly split the observations ${\mathcal{D}}$ into two datasets, ${\mathcal{D}}^{(1)}$ and ${\mathcal{D}}^{(2)}$, such that ${\mathcal{D}} = {\mathcal{D}}^{(1)} \cup {\mathcal{D}}^{(2)}$. For simplicity, let $n$ be even and set $m = n/2$. Then, for each $\ell \in \{1, 2\}$, we define
\[
{\mathcal{D}}^{(\ell)} \coloneqq \big\{(X_i, Y_i)\big\}_{i \in {\mathcal{I}}^{(\ell)}},
\]
where ${\mathcal{I}}^{(\ell)}$ is the index set of ${\mathcal{D}}^{(\ell)}$, i.e., ${\mathcal{I}}^{(1)} \cup {\mathcal{I}}^{(2)} = \{1, 2, \dots, n\}$.
We consider the following three-stage estimation for each point $x_0 \in {\mathcal{X}}$ of interest:
\begin{description}
\item[First-stage:] Estimate $f_0$ using \emph{any} smooth model, and denote this estimator by $\widehat{f}_n$.
\item[Second-stage:] Estimate the conditional expected residual of the first-stage estimator, ${\mathbb{E}}\big[Y - f(X)\mid X = x_0\big]$, or equivalently, the estimation error $f_0(x_0) - \widehat{f}_n(x_0)$. Denote this estimator by $\widehat{b}_n(x_0)$.
\item[Third-stage:] Sum the first- and second-stage estimators to obtain
\begin{align}
\label{eq:debiased}
&\widetilde{f}_n(x_0) \coloneqq \widehat{b}_n(x_0) + \widehat{f}_n(x_0),
\end{align}
which serves as a debiased estimator of $f_0$.
\end{description}
The term debiased estimator refers to our objective of correcting bias arising from the first-stage estimator. In many nonparametric and machine learning methods, bias can be problematic. By subtracting the estimated residual, ${\mathbb{E}}\big[Y - \widehat{f}_n(X) \mid X = x_0\big]$, we compensate for the difference between $f_0$ and $\widehat{f}_n$, resulting in an estimator with improved convergence and inferential properties.
Our estimator (\ref{eq:debiased}) is closely related to the influence function in conditional mean estimation \citep{kennedy2023semiparametricdoublyrobusttargeted}. It is also connected to the Neyman orthogonal score \citep{ChernozhukovVictor2018Dmlf}, which is essentially equivalent to the canonical gradient in one-step bias correction \citep{Schuler2024}. Notably, in conditional mean estimation, constructing the influence function is more intricate due to the conditioning. For an influence function that involves parameters represented via the conditional expected value, see \citet{Ichimura2022}. Independently and concurrently, \citet{chernozhukov2024conditionalinfluencefunctions} investigates a related topic with different motivations and methods.
\section{The debiased nonparametric regression with local polynomial conditional expected residual estimation}
\label{sec:local_polynomial}
In this section, we present an example of our debiased estimator. We employ local polynomial regression to estimate the conditional expected residual. For simplicity, we focus on the case where the covariate is one-dimensional ($d = 1$) and ${\mathcal{X}} = [0, 1]$.
\subsection{First-stage nonparametric regression}
For the first-stage estimator $\widehat{f}_n$, any regression method can be used, provided that smoothness conditions are met. As we show later, our key assumption is that the difference $f_0 - \widehat{f}_n$ is smooth. While various notions of smoothness exist, for simplicity, we focus on smoothness in the sense of the H\"older class (Definition~\ref{def:holder}). Notably, neither $f_0$ nor $\widehat{f}_n$ needs to be smooth individually, as long as their difference $f_0 - \widehat{f}_n$ is smooth.
We emphasize that our analysis does not impose any specific requirements on $\widehat{f}_n$ beyond smoothness. The first-stage estimator $\widehat{f}_n$ may converge to $f_0$ at an arbitrarily slow rate or even be inconsistent. Nevertheless, as long as the smoothness conditions hold, our asymptotic theoretical results remain valid. However, to achieve strong finite-sample performance, a high-quality first-stage estimator is typically preferred. In practice, methods such as random forests, neural networks, or other advanced machine learning techniques can be employed, provided they satisfy the smoothness conditions to a reasonable extent.
\subsection{Second-stage conditional expected residual estimation}
In the second-stage estimation of the conditional expected residual, we employ local polynomial regression. We define the $\ell$-th order polynomial basis as follows:
\[
\rho(u) \coloneqq \begin{pmatrix}
1 & u & u^2/2! & \cdots & u^\ell / \ell!
\end{pmatrix}^\top.
\]
In the theoretical analysis, the value of $\ell$ is chosen based on the smoothness of $f_0 - \widehat{f}_n$.
Using this polynomial basis, we define the local polynomial regression estimator as follows:
\begin{align*}
&\widehat{\beta}_n(x) \coloneqq \begin{pmatrix}
\widehat{\beta}_{0, n}(x) & \widehat{\beta}_{1, n}(x) & \widehat{\beta}_{2, n}(x) & \cdots & \widehat{\beta}_{\ell, n}(x)
\end{pmatrix}^\top\\
&\coloneqq \operatorname*{arg\,min}_{\beta \in {\mathbb{R}}^{\ell + 1}}\sum_{i\in{\mathcal{I}}^{(2)}}\left(Y_i - \widehat{f}_n(X_i) - \beta^\top \rho\left(\frac{X_i - x}{h_n}\right)\right)^2 K\left(\frac{X_i - x}{h_n}\right),
\end{align*}
where $K\colon {\mathcal{X}} \to {\mathbb{R}}$ is a kernel function defined by
\[
K(u) \coloneqq K_h(u) \coloneqq \mathbbm{1}\left[\big|u\big| \leq h\right].
\]
The conditional expected residual estimator is then given by the first component of $\widehat{\beta}_n(x)$:
\[
\widehat{b}_n(x) = \widehat{\beta}_{0, n}(x).
\]
\section{Convergence analysis}
\label{sec:conv_analysis}
This section presents a convergence analysis of $\widetilde{f}_n$ defined in Section~\ref{sec:local_polynomial}. We begin by defining the H\"older class.
\begin{definition}[H\"older class]
\label{def:holder}
Given an interval ${\mathcal{T}} \subset {\mathbb{R}}$, and positive constants $\beta$ and $L$, the H\"older class $\Sigma(\beta, L)$ on ${\mathcal{T}}$ is the set of $\ell = \floor{\beta}$ times differentiable functions $f \colon {\mathcal{T}} \to {\mathbb{R}}$ whose $\ell$-times derivative $f^{(\ell)}$ satisfies
\[
\Big|f^{(\ell)}(x) - f^{(\ell)}(z)\Big|\leq L \big|x - z\big|^{\beta - \ell}\quad \forall x, z \in {\mathcal{T}}.
\]
\end{definition}
\subsection{Closed-form solution}
The debiased estimator with local polynomial regression admits a closed-form solution. We define
\begin{align*}
\widehat{{\mathcal{B}}}_{n}(x_0)
&\coloneqq \frac{1}{mh_n}\sum_{i\in{\mathcal{I}}^{(2)}}
\rho\left(\frac{X_i - x_0}{h_n}\right)
\rho\left(\frac{X_i - x_0}{h_n}\right)^\top
K\left(\frac{X_i - x_0}{h_n}\right),\\
\widehat{w}_{h}(x, x_0)
&\coloneqq \frac{1}{h_n}\rho(0)^\top
\widehat{{\mathcal{B}}}_{n}(x_0)^{-1}
\rho\left(\frac{x - x_0}{h_n}\right)
K\left(\frac{x - x_0}{h_n}\right).
\end{align*}
Then, the second-stage estimator is given by
\[
\widehat{b}_n(x_0)
= \frac{1}{m}\sum_{i\in{\mathcal{I}}^{(2)}}
\left(Y_i - \widehat{f}_n(X_i)\right)
\widehat{w}_{h}(X_i, x_0).
\]
Thus, the debiased estimator is expressed as
\[
\widetilde{f}_n(x_0) = \frac{1}{m}\sum_{i\in{\mathcal{I}}^{(2)}}
\left(Y_i - \widehat{f}_n(X_i)\right)
\widehat{w}_{h}(X_i, x_0) + \widehat{f}_n(x_0).
\]
\subsection{Bias and variance decomposition}
We begin by analyzing the bias and variance of $\widetilde{f}_n$, defined as
\begin{align*}
\mathrm{Bias}(x_0) &\coloneqq {\mathbb{E}}\left[\widetilde{f}_n(x_0)\right] - f_0(x_0),\\
\mathrm{Variance}(x_0) &\coloneqq {\mathbb{E}}\left[\left(\widetilde{f}_n(x_0) - {\mathbb{E}}[\widetilde{f}_n(x_0)]\right)^2\right].
\end{align*}
The bias term quantifies how closely the estimator $\widetilde{f}_n$ centers around the true regression function $f_0$ at $x_0$, while the variance term measures the fluctuation of $\widetilde{f}_n$ around its own mean. Balancing bias and variance is a fundamental concern in statistical estimation. Our approach, which employs local polynomial smoothing in the second stage, aims to systematically reduce the bias introduced by the first-stage estimator while controlling variance through appropriate bandwidth selection.
To facilitate our analysis, we introduce the population version of $\widehat{{\mathcal{B}}}_{n}(x_0)$:
\[
\widetilde{{\mathcal{B}}}_{n}(x_0)
\coloneqq \int
\rho\left(\frac{x - x_0}{h_n}\right)
\rho\left(\frac{x - x_0}{h_n}\right)^\top
K\left(\frac{x - x_0}{h_n}\right)
{\mathrm{d}} F(x).
\]
We then establish bounds for the bias and variance terms as follows:
\begin{theorem}[Bias and variance decomposition]
\label{thm:bias_variance}
Let $s, L, C, C_1, C_2, C_3 > 0$ be constants independent of $f_0$ and $n$. Let $h_n$ be the bandwidth of the local polynomial estimator. For every $x \in {\mathcal{X}}$, the following hold:
\begin{itemize}
\item The eigenvalues of $\widetilde{{\mathcal{B}}}_{n}(x)$ are uniformly bounded above and below away from zero.
\item There exists a constant $C > 0$ such that
\[
\int \mathbbm{1}\big[| x - x_0 | \leq h_n\big]{\mathrm{d}} P(x)
\leq C h_n,
\]
where $P(x)$ denotes the distribution of $x$.
\end{itemize}
If $f_0 - \widehat{f}_n$ belongs to the H\"older class $\Sigma(s, L)$ almost surely as $n\to \infty$, then for any $\varepsilon > 0$ and for all $x_0 \in {\mathcal{X}}$, there exists $n_0 > 0$ such that for all $n \geq n_0$, with probability at least $1 - \varepsilon$, it holds that
\begin{align*}
\bigl|\mathrm{Bias}(x_0)\bigr| &\leq C_1h^s_n,\\
\mathrm{Variance}(x_0) &\leq \frac{C_2}{nh_n} + C_3h^{2s}_n.
\end{align*}
\end{theorem}
The smoothness condition on $f_0 - \widehat{f}_n$ can be satisfied in various ways; see Section~\ref{sec:allyouneed} for further discussion. The above result implies that as long as the smoothness of $f_0 - \widehat{f}_n$ is well-controlled and the bandwidth $h_n$ is chosen appropriately to shrink with $n$, both bias and variance remain bounded and ultimately vanish as $n$ grows.
\subsection{Pointwise convergence of the MSE}
From Theorem~\ref{thm:bias_variance}, the bounds on bias and variance imply the convergence of the MSE, noting that
\[
{\mathbb{E}}\left[\left(\widehat{f}_n(x_0) - f_0(x_0)\right)^2\right]
= \mathrm{Bias}(x_0)^2 + \mathrm{Variance}(x_0).
\]
To minimize the MSE, we set
\[
h_n = \alpha n^{-\frac{1}{2s + 1}},
\]
where $\alpha$ is a constant independent of $f_0$ and $n$, $s$ is the smoothness parameter, and $d$ is the dimension of $X$. We then establish the pointwise MSE convergence result.
\begin{theorem}[Pointwise MSE convergence]
\label{thm:pointwise}
Let $s, L, C, \alpha > 0$ be constants independent of $f_0$ and $n$. Let
\[
h_n = \alpha n^{-\frac{1}{2s + 1}}
\]
be the bandwidth of the local polynomial estimator. For every $x_0 \in {\mathcal{X}}$, the following hold:
\begin{itemize}
\item The eigenvalues of $\widetilde{{\mathcal{B}}}_{n}(x)$ are uniformly bounded above and below away from zero.
\item There exists a constant $C > 0$ such that
\[
\int \mathbbm{1}\big[| x - x_0 | \leq h_n\big]{\mathrm{d}} P(x)
\leq C h_n,
\]
where $P(x)$ denotes the distribution of $x$.
\end{itemize}
If $f_0 - \widehat{f}_n$ belongs to the H\"older class $\Sigma(s, L)$ almost surely as $n\to \infty$, then for all $x_0 \in {\mathcal{X}}$, it holds that
\[
{\mathbb{E}}\left[\left(\widetilde{f}_n(x_0) - f_0(x_0)\right)^2\right]
= O\left(n^{-\frac{s}{2s + 1}}\right),
\]
as $n\to\infty$.
\end{theorem}
Theorem~\ref{thm:pointwise} can be extended to establish the convergence of the MSE over the distribution of $X$.
\begin{corollary}[MSE over the distribution of $X$]
Suppose the conditions in Theorem~\ref{thm:pointwise} hold. Then, we have
\[
{\mathbb{E}}\left[\left(\widetilde{f}_n(X) - f_0(X)\right)^2\right] = O\left(n^{-\frac{2s}{2s + 1}}\right),
\]
where the expectation is taken over both $\widetilde{f}_n$ and $X$, as $n \to \infty$.
\end{corollary}
\subsection{Pointwise asymptotic normality}
We now establish the asymptotic normality of $\widetilde{f}(x_0)$. This result follows directly from Theorem~\ref{thm:bias_variance}.
\begin{theorem}[Asymptotic normality]
\label{thm:fastconvergencerates}
Let $s, L, C > 0$ be constants independent of $f_0$ and $n$. Let
$h_n$ be the bandwidth of the local polynomial estimator such that
$h_n \to 0$, $nh_n \to \infty$, and $n h^{2s + 1}_n \to 0$ hold as $n\to \infty$. For every $x \in {\mathcal{X}}$, the following hold:
\begin{itemize}
\item The eigenvalues of $\widetilde{{\mathcal{B}}}_{n}(x)$ are uniformly bounded above and below away from zero.
\item There exists a constant $C > 0$ such that
\[
\int \mathbbm{1}\big[| x - x_0 | \leq h\big]{\mathrm{d}} P(x)
\leq C h_n,
\]
where $P(x)$ is the distribution of $x$.
\end{itemize}
If $f_0 - \widehat{f}_n$ belongs to the H\"older class $\Sigma(s, L)$ almost surely as $n\to \infty$, then for every $x_0 \in {\mathcal{X}}$, it holds that
\[
\sqrt{nh_n}\left(\widetilde{f}(x_0) - f_0(x_0)\right)
\xrightarrow{{\mathrm{d}}} {\mathcal{N}}\big(0, V(x_0)\big),
\]
as $n \to \infty$, where $V(x_0) \geq 0$ is a constant depending on $x_0$ but independent of $n$.
\end{theorem}
Here, note that the following holds from the bias and variance decomposition:
\[
V(x_0) \coloneqq nh_n\mathrm{Variance}(x_0) \leq C_2 + o(h_n).
\]
The asymptotic normality result facilitates statistical inference, such as constructing confidence intervals or hypothesis tests for $f_0(x_0)$. In practice, knowing that $\sqrt{nh_n}(\widetilde{f}(x_0) - f_0(x_0))$ converges in distribution to a normal allows for classical inferential techniques, provided that the variance $V(x_0)$ can be consistently estimated. However, constructing a confidence interval can be challenging, because $s$ (the smoothness parameter) is often unknown and affects the choice of $h_n$. In real applications, data-driven methods, including cross-validation, can be employed to select $h_n$, and bootstrap or other resampling approaches can help in variance estimation.
\subsection{Uniform convergence}
Finally, we establish the uniform (sup-norm) convergence of $\widetilde{f}_n(x_0)$. Unlike mean squared convergence, which depends on the distribution of $X$, uniform convergence bounds the supremum of the estimation error over all $x \in {\mathcal{X}}$. This property is particularly valuable for addressing the covariate shift problem \citep{Shimodaira2000,SchmidtHieber2024}, as uniform convergence remains unaffected by changes in the distribution of $X$ between training and test datasets.
\begin{theorem}[Uniform convergence]
\label{thm:uniform}
Let $s, L_K, L, C, \alpha > 0$ be constants independent of $f_0$ and $n$. Let
\[
h_n \coloneqq \alpha\left(\frac{\log(n)}{n}\right)^{\frac{1}{2s + 1}}
\]
be the bandwidth of the local polynomial estimator. For every $x_0 \in {\mathcal{X}}$, the following hold:
\begin{itemize}
\item The eigenvalues of $\widetilde{{\mathcal{B}}}_{n}(x)$ are uniformly bounded above and below away from zero.
\item There exists a constant $C > 0$ such that
\[
\int \mathbbm{1}\big[| x - x_0 | \leq h_n\big]{\mathrm{d}} P(x)
\leq C h_n,
\]
where $P(x)$ denotes the distribution of $x$.
\item The noise $\varepsilon_i$ is sub-Gaussian, satisfying
\[
{\mathbb{E}}\left[\exp(\lambda \varepsilon_i)\right]
\leq \exp\left(K^2\lambda^2\right)
\quad \forall \lambda \in {\mathbb{R}},
\]
for some constant $K > 0$.
\item The kernel function $K$ is Lipschitz: $K \in \Sigma(1, L_K)$ on ${\mathcal{X}}$.
\end{itemize}
If $f_0 - \widehat{f}_n$ belongs to the H\"older class $\Sigma(s, L)$ almost surely as $n\to \infty$, then the following holds:
\[
{\mathbb{E}}\left[\left\|\widetilde{f}_n - f_0\right\|^2_\infty\right] = O\left(\frac{\log(n)}{n}\right)^{\frac{2s}{2s + 1}}.
\]
\end{theorem}
Although uniform convergence ensures robust performance under distribution shifts, it has not been established for many nonparametric regression estimators derived from modern machine learning methods. One reason for this gap is their data-adaptive nature and the reliance on empirical process techniques, which primarily focus on population risk measures whose expectation is taken over the covariate distribution rather than on sup-norm bounds. For example, \citet{SchmidtHieber2024} analyzes a neural network-based estimator in a one-dimensional covariate setting, demonstrating restricted uniform optimality. In contrast, our approach is model-free and accommodates multi-dimensional covariates, offering broader applicability and greater robustness.
\subsection{Double robustness}
Our proposed debiased estimator also possesses the property of double robustness: if either $\widehat{b}_n$ or $\widehat{f}_n$ is consistent, then $\widetilde{f}_n$ remains consistent.
Recall that
\[
\widetilde{f}_n(x_0)
= \frac{1}{m}\sum_{i \in {\mathcal{I}}^{(2)}}
\left(Y_i - \widehat{f}_n(X_i)\right)\widehat{w}_h(X_i, x_0)
+ \widehat{f}_n(x_0).
\]
Suppose $\widehat{f}_n(x_0)$ converges in probability to $f^\dagger(x_0)$ and $\widehat{b}_n(x_0)$ converges in probability to $b^\dagger(x_0)$. We analyze two cases:
\begin{itemize}
\item If $f^\dagger(x_0) = f_0(x_0)$ but $b^\dagger(x_0) \neq f_0(x_0) - f^\dagger(x_0)$, then
\[
\widetilde{f}_n(x_0)
= \frac{1}{m}\sum_{i \in {\mathcal{I}}^{(2)}}\Big(Y_i - f_0(X_i)\Big)w^\dagger(X_i, x_0)
+ f_0(x_0) + o_p(1),
\]
where $w^\dagger(X_i, x_0)$ depends only on $X_i$ and $x_0$. Since the first term converges in probability to zero, it follows that $\widetilde{f}_n(x_0) \xrightarrow{{\mathrm{p}}} f_0(x_0)$.
\item If $f^\dagger(x_0) \neq f_0(x_0)$ but $b^\dagger(x_0) = f_0(x_0) - f^\dagger(x_0)$, then
\[
\widetilde{f}_n(x_0)
= \left(f_0(x_0) - f^\dagger(x_0)\right) + f^\dagger(x_0) + o_p(1)
= f_0(x_0) + o_p(1).
\]
\end{itemize}
In both cases, $\widetilde{f}_n(x_0) \xrightarrow{{\mathrm{p}}} f_0(x_0)$. We summarize this result as follows.
\begin{theorem}[Double robustness]
For any $x_0 \in {\mathcal{X}}$, if either $\widehat{b}_n(x_0)$ or $\widehat{f}_n(x_0)$ is consistent, then
\[
\widetilde{f}_n(x_0) \xrightarrow{{\mathrm{p}}} f_0(x_0)\qquad \text{as } n\to \infty.
\]
\end{theorem}
Double robustness thus provides an additional safeguard against model misspecification or suboptimal performance of one of the estimators. In practice, this property is particularly valuable when the first-stage model $\widehat{f}_n$ fails to accurately capture the true function, but the second-stage residual model $\widehat{b}_n$ remains reliable—or vice versa.
\section{Smoothness is all you need}
\label{sec:allyouneed}
As demonstrated above, our method relies solely on the smoothness of $\widehat{f} - f_0$. This assumption holds when both $f_0$ and $\widehat{f}$ are individually smooth, but it also remains valid in cases where neither $f_0$ nor $\widehat{f}$ is smooth, as long as their difference is smooth. Neural networks satisfy this condition under the analysis in \citet{SchmidtHieber2020}.
Notably, our results require only the smoothness of the first-stage estimator (i.e., the difference $\widehat{f}_n - f_0$). This assumption significantly broadens the applicability of our method. Even when $f_0$ and $\widehat{f}_n$ are not individually smooth, our approach can still ensure MSE convergence, asymptotic normality, and uniform convergence, provided that the difference is smooth.
From a different perspective, this flexibility allows one to use state-of-the-art or specialized methods in the first stage without compromising the theoretical guarantees of the final estimator, as long as the first-stage error exhibits a level of smoothness akin to the H\"older condition. This feature is particularly useful when the true regression function $f_0$ has a complex or partially known structure, yet the difference $f_0 - \widehat{f}_n$ remains more regular due to design choices or inductive biases in $\widehat{f}_n$.
\section{Conclusion}
This study developed a debiased estimator for nonparametric regression. For any smooth regression estimator, our debiased approach ensures pointwise and uniform MSE convergence, asymptotic normality, and double robustness. The asymptotic normality facilitates statistical inference, while uniform convergence provides distributional robustness for nonparametric regression. Additionally, if either the first-stage or the second-stage nonparametric estimator is consistent, then the resulting regression estimator remains consistent, exhibiting the doubly robust property.
Notably, our estimator requires only the smoothness of the first-stage regression estimator. Specifically, if either (a) both the first-stage estimator $\widehat{f}$ and the true regression function $f_0$ are smooth, or (b) the difference $\widehat{f} - f_0$ is smooth, then the desirable properties of the debiased estimator hold. Thus, our method imposes minimal structural constraints while guaranteeing strong theoretical results, making it an attractive and flexible approach for a wide range of nonparametric regression problems.
\bibliography{Bibtex/citation,Bibtex/causalinference,Bibtex/rdd,Bibtex/semiparametric,Bibtex/nonparametric,Bibtex/highdimensional,Bibtex/experimentaldesign,Bibtex/reinforcement,Bibtex/machinelearning,Bibtex/neuralnet,Bibtex/statistics}
\bibliographystyle{tmlr}
\clearpage