Minimum discrepancy principle strategy for choosing $k$ in $k$-NN regression
We present a novel data-driven strategy to choose the hyperparameter $k$ in the $k$-NN regression estimator. We treat the problem of choosing the hyperparameter as an iterative procedure (over $k$) and propose using an easily implemented in practice strategy based on the idea of early stopping and the minimum discrepancy principle. This model selection strategy is proven to be minimax-optimal, under the fixed-design assumption on covariates, over some smoothness function classes, for instance, the Lipschitz functions class on a bounded domain. The novel method often improves statistical performance on artificial and real-world data sets in comparison to other model selection strategies, such as the Hold-out method and 5-fold cross-validation. The novelty of the strategy comes from reducing the computational time of the model selection procedure while preserving the statistical (minimax) optimality of the resulting estimator. More precisely, given a sample of size $n$, assuming that the nearest neighbors are already precomputed, if one should choose $k$ among $\left\{ 1, \ldots, n \right\}$, the strategy reduces the computational time of the generalized cross-validation or Akaike's AIC criteria from $\mathcal{O}\left( n^3 \right)$ to $\mathcal{O}\left( n^2 (n - k) \right)$, where $k$ is the proposed (minimum discrepancy principle) value of the nearest neighbors. Code for the simulations is provided atthis https URL.
\section{Introduction}
\label{sec:intro}
Nonparametric regression estimation is a fundamental problem in statistics and machine learning. The $k$-NN regression estimator (\cite{gyorfi2006distribution, biau2015lectures}) is a very simple and popular choice in practice. For this estimator, the central issue is choosing properly the number of neighbors $k$.

The theoretical performance of the $k$-NN regression estimator has been widely studied since the 1970s (\cite{devroye1978uniform, collomb1979estimation, devroye1981almost, bhattacharya1987weak, biau2010rates, kpotufe2011k, biau2015lectures, zhao2021minimax}). For example, in (\cite[Chapter 12]{biau2015lectures}) the uniform consistency of the $k$-NN estimator is proved under the condition that $k(n) / n \to 0$ as $n \to \infty$, where $n$ is the sample size. However, as it was shown in \cite{gyorfi2006distribution}, the nearest neighbor estimator ($k = 1$) is proved to be consistent only in the noiseless case. 

Recently, researchers started to be interested in choosing $k$ from the data (\cite{gyorfi2006distribution, arlot2009data, kpotufe2011k, azadkia2019optimal}). Apparently, the most common and simplest strategy to choose $k$ is to assume some smoothness assumption on the regression function (e.g., the Lipschitz condition \cite{gyorfi2006distribution}) and to find $k$ that makes an upper bound on the bias and the variance of the $k$-NN regression estimator equal. This method has a clear lack: one needs to know the smoothness of the regression function (e.g., the Lipschitz constant). \cite{2009arXiv0909.1884A} gave a \textit{data-driven strategy} for choosing a hyperparameter for different linear estimators (e.g., the $k$-NN estimator) based on the idea of minimal penalty, introduced previously in \cite{birge2007minimal}. The main inconvenience of this strategy is that one needs to compute all the estimators $\mathbf{F}_n = \{ f^k, \ k = 1, \ldots, n \}$ of the regression function in order to choose the optimal one among them by comparing them via a special criterion that involves the empirical error (least-squares loss). To list other (similar) strategies, one can think about the Akaike's AIC (\cite{akaike1998information}), Mallows's $C_p$ (\cite{mallows2000some}) criteria or generalized cross-validation (\cite{li1987asymptotic, hastie2009elements}), where one has to compute the empirical risk error plus a penalty term for any $k = 1, \ldots, n$. Often it is computationally expensive and restricts the use in practice. This gives rise to the problem of choosing the hyperparameter "in real-time", meaning that the practitioner should compute iteratively $f^k \in \mathbf{F}_n$. Eventually, this iterative process has to be stopped. This problem can be solved by applying the \textit{early stopping rule} (\cite{engl1996regularization, zhang2005boosting, Yao2007, raskutti2014early, blanchard2018optimal, wei2017early}).





The first early stopping rule that could be potentially data-driven was proposed by (\cite{blanchard2018optimal, blanchard2018early, celisse2021analyzing}) for spectral filter iterative algorithms (see, e.g., \cite{bauer2007regularization, gerfo2008spectral} for examples of such algorithms). The idea behind the construction of this early stopping rule is the so-called \textit{minimum discrepancy principle} (MDP) that is based on finding a first iteration for which a learning algorithm starts to fit the noise. The key quantity for the analysis of the minimum discrepancy principle is the \textit{empirical risk error} (the train error in the terminology of the machine learning community \cite{goodfellow2016deep}), which is monitored throughout the whole learning process. The process thus is stopped if the empirical risk starts to fit the noise.


\textbf{Contribution.} In the present paper, we propose applying the minimum discrepancy principle stopping rule for the $k$-NN regression estimator in order to select $k$. We prove a non-asymptotic upper bound on the performance of the minimum dicrepancy principle stopping rule measured in the empirical $L_2(\mathbb{P}_n)$ norm. This bound implies that, under a quite mild assumption on the regression function, the minimum discrepancy principle stopping rule provides a statistically (minimax) optimal functional estimator, in particular, over the class of Lipschitz functions on a bounded domain. In Section 5, we show that our method improves the performances of classical selection procedures, such as 5-fold cross-validation and Hold-out, on artificial and real data. Besides that, the proposed strategy lowers the computational time of the selection procedure compared to some well-known model selection strategies, such as the generalized cross-validation and Akaike's AIC, i.e., $\mathcal{O}\left( n^2 \left( n - k \right) \right)$ for MDP($k$) and $\mathcal{O}\left( n^3 \right)$ for Akaike's AIC and generalized cross-validation, respectively. We emphasize that the proposed strategy for choosing $k$ is \textit{easy-to-implement} in practice since it involves only computing the empirical risk and estimating the noise level of the regression model.   


\textbf{Outline of the paper.} The organization of the paper is as follows. Section \ref{sec:model} describes the statistical model, its main assumption and introduces the notation that will be used along the paper. In Section \ref{sec:estimator_stopping_rule}, we introduce the $k$-NN estimator and explain how to compute the minimum discrepancy early stopping rule. Section \ref{sec:optimality_result} provides the main theoretical result that shows that the proposed rule is statistically optimal for some classes of functions (e.g., the well-known class of Lipschitz functions on a bounded domain). In Section \ref{sec:simulations}, one can find simulation results for the proposed stopping rule. To be precise, we compare this rule to the generalized cross-validation estimator, Akaike's AIC criterion, $V$--fold and Hold-out cross-validation stopping rules (\cite{arlot2010survey}) tested on some artificial and real-world data sets. Section \ref{sec:conc} concludes the paper. All the technical proofs are in the supplementary material.


\section{Statistical model, main assumption, and notation} \label{sec:model}

In the nonparametric regression setting, one works with a sample $(x_1, y_1), \ldots, \\ (x_n, y_n) \in \mathcal{X}^n \times \mathbb{R}^n$ that satisfies the statistical model
\begin{equation} \label{nonparametric_regression_model}
    y_i = f^*(x_i) + \varepsilon_i, \ \ i = 1,\ldots, n,
\end{equation}
where $f^*:\mathcal{X} \mapsto \mathbb{R}$, $\mathcal{X} \subseteq \mathbb{R}^d$, is a measurable function on some set $\mathcal{X}$, and $\{ \varepsilon_i \}_{i=1}^n$ are i.i.d. Gaussian noise variables $\mathcal{N}(0, \sigma^2)$. On the one hand, in the theoretical analysis in Section \ref{sec:optimality_result}, we assume that the parameter $\sigma^2 > 0$ is fixed and known. On the other hand, from the practical perspective in Sections \ref{artificial_data} and \ref{real_data}, we estimate $\sigma^2$ and provide experimental results on artificial and real-world data sets. One should point out here that the assumption of known variance $\sigma^2$ is quite typical in the (theoretical) model selection literature with nonparametric regression (see, e.g., \cite{li1986asymptotic, li1987asymptotic, yang2005can, cao2006oracle}). In addition to that, we assume that $\{ x_i \in \mathcal{X} \}_{i=1}^n$ are \textit{fixed} covariates (corresponds to the so-called fixed design setting), thus we observe noise only in the responses $\{ y_i \}_{i=1}^n$.  The goal of the present paper is to estimate optimally the regression function $f^*$. The term "optimally" will be explained in Section \ref{sec:estimator_stopping_rule}.


In the context of the \textit{fixed design} setting, the performance of an estimator $\widehat{f}$ of $f^*$ is measured in terms of the so-called \textit{empirical norm} defined as 
\begin{equation}
    \lVert \widehat{f} - f^* \rVert_n^2 \coloneqq \frac{1}{n}\sum_{i=1}^n \left[ \widehat{f}(x_i) - f^*(x_i) \right]^2,  
\end{equation}
where $\lVert h \rVert_n \coloneqq \sqrt{1/n\sum_{i=1}^n h(x_i)^2}$ for any bounded on $\mathcal{X}$ function $h$. We denote the empirical norm as $L_2(\mathbb{P}_n)$. For each bounded over $\mathcal{X}$ functions $h_1, h_2$, $\langle h_1, h_2 \rangle_n$ denotes the related inner product defined as $\langle h_1, h_2 \rangle_n \coloneqq 1/n \sum_{i=1}^n h_1(x_i)h_2(x_i)$.
Further, $\mathbb{P}_{\varepsilon}$ and $\mathbb{E}_{\varepsilon}$ denote the probability and expectation with respect to $\{ \varepsilon_i \}_{i=1}^n$.

\textbf{Notation.}
Throughout the paper, $\lVert \cdot \rVert$ and $\langle \cdot, \cdot \rangle$ are the usual Euclidean norm and related inner product. $\lVert M \rVert_2$ and $\lVert M \rVert_{F}$ signify the operator and Frobenius norms of the matrix $M \in \mathbb{R}^{n \times n}$, respectively. We denote the trace of the matrix $M$ by $\textnormal{tr}(M)$. In addition to that, $\mathbb{I}\left\{ \mathcal{E} \right\}$ is equal to $1$ if the probabilistic event $\mathcal{E}$ holds true, otherwise it is equal to $0$. For $a \geq 0$, we denote by $\floor*{a}$ the largest natural number that is smaller than or equal to $a$. We denote by $\ceil*{a}$ the smallest natural number that is greater than or equal to $a$. Along the paper, $I_n$ is the identity matrix of size $n \times n$.


We make the following assumption on the regression function $f^*$ introduced earlier in Eq. (\ref{nonparametric_regression_model}).
\begin{assumption}[Boundness of the r.f.] \label{assumption_boundness}
$f^*$ is bounded on $\mathcal{X}$, meaning that there exists a constant $\mathcal{M} > 0$ such that
\begin{equation}
    \left| f^*(x) \right| \leq \mathcal{M} \qquad \textnormal{ for all }x \in \mathcal{X}.
\end{equation}
\end{assumption}

Assumption \ref{assumption_boundness} is quite standard in the nonparametric regression literature (\cite{gyorfi2006distribution, zhao2021minimax}). In particular, Assumption \ref{assumption_boundness} holds when the set $\mathcal{X}$ is bounded, and the regression function $f^*$ is $L$-Lipschitz with some positive constant $L$ (cf. \cite{gyorfi2006distribution}).




Along the paper, we use the notation $c, c_1, C, \widetilde{c}, \widetilde{C}, \ldots$ to show that numeric constants $c, c_1, C, \widetilde{c}, \widetilde{C}, \ldots$ can depend only on $d, \sigma$, and $\mathcal{M}$. The values of all the constants may change from line to line or even in the same line. 




\section{$k$-NN estimator and minimum discrepancy stopping rule} \label{sec:estimator_stopping_rule}

\subsection{$k$-NN regression estimator}

Let us transform the initial nonparametric regression model (\ref{nonparametric_regression_model}) into its vector form 
\begin{equation} \label{nonparametric_regression_vector}
    Y \coloneqq \left[y_1, \ldots, y_n\right]^\top = F^* + \varepsilon \in \mathbb{R}^n, 
\end{equation}
where the vectors $F^* \coloneqq [f^*(x_1), \ldots, f^*(x_n)]^\top$ and $\varepsilon \coloneqq [\varepsilon_1, \ldots, \varepsilon_n]^\top$.

Define a $k$-nearest neighbor estimator $f^k$ of $f^*$ from (\ref{nonparametric_regression_model}) at the point $x_i, \ i = 1, \ldots, n,$ as
\begin{equation} \label{estimator}
    f^k(x_i) \coloneqq F^k_i = \frac{1}{k} \sum_{j \in \mathcal{N}_k(i)} y_j, \qquad k = 1, \ldots, n,
\end{equation}
where $\mathcal{N}_k(i)$ denotes the indices of the $k$ nearest neighbors of $x_i$ among $\{1, \ldots, n \}$ in the usual Euclidean norm in $\mathbb{R}^d$, where ties are broken at random. In words, in Eq. (\ref{estimator}) one weights by $1/k$ the response $y_j$ if $x_j$ is a $k$ nearest neighbor of $x_i$ measured in the Euclidean norm. Note that other adaptive metrics (instead of the Euclidean one) have been also considered in the literature (\cite[Chap. 14]{hastie2009elements}).


One can notice that the $k$-NN regression estimator (\ref{estimator}) belongs to the class of linear estimators (\cite{2009arXiv0909.1884A, hastie2009elements}), meaning that $F^k \in \mathbb{R}^n$ estimates the vector $F^*$ as it follows.
\begin{equation} \label{main_estimator}
    F^k \coloneqq \left( f^k(x_1), \ldots, f^k(x_n) \right)^\top = A_k Y,
\end{equation}
where $A_k \in \mathbb{R}^{n \times n}$ is the matrix described below.
\begin{equation} \label{a_k_matrix}
    \begin{cases}
    \forall 1 \leq i, j \leq n, \ \left( A_k \right)_{ij} \in \{0, 1/k\} \textnormal{ with } k \in \{1, \ldots, n \},\\
    \forall 1 \leq i \leq n, \ \left( A_k \right)_{ii} = 1/k \textnormal{ and } \sum_{j=1}^n \left( A_k \right)_{ij} = 1.
    \end{cases}    
\end{equation}
\vspace{0.05cm}
Saying differently, $(A_k)_{ij} = 1/k$ if $x_j$ is a $k$ nearest neighbor of $x_i$, otherwise $(A_k)_{ij} = 0, \ i, j \in \{1, \ldots, n\}$.

Define the mean-squared error (the risk error) of the estimator $f^k$ as
\begin{equation} \label{mean-squared-error}
    \textnormal{MSE}(k) \coloneqq \mathbb{E}_{\varepsilon} \lVert f^k - f^* \rVert_n^2 = \frac{1}{n} \mathbb{E}_{\varepsilon} \sum_{i=1}^n \Big( \frac{1}{k}\sum_{j \in \mathcal{N}_k(i)}y_j - f^*(x_i) \Big)^2.
\end{equation}
Further, we will introduce the (squared) bias and variance of the functional estimator $f^k$ (see, e.g., \cite[Eq. (7)]{2009arXiv0909.1884A}),
\begin{equation}
    \textnormal{MSE}(k) = B^2(k) + V(k),
\end{equation}
where 
\begin{equation*}
    B^2(k) = \lVert (I_n - A_k)F^* \rVert_n^2, \qquad V(k) = \frac{\sigma^2}{n} \textnormal{tr}\left( A_k^\top A_k \right) = \frac{\sigma^2}{k}.
\end{equation*}

Thus, the variance term $\sigma^2/k$ is a decreasing function of $k$. Note that $B^2(1) = 0, \ V(1) = \sigma^2$, and $B^2(n) = (1 - 1/n)^2\lVert f^* \rVert_n^2, \ V(n) = \sigma^2/n$. Importantly, the bias term $B^2(k)$ can have \textit{arbitrary} behavior on the interval $[1, n]$. 

\vspace{0.05cm}

Ideally, we would like to minimize the mean-squared error (\ref{mean-squared-error}) as a function of $k$. However, since the bias term is not known (it contains the unknown regression function), one should introduce other quantities that will be related to the bias. In our case, this quantity will be the \textit{empirical risk} at $k$:
\begin{equation} \label{empirical_risk}
    R_k \coloneqq \lVert (I_n - A_k) Y \rVert_n^2.
\end{equation}

$R_k$ measures how well the estimator $f^k$ fits $Y$. Remark that $R_1 = 0$ (corresponds to the "overfitting" regime) and $R_n = (1 - 1/n)^2 \frac{1}{n} \sum_{i=1}^n y_i^2$ (corresponds to the "underfitting" regime), but there is no information about the monotonicity of $R_k$ on the interval $[1, n]$. 

Furthermore, some information about the bias is contained in the expectation (over the noise $\{ \varepsilon_i \}_{i=1}^n$) of the empirical risk. To be precise, for any $k \in \{ 1, \ldots, n\}$,
\begin{align} \label{exp_empirical_risk}
    \begin{split}
    \mathbb{E}_{\varepsilon}R_k &= \sigma^2 + B^2(k) - \frac{\sigma^2(2 \textnormal{tr}(A_k) - \textnormal{tr}(A_k^\top A_k))}{n}\\
   
    &= \sigma^2 + B^2(k) - V(k).
    \end{split}
\end{align}

Let us illustrate all the mentioned quantities in one example in Fig \ref{fig:bvr}. We take the regression function equal to $f^*(x) = \lVert x - 0.5 \rVert / \sqrt{3} - 0.5$ and the noise variance $\sigma^2 = 0.01$. We take $n = 50$, the uniform covariates $x_i \overset{\text{i.i.d.}}{\sim} \mathbb{U}[0, 1]^3$, and plot the bias term $B^2(k)$, the variance term $V(k)$, risk error $\textnormal{MSE}(k)$, empirical risk $R_k$, and its expectation $\mathbb{E}_{\varepsilon}R_k$ versus the number of neighbours $k$.


Note that among all defined quantities, only the variance term $V(k)$ can be proved monotonic (without an additional assumption on the smoothness of $f^*$). Importantly, Fig \ref{fig:bvr} indicates that choosing $k = 6$ will provide the user with the global optimum of the risk (the mean-squared error) curve. Thus, for instance, it would be meaningless (according to the risk curve) to compute all the estimators $f^k$ (\ref{estimator}) for $k = 1, \ldots, 6$. 


\begin{figure} 
\centering
\begin{tikzpicture}[scale=0.85]
\tikzstyle{every node}=[font=\normalsize]
\begin{axis}[
    title={},
    xlabel={Number of neighbors},
    ylabel={Value},
    xmin=1, xmax=12,
    ymin=0, ymax=0.02,
    xtick={1,3,5,7, 10},
    ytick={0.002, 0.006, 0.01, 0.02},
    legend pos=north west,
    ymajorgrids=true,
    xmajorgrids=true,
    grid style=dashed,
    legend style={nodes={scale=0.6, transform shape}}, 
]

\addplot[
    color=blue,
    mark=square,
    ]
    coordinates {
    (1,0)(2,6.29080327e-03)(3,8.47616891e-03)(4,1.08716416e-02)(5,1.12560012e-02)(6,1.14496210e-02)(7,1.37171039e-02)(8,1.57263962e-02)(9,1.83168233e-02)(10,1.78368289e-02)(11,1.94323757e-02)(12,1.98591401e-02)
    };
    \addlegendentry{Empirical risk}
   
\addplot[
    color=red,
    mark=+,
    ]
    coordinates {
    (1,0)(2,1.48163079e-03)(3,1.91217442e-03)(4,1.90325965e-03)(5,2.31901908e-03)(6,2.42449703e-03)(7,3.02981110e-03)(8,3.27183370e-03)(9,4.00873281e-03)(10,4.53268739e-03)(11,5.57342952e-03)(12,6.21032204e-03)
    };
    \addlegendentry{Bias}
\addplot[
    color=violet,
    mark=star,
    ]
    coordinates {
    (1,0.01)(2,0.005)(3,0.00333333)(4,0.0025)(5,0.002)(6,0.00166667)(7,0.00142857)(8,0.00125)(9,0.00111111)(10,0.001)(11,0.00090909)(12,0.00083333)
    };
    \addlegendentry{Variance}
\addplot[
    color=orange,
    mark=x,
    ]
    coordinates {
    (1,0.01)(2,0.00648163)(3,0.00524551)(4,0.00440326)(5,0.00431902)(6,0.00409116)(7,0.00445838)(8,0.00452183)(9,0.00511984)(10,0.00553269)(11,0.00648252)(12,0.00704366)
    };
    \addlegendentry{Risk}
\addplot[
    color=black,
    mark=*,
    ]
    coordinates {
    (1,0)(2,6.48163079e-03)(3,8.57884109e-03)(4,9.40325965e-03)(5,1.03190191e-02)(6,1.07578304e-02)(7,1.16012397e-02)(8,1.20218337e-02)(9,1.28976217e-02)(10,1.35326874e-02)(11,1.46643386e-02)(12,1.53769887e-02)
    };
    \addlegendentry{Exp. empirical risk}
\end{axis}
\end{tikzpicture}
\caption{Sq. bias, variance, risk and (expected) empirical risk behaviour.}
\label{fig:bvr}
\end{figure}




Our main concern is to design a data-driven strategy to choose $\widehat{k} \in \{1, \ldots, n \}$, which can be seen as a mapping from the data $\{ (x_i, y_i) \}_{i=1}^n$ to a positive integer so that the $L_2(\mathbb{P}_n)$-error $\lVert f^{\widehat{k}} - f^* \rVert_n^2$ is as small as possible. More precisely, the goal is to define a data-driven $\widehat{k}$ such that it satisfies the following non-asymptotic upper bound ("oracle-type inequality" \cite{wainwright2019high}):
\begin{equation} \label{or_type_ineq}
     \lVert f^{\widehat{k}} - f^* \rVert_n^2 \leq C_n \mathbb{E}_{\varepsilon} \lVert f^{k_{\textnormal{opt}}} - f^* \rVert_n^2 + r_n
\end{equation}
with high (exponential) probability over $\{ \varepsilon_i \}_{i=1}^n$, where $f^{k_{\textnormal{opt}}}$ is a minimax-optimal estimator of the regression function $f^* \in \mathcal{F}$, $\mathcal{F}$ is some a priori chosen function space. The leading constant $C_n$ should be bounded, and the remainder term $r_n$ is negligible (smaller) with respect to the "optimal risk error" $\mathbb{E}_{\varepsilon} \lVert f^{k_{\textnormal{opt}}} - f^* \rVert_n^2$. 

\subsection{Related work} \label{related_work}
The idea of choosing the hyperparameter $k$ from the data has been already considered in the literature. For example, the classical procedures such as generalized cross-validation (\cite{craven1978smoothing, li1987asymptotic, cao2006oracle}), penalized estimators (\cite{li1987asymptotic, mallows2000some, arlot2009data, arlot2009data_}), and different cross-validation methods (\cite{arlot2010survey}) are popular choices for linear estimators. Let us consider them in more detail.


\textbf{Generalized CV (\cite{li1987asymptotic, cao2006oracle, hastie2009elements}).} This model selection method has been widely studied in the case of (kernel) ridge regression (\cite{craven1978smoothing}) and smoothing splines (\cite{cao2006oracle}). In particular, \cite{cao2006oracle} proved a non-asymptotic oracle inequality for the generalized CV estimator when the variance $\sigma^2$ is known. However, in a more general case, GCV estimates $\sigma^2$ implicitly, which is an advantage of the method. In addition to that, GCV for $k$-NN regression is proved by \cite{li1987asymptotic} to be an asymptotically optimal model selection criterion (i.e., $\lVert f^{k_{\textnormal{GCV}}} - f^* \rVert_n^2 \ / \ \underset{k}{\inf}\lVert f^{k} - f^* \rVert_n^2 \to 1$ in probability when $n \to + \infty$) under the assumption $\lVert A_k \rVert_2 \leq c, \ k = 1, \ldots, n,$ for some constant $c$. It is worth to mention that generalized cross-validation provides an approximation to the so-called leave-one-out cross-validation (\cite{arlot2010survey, celisse2018theoretical}), which is an exhaustive model selection procedure. In the case of GCV, if the nearest neighbors' matrices are already precomputed, its computational time is $\mathcal{O}\left( n^3 \right)$ elementary operations. The GCV strategy will be later considered in our simulations (see Section \ref{sec:simulations}).





\textbf{Penalized estimators} date back to the work on AIC (\cite{akaike1998information}) or  Mallow's $C_p$ (\cite{mallows2000some}) criteria, where a penalty proportional to the dimension of the model is added to the quadratic loss (i.e., \textit{the empirical risk error} in our notation (\ref{empirical_risk})) when the noise level $\sigma^2$ is known. As it was for the GCV strategy, the computational time of AIC and Mallow's $C_p$ are $\mathcal{O}\left( n^3 \right)$. After that, a new approach was developed by \cite{birge2007minimal}, where the authors introduced the so-called "slope heuristics" for projection matrices. This notion was based on the introduction of the penalty $\textnormal{pen}(k) = K \textnormal{tr}(A_k)$, where $\textnormal{tr}(A_k)$ is the dimension of the model, and $K$ is a constant that can depend on $\sigma^2$, in particular. It appeared that there exists a constant $K_{\textnormal{min}}$ such that $2 K_{\min} \textnormal{tr}(A_k)$ yields an asymptotically optimal model selection procedure. This gives rise to some strategies for the estimation of constant $K_{\textnormal{min}}$ from the data, as it was done, for instance, by \cite{2009arXiv0909.1884A} for some linear estimators when $\sigma^2$ is unknown. 


\textbf{Cross-validation methods (\cite{arlot2010survey}).} These model selection methods are the most used in practice. Compared to generalized cross-validation, for instance, $V$--fold cross-validation method (\cite{geisser1975predictive, arlot2010survey}) incurs a large computational cost (with $V$, which is not too small). To be precise, the $V$--fold cross-validation requires the model selection procedure to be performed $V$ times for each value of $k \in \{1, \ldots, n\}$. Another alternative could be the Hold-out method (\cite{wegkamp2003model, arlot2010survey}), which consists in randomly splitting the data into two parts for each value $k \in \{1, \ldots, n\}$:  one is dedicated for training the estimator (\ref{estimator}) and the other one is dedicated for testing (see, e.g., Section \ref{sec:simulations} for more details in a simulated example).


\subsection{Bias-variance trade-off and minimum discrepancy principle rule}

We are at the point to define our first "reference rule". Based on the nonparametric statistics literature (\cite{wasserman2006all, tsybakov2008introduction}), the bias-variance trade-off usually provides an optimal estimator in the minimax sense:
\begin{equation} \label{k_star}
    k^* = \inf \left\{ k \in \{ 1, \ldots, n \} \mid B^2(k) \geq V(k) \right\}.
\end{equation}

The bias-variance trade-off stopping rule $k^*$ does not always exist due to arbitrary behavior of the bias term $B^2(k)$. Thus, if no such $k^*$ exists, set $k^* = n$. If it exists, then $k^* \geq 2$ since $V(1) > B^2(1)$.



Notice that the stopping rule $k^*$ is \textit{not computable} in practice, since it depends on the unknown bias $B^2(k)$. Nevertheless, we can create a data-driven version of $k^*$ by means of the empirical risk $R_k$ from Eq. (\ref{empirical_risk}). 

Eq. (\ref{exp_empirical_risk}) gives us that the event $\{ B^2(k) \geq V(k) \}$ is equivalent to the event $\{ \mathbb{E}_{\varepsilon}R_k \geq \sigma^2 \}$, so we conclude that $k^* = \inf \{ k \in \{ 1, \ldots, n \} \mid \mathbb{E}_{\varepsilon}R_k \geq \sigma^2 \}$. This gives rise to an estimator of $k^*$ that we denote as $k^{\tau}$. This stopping rule is called \textit{the minimum discrepancy principle} stopping rule and is defined as
\begin{equation} \label{k_tau}
    k^{\tau} = \sup \left\{ k \in \{ 1, \ldots, n \} \mid R_k \leq \sigma^2 \right\}. 
\end{equation}


\textbf{Remarks.} Note that in Eq. (\ref{k_tau}), we introduced a supremum instead of the infimum from Eq. (\ref{k_star}). That was done on purpose since there could be several points of the bias-variance trade-off, and apparently, the bias (and the empirical risk) could behave badly in the area "in-between". In order to calculate $k^{\tau}$, the user should, first, compute the empirical risk $R_k$ at $k = n$ (thus, the matrix $A_n$ of $n$ nearest neighbors). After that, one needs to decrease $k$ until the event $\{ R_k \leq \sigma^2 \}$ holds true. It is worth mentioning that it is not necessary to compute explicitly all the matrices $A_k, \ k = n, \ n-1, \ldots$, since, for instance, the matrix $A_{n-1}$ could be easily derived from the matrix $A_n$ (assuming that one has already arranged the neighbors and removed the $n^{\text{th}}$ neighbors from the matrix $A_n$), i.e., 
\begin{equation} \label{iteration_for_matrix}
    [A_{n-1}]_{ij} = \frac{n}{n-1}[A_n]_{ij}, \ \forall i, j \in \{1, \ldots, n\}.
\end{equation}
It is one of the main computational advantages of the proposed rule (\ref{k_tau}). For more details on the efficient computation of the nearest neighbors' matrices, see, e.g., \cite{bentley1975multidimensional, omohundro1989five}. In addition to all of that, we emphasize that the definition (\ref{k_tau}) of $k^{\tau}$ does not require the knowledge of the constant $\mathcal{M}$ from Assumption \ref{assumption_boundness}, and $k^{\tau}$ does not require computing the empirical risk $R_k$ for all values $k = 1, \ldots, n$, while it is the case, for instance, for generalized cross-validation or Mallow's $C_p$ (see Section \ref{sec:simulations}). Moreover, we need to point out that the stopping rule (\ref{k_tau}) depends on the noise level $\sigma^2$, which should be estimated in practice, as for the AIC or Mallow's $C_p$ criteria (\cite{akaike1974new, mallows2000some, hastie2009elements}). We will provide a consistent estimator of $\sigma^2$ in Section \ref{real_data} while dealing with real-world data. Regarding the computational time of $k^{\tau}$, if the nearest neighbors' matrices are already computed, it is of the order $\mathcal{O}\left( n^2 \left( n - k^{\tau} \right) \right)$, which is less than $\mathcal{O}\left( n^3 \right)$ for the AIC/Mallow's $C_p$ criteria, or GCV.


\vspace{0.2cm}


There is a large amount of literature (\cite{engl1996regularization, bauer2007regularization, blanchard2018optimal, blanchard2018early, celisse2021analyzing}) on the minimum discrepancy principle for spectral filter algorithms such as gradient descent, ridge regularization, and spectral cut-off regression, and providing an exhaustive review on this strategy is out of the scope of the paper (e.g., \cite{blanchard2018optimal, celisse2021analyzing} provide a thorough review). We should, however, emphasize that intuitively the minimum discrepancy principle determines the first time at which a learning algorithm starts to fit noise, which is measured by $\sigma^2$ in the present context. 

Moreover, one is able to notice that, if the empirical risk is close to its expectation, $k^{\tau}$ should produce an optimal estimator in some sense. The main question that should be asked is "In which setting is it possible to quantify this gap between $R_k$ and $\mathbb{E}_{\varepsilon}R_k$ that will not be statistically large?". This question is the main technical obstacle of the present paper. In what follows, we show that for a quite large class of regression functions, $k^{\tau}$ is optimal in the sense of Ineq. (\ref{or_type_ineq}). 




\section{Theoretical optimality result} \label{sec:optimality_result} 

Let us start to describe the main theoretical result of the present paper. The following theorem applies to the estimator defined in Eq. (\ref{main_estimator}). 
\begin{theorem} \label{main_th}
Under Assumption \ref{assumption_boundness}, for arbitrary $u \geq 0$,
    \begin{equation} \label{main_ineq}
        \lVert f^{k^{\tau}} - f^* \rVert_n^2 \leq 8 V(k^*) + C_1 \left( \frac{u}{n} + \frac{\sqrt{u}}{\sqrt{n}} \right) + C_2 \sqrt{\frac{\log n}{n}}
\end{equation}
    with probability at least $1 - 16 \exp \left(- u \right)$, where positive constants $C_1, C_2$ can depend on $d$, $\sigma$, and $\mathcal{M}$.
    
Moreover, if $k^*$ from Eq. (\ref{k_star}) exists, then for arbitrary $u \geq 0$, 
\begin{equation} \label{corollary_from_the_main}
    \lVert f^{k^{\tau}} - f^* \rVert_n^2 \leq \underbrace{4 \ \textnormal{MSE}(k^*)}_{\textnormal{Main term}} + \underbrace{C_1 \left( \frac{u}{n} + \frac{\sqrt{u}}{\sqrt{n}} \right) + C_2 \sqrt{\frac{\log n}{n}}}_{\textnormal{Remainder term}}
\end{equation}
with probability at least $1 - 16 \exp (-u)$, where constants $C_1, C_2$ are from Ineq. (\ref{main_ineq}).
\end{theorem}

\begin{proof}
The full proof is deferred to the supplementary material. Let us provide a sketch of the proof here.


The main ingredients of the proof are two deviation inequalities (cf. Corollary \ref{corollary_for_variance} and Lemma \ref{deviation_bias_lemma} in the supplementary material): for any $x \geq 0$,
\begin{equation}
    \mathbb{P}_{\varepsilon} \left( V(k^{\tau}) > 2 V(k^*) + x \right) \leq 2 \exp \left( -c_d n \min \left( \frac{x}{\sigma^2}, \frac{x^2}{\sigma^4} \right) \right),
\end{equation}
and 
\begin{equation} \label{bias_k_tau_variance_k_star}
    B^2(k^{\tau}) \leq 2 V(k^*) + c_1 \sqrt{\frac{\log n}{n}} + 2x,
\end{equation}
where Ineq. (\ref{bias_k_tau_variance_k_star}) holds with probability at least $1 - 12 \exp \left( -c n \min \left( x^2, x \right) \right)$.


After that, one can split the $L_2(\mathbb{P}_n)$-error at $k^{\tau}$ into two parts:
\begin{equation*}
    \lVert f^{k^{\tau}} - f^* \rVert_n^2 \leq 2 B^2(k^{\tau}) + 2 \lVert A_{k^{\tau}}\varepsilon \rVert_n^2.    
\end{equation*}

It is sufficient to derive high probability control of $\underset{k}{\sup} \left| \lVert A_{k}\varepsilon \rVert_n^2 - V(k)   \right|$ for $k = 1, \ldots, n$ (see Appendix \ref{var_control} in the supplementary material). That was the reason why the term $\mathcal{O}\left(\sqrt{\frac{\log n}{n}}\right)$ appeared in Eq. (\ref{main_ineq}).


Finally, one can apply $V(k^*) \leq \frac{1}{2}\textnormal{MSE}(k^*)$, if $k^*$ exists, and $u = c n \min \left( x^2, x \right)$. The claim follows.

\end{proof}
In order to gain some intuition of the claim of Theorem \ref{main_th}, let us make some comments.


First of all, Ineq. (\ref{corollary_from_the_main}) is non-asymptotic, meaning that it holds true for any sample size $n \geq 1$. Second, Ineq. (\ref{corollary_from_the_main}) holds "with high probability", which is a stronger result than in expectation since (\cite{li1987asymptotic}) there are model selection procedures that are asymptotically optimal (when $n \to +\infty$) with high probability but not in expectation. 

Third, the main term in Ineq. (\ref{corollary_from_the_main}) is the risk error at the bias-variance trade-off times 4 (this constant could be improved). Ideally, one should rather introduce the oracle risk $\underset{k=1, \ldots, n}{\inf} \left[ \mathbb{E}_{\varepsilon} \lVert f^k - f^* \rVert_n^2 \right]$ and compare $\lVert f^{k^{\tau}} - f^* \rVert_n^2$ to it. However, to the best of our knowledge, a smoothness assumption is needed to connect the bias-variance trade-off risk and the oracle risk. That was the reason to keep the main term as it was stated. Fourth, the right hand side term of Ineq. (\ref{corollary_from_the_main}) is of the order $\mathcal{O}\left(\sqrt{\frac{\log n}{n}}\right)$. Notice that the same rate for this term was achieved in \cite{azadkia2019optimal} but in terms of the expectation over the noise. 




A natural question would be to understand if the rate $\mathcal{O}\left(\sqrt{\frac{\log n}{n}}\right)$ in Ineq. (\ref{main_ineq}) or (\ref{corollary_from_the_main}) is sufficiently fast. In order to do that, one should precise the function space $\mathcal{F}$, where $f^*$ lies in. In what follows, we will mention one famous example (among others) of a such function space $\mathcal{F}$.
\begin{example}
Consider the class of functions 
\begin{equation} \label{lipschitz}
    \mathcal{F}_{\textnormal{Lip}}(L) \coloneqq \left\{ f: [0, 1]^d \mapsto \mathbb{R} \mid f(0) = 0, \ f \textnormal{ is } L-\textnormal{Lipschitz} \right\},
\end{equation}
where $f$ is $L$-Lipschitz means that $|f(x) - f(x^\prime)| \leq L \lVert x - x^{\prime} \rVert$ for all $x, x^{\prime} \in [0, 1]^d$. In this case (see, e.g., \cite[Theorem 3.2]{gyorfi2006distribution} with $p = 1$), 
\begin{equation} \label{minimax_lipschitz}
    \mathbb{E}_{\varepsilon} \lVert \widehat{f} - f^* \rVert_n^2 \geq c_l n^{-\frac{2}{2 + d}},
\end{equation}
for some positive constant $c_l$, for any measurable of the input data $\widehat{f}$.
\end{example}

Therefore, for the class of $L$-Lipshitz functions, the rate $\mathcal{O}(\sqrt{\log n / n})$ is faster than the minimax-optimal rate $\mathcal{O}(n^{-\frac{2}{2+d}})$ for any $d > 2$.


As for the main term $8 V(k^*)$ in Ineq. (\ref{main_ineq}), it should be of a minimax-optimal order since the common strategy for obtaining optimal rates for the $k$-NN regression estimator is two-fold. First, one should derive a uniform (over $k$) upper bound on the bias term (knowing the smoothness of the regression function), which is a non-decreasing function of $k$. After that, this upper bound is made equal to the variance term, which results in the optimal $k^{\textnormal{b/v}}$. Following this argument, one can conclude that $k^{\textnormal{b/v}} \leq k^*$, which implies $V(k^*) \leq V(k^{\textnormal{b/v}})$. We summarize our findings in the theorem  and corollary below.



\begin{theorem}[Theorem 6.2 in \cite{gyorfi2006distribution}] \label{theorem_gyorfi}
Under the Lipschitz condition (\ref{lipschitz}) on the regression function $f^*$, for any $k \in \{ 1, \ldots, n \}$,
\begin{equation} \label{upper_bound_bias_ineq}
    \mathbb{E}_{\varepsilon} \lVert f^k - f^* \rVert_n^2 \leq C \left( \frac{k}{n} \right)^{2/d} + \frac{\sigma^2}{k},
\end{equation}
where $C$ may depend on $d$ and $L$. Thus, Ineq. (\ref{upper_bound_bias_ineq}) gives $k^{\textnormal{b/v}} = \ceil*{\left( \frac{\sigma^2}{C} \right)^{d/(2+d)} n^{\frac{2}{2 + d}}}$.  
\end{theorem}


\begin{corollary} \label{main_corollary}
    Set $u = \log n$ in Ineq. (\ref{main_ineq}), then, under the $L$-Lipschitz condition (\ref{lipschitz}) on the regression function $f^*$, early stopping rule $k^{\tau}$ from Eq. (\ref{k_tau}) satisfies
    \begin{equation}
        \mathbb{E}_{\varepsilon} \lVert f^{k^{\tau}} - f^* \rVert_n^2 \leq c_u n^{-\frac{2}{2+d}},
    \end{equation}
    where positive constant $c_u$ depends on $d, \sigma$, and $L$; $d > 2$.
\end{corollary}

\begin{proof}
First, by taking the expectation of Ineq. (\ref{main_ineq}), it gives
\begin{align} \label{to_exp_from_proba}
    \begin{split}
    \mathbb{E}_{\varepsilon} \lVert f^{k^{\tau}} - f^* \rVert_n^2 &= \mathbb{E}_{\varepsilon} \left[ \lVert f^{k^{\tau}} - f^* \rVert_n^2 \mathbb{I}\left\{ \lVert f^{k^{\tau}} - f^* \rVert_n^2 \leq 8 V(k^*) + C_1 \frac{\sqrt{\log n}}{\sqrt{n}} + C_2 \frac{\log n}{n} \right\} \right] \\
    &+ \mathbb{E}_{\varepsilon} \left[ \lVert f^{k^{\tau}} - f^* \rVert_n^2 \mathbb{I}\left\{ \lVert f^{k^{\tau}} - f^* \rVert_n^2 > 8 V(k^*) + C_1 \frac{\sqrt{\log n}}{\sqrt{n}} + C_2 \frac{\log n}{n} \right\}  \right].
    \end{split}
\end{align}
After that, due to Lemma \ref{operator_norm_lemma} from the supplementary material, $\lVert I_n - A_k \rVert_2 \leq c$ for any $k \in \{1, \ldots, n\}$,  and $| f^*(x_i) | \leq \mathcal{M}$ for $i \in \{1, \ldots, n\}$ due to the Lipschitz condition (\ref{lipschitz}), which implies that 
\begin{align*}
    \lVert f^{k^{\tau}} - f^* \rVert_n^2 &\leq 2 \lVert (I_n - A_{k^{\tau}})F^* \rVert_n^2 + 2 \lVert A_{k^{\tau}}\varepsilon \rVert_n^2\\
    &\leq 2 \lVert I_n - A_{k^{\tau}} \rVert_2^2 \lVert f^* \rVert_n^2 + 2 \lVert A_{k^{\tau}} \rVert_2^2 \lVert \varepsilon \rVert_n^2\\
    &\leq c_1 + c_2 \lVert \varepsilon \rVert_n^2, 
\end{align*}
where constants $c_1$ and $c_2$ depend only on $\mathcal{M}$ and $d$. Thus,
\begin{equation} \label{aux_inequal}
    \Vert f^{k^{\tau}} - f^* \rVert_n^4 \leq c_1 + c_2 \lVert \varepsilon \rVert_n^4 + c_3 \lVert \varepsilon \rVert_n^2.    
\end{equation}
From Ineq. (\ref{to_exp_from_proba}) and Cauchy-Schwarz inequality, it comes
\begin{align*}
    &\mathbb{E}_{\varepsilon} \lVert f^{k^{\tau}} - f^* \rVert_n^2 \leq 8V(k^*) + C_1 \frac{\sqrt{\log n}}{\sqrt{n}} + C_2 \frac{\log n}{n} \\
    &+ \sqrt{\mathbb{E}_{\varepsilon} \lVert f^{k^{\tau}} - f^* \rVert_n^4} \sqrt{\mathbb{P}_{\varepsilon} \left( \lVert f^{k^{\tau}} - f^* \rVert_n^2 > 8 V(k^*) + C_1 \frac{\sqrt{\log n}}{\sqrt{n}} + C_2 \frac{\log n}{n} \right) }.
\end{align*}
Applying Ineq. (\ref{main_ineq}) and Ineq. (\ref{aux_inequal}), we obtain
\begin{equation*}
    \mathbb{E}_{\varepsilon} \lVert f^{k^{\tau}} - f^* \rVert_n^2 \leq 8 V(k^*) + C_1 \frac{\sqrt{\log n}}{\sqrt{n}} + C_2 \frac{\log n}{n} + \sqrt{c_1 + c_3 \sigma^4 + c_2 \sigma^2}\frac{4}{\sqrt{n}}.  
\end{equation*}
The claim follows from $V(k^*) \leq V(k^{\textnormal{b/v}}) = \sigma^2 / k^{b/v}$, for $k^{\textnormal{b/v}}$ defined in Theorem \ref{theorem_gyorfi}.
\end{proof}


Therefore, the function estimator $f^{k^{\tau}}$ achieves (up to a constant) the minimax bound presented in Eq. (\ref{minimax_lipschitz}), thus non-improvable in general for the class of Lipschitz functions on a bounded domain. 




\section{Empirical comparison to other model selection rules} \label{sec:simulations}

The present section aims at comparing the practical behaviour of our stopping rule $k^{\tau}$ from Eq. (\ref{k_tau}) with other existing and the most used in practice model selection rules. We split the section into three parts: Subsection \ref{description_early_stopping_rules} defines the competitive stopping rules. Subsection \ref{artificial_data} presents experiments on some artificial data sets, while Subsection \ref{real_data} presents experiments on some real data sets. 

\subsection{Description of the model selection rules to compare} \label{description_early_stopping_rules}
In what follows, we will briefly describe five competitive model selection rules. 

\subsection*{Akaike's AIC criterion}
The Akaike's information criterion (\cite{akaike1974new, hastie2009elements}) estimates the risk error by means of a log-likelihood loss function. In the case of $k$-NN regression with Gaussian noise, the maximum likelihood and least-squares are essentially the same things. This gives the AIC criterion as
\begin{equation*}
    R_{\textnormal{AIC}}(f^k) = \frac{1}{n \widehat{\sigma}^2} \left( \lVert Y - A_k Y \rVert^2 + 2 \textnormal{tr}(A_k)\widehat{\sigma}^2 \right), \qquad k = 1, \ldots, n,
\end{equation*}
where $\widehat{\sigma}^2$ is an estimator of $\sigma^2$ obtained from a low-bias model (i.e. with $k = 2$). Using this criterion,
we adjust the training error by a factor proportional to the "degree of freedom". Then, the AIC choice for the optimal $k$ is
\begin{equation} \label{k_aic}
    k_{\textnormal{AIC}} \coloneqq \underset{k = 2, \ldots, n}{\textnormal{argmin}} \left\{ R_{\textnormal{AIC}}(f^k) \right\} - 1.
\end{equation}
Notice that in the mentioned case, AIC criterion is equivalent to Mallow's $C_p$ criterion (\cite{mallows2000some}). AIC criterion has been widely criticized in the literature, especially for the constant $2$ in the definition and/or its asymptotic nature. That is why some authors proposed corrections for the criterion (\cite{schwarz1978estimating, yang1999model}). Nevertheless, AIC, $C_p$, and other related criteria have been proved to satisfy non-asymptotic oracle-type inequalities (see \cite{birge2007minimal} and references therein). Notice that the computational time of the AIC criterion is $\mathcal{O}\left( n^3 \right)$, which is higher than $\mathcal{O}\left( n^2 \left( n - k^{\tau} \right) \right)$ for the minimum discrepancy principle stopping rule $k^{\tau}$ (\ref{k_tau}).

\subsection*{Generalized cross-validation.}
The generalized (GCV) cross-validation strategy (\cite{craven1978smoothing, arlot2010survey}) was introduced in least-squares regression as a rotation-invariant version of the leave-one-out cross-validation procedure. The GCV estimator of the risk error of the linear estimator $A_k Y, \ k = 1, \ldots, n,$ is defined via
\begin{equation*}
    R_{\textnormal{GCV}}(f^k) = \frac{n^{-1}\lVert Y - A_k Y \rVert^2}{(1 - n^{-1}\textnormal{tr}(A_k))^2},    
\end{equation*}
The final model selection rule is
\begin{equation} \label{k_1_out}
    k_{\text{GCV}} \coloneqq \underset{k = 2, \ldots, n}{\textnormal{argmin}} \left\{ R_{\textnormal{GCV}}(f^k) \right\} - 1.
\end{equation}

GCV should be close to the AIC model selection procedure when the sample size $n$ is large. The asymptotic optimality of GCV, meaning that $\lVert f^{k_{\textnormal{GCV}}} - f^* \rVert_n^2 \ / \ \underset{k}{\inf}\lVert f^{k} - f^* \rVert_n^2 \to 1$ in probability, has been proved for the $k$-NN estimator in (\cite{li1987asymptotic}) under mild assumptions. As for the AIC criterion, the computational time of generalized cross-validation is $\mathcal{O}\left( n^3 \right)$ elementary operations.



\subsection*{Hold-out cross-validation stopping rule.}
The Hold-out cross-validation strategy (\cite{geisser1975predictive, arlot2010survey}) is described as follows. The data $\{ x_i, y_i \}_{i=1}^n$ are randomly split into two parts of equal size: the training sample $S_{\textnormal{train}} = \{ x_{\textnormal{train}}, y_{\textnormal{train}} \}$ and the test sample $S_{\textnormal{test}} = \{ x_{\textnormal{test}}, y_{\textnormal{test}} \}$ so that the training and test samples represent a half of the whole data set $\approx n/2$. For each $k = 1, \ldots, n$, one trains the $k$-NN estimator (\ref{estimator}) and evaluates its performance by $R_{\text{HO}}(f^k) = \frac{1}{n}\sum_{i \in S_{\textnormal{test}}}(f^k(x_i)- y_i)^2$, where $f^k(x_i)$ denotes the output of the algorithm trained for $k$ and evaluated at the point $x_i \in x_{\textnormal{test}}$. Then, the Hold-out CV stopping rule is defined via
\begin{equation} \label{k_ho}
    k_{\text{HO}} \coloneqq \underset{k = 2, \ldots, n}{\textnormal{argmin}} \left\{ R_{\text{HO}}(f^k) \right\} - 1. 
\end{equation}
The main inconvenience of this model selection rule is the fact that a part of the data is lost, which increases the risk error. Besides that, the Hold-out strategy is not stable (\cite{arlot2010survey}), which often requires some aggregation of it. The (asymptotic) computational time of the Hold-out strategy is $\mathcal{O}\left( n^3 \right)$ elementary operations.


\subsection*{$V$--fold cross-validation}
$V$--fold cross-validation is certainly the most used cross-validation procedure: the data $\{ (x_i, y_i) \}_{i=1}^n$ are randomly split into $V=5$ equal sized blocks, and at each round (among the $V$ ones), $V-1$ blocks are devoted to training $S_{\text{train}} = (x_{\text{train}}, y_{\text{train}})$, and the remaining one is used for the evaluation $S_{\text{test}} = (x_{\text{test}}, y_{\text{test}})$. The risk error of the $k$-NN estimator is estimated by $R_{\text{VFCV}}(f^k) = \frac{1}{V}\sum_{j=1}^{V} \frac{1}{n/V}\sum_{i \in S_{\text{test}}(j)}\left( f^k(x_i) - y_i \right)^2$, where $f^k(x_i)$ denotes the output of the algorithm trained for $k$ and evaluated at the point $x_i \in S_{\text{test}}(j)$, thus 
\begin{equation} \label{k_vfcv}
    k_{\text{VFCV}} \coloneqq \underset{k = 2, \ldots, n}{\textnormal{argmin}} \left\{ R_{\text{VFCV}}(f^k) \right\} - 1.
\end{equation}
$V$--fold cross-validation is a more computationally tractable solution than other splitting-based model selection methods, such as the leave-one-out or leave-$p$-out (\cite{hastie2009elements, arlot2010survey}). Usually, the optimal $V$ is equal to $5$ or $10$ due to the fact that the statistical error does not increase a lot for larger values of $V$ whereas averaging over more than $10$ folds becomes infeasible. To the best of our knowledge, there are no theoretical results for the $V$--fold cross validation model selection strategy with the $k$-NN regression estimator.  


\subsection*{Theoretical bias-variance trade-off stopping rule}

The fourth stopping rule is the one introduced in Eq. (\ref{k_star}). This stopping rule is the classical bias-variance trade-off stopping rule that provides minimax-optimal rates (see the monographs \cite{wasserman2006all, tsybakov2008introduction}):
\begin{equation}
    k^* = \inf \{ k \in \{ 1, \ldots, n \} \mid B^2(k) \geq V(k) \}. 
\end{equation}
The stopping rule $k^*$ is introduced for comparison purposes only because it \textit{cannot be computed} in practice. One can say that this stopping rule is minimax-optimal if $f^*$ belongs, for instance, to the class of Lipschitz functions on a bounded domain (\ref{lipschitz}). Therefore, it could serve as a (lower bound) reference in the present simulated experiments with artificial data.





\subsection{Artificial data} \label{artificial_data}

First, the goal is to perform some simulated experiments (a comparison of mentioned stopping rules) on artificial data.

\vspace{0.2cm}

\textbf{Description of the simulation design}

The data in this case is generated according to the regression model $y_j = f^*(x_j) + \varepsilon_j$, where $\varepsilon_j \overset{i.i.d.}{\sim} \mathcal{N}(0, \sigma^2)$ (Gaussian),  $j = 1, \ldots, n$. We choose the uniform covariates $x_j \overset{i.i.d.}{\sim} \mathbb{U}[0, 1]^3$ (uniform), $j = 1, \ldots, n,$ and $\sigma = 0.15$. Consider two regression functions with different smoothness: a "smooth" $f_1^*(x) = 1.5 \cdot \left[ \lVert x - 0.5 \rVert / \sqrt{3} - 0.5 \right]$ and a "sinus" $f_2^*(x) = 1.5 \cdot \sin (\lVert x \rVert / \sqrt{3})$ for any $x \in [0, 1]^3$. Notice that both functions belong to the class of Lipschitz functions (\ref{lipschitz}) on $[0, 1]^3$. The sample size $n$ varies from $50$ to $250$.

The $k$-NN algorithm (\ref{estimator}) is trained first for $k = \floor*{\sqrt{n}}$, after that we decrease the value of $k$ until $k = 1$ such that at each step of the iteration procedure we increase the variance of the $k$-NN estimator $V(k) = \sigma^2/k$ (cf. Fig. \ref{fig:bvr}). In other words, the model becomes more complex successively due to the increase of its "degree of freedom" measured by $\textnormal{tr}(A_k) = n/k$. If the condition in Eq. (\ref{k_tau}) is satisfied, the learning process is stopped, and it outputs the stopping rule $k^{\tau}$. 

The performance of the stopping rules is measured in terms of the empirical $L_2(\mathbb{P}_n)$-norm $\lVert f^k - f^* \rVert_n^2$ averaged over $N = 1000$ repetitions (over the noise $\{ \varepsilon_j \}_{j=1}^n$).  

For our simulations, we use a consistent (low-bias) estimator of $\sigma^2$ described in Section \ref{real_data} (see Eq. (\ref{var_est})).

\begin{figure}
\begin{subfigure}[b]{0.45\textwidth}
\begin{tikzpicture}[scale=0.7]
\tikzstyle{every node}=[font=\normalsize]
\begin{axis}[
    title={$k$-NN regressor, $\sigma = 0.15$, smooth r.f.},
    xlabel={Sample size},
    ylabel={Average loss},
    xmin=50, xmax=250,
    ymin=0.004, ymax=0.0115,
    ytick={0.005,0.007,0.009},
    xtick={50, 100, 150, 200, 250},
    legend pos=north east,
    ymajorgrids=true,
    xmajorgrids=true,
    grid style=dashed,
]
\addplot[
    color=blue,
    mark=star,
    ]
    coordinates {
    (50,0.00959589)(80,0.00723882)(100,0.0064927)(160,0.00538029)(200,0.00527075)(250,0.00440348)
    };
    \addlegendentry{$k^{\tau}$}
   
\addplot[
    color=red,
    mark=+,
    ]
    coordinates {
    (50, 0.01117475)(80, 0.00820779)(100, 0.00692115)(160, 0.00578596)(200, 0.00556723)(250,0.00480027)
    };
    \addlegendentry{$k_{\textnormal{HO}}$}
\addplot[
    color=black,
    mark=x,
    ]
    coordinates {
    (50, 0.00960814)(80, 0.00708104)(100, 0.0064098)(160, 0.00534147)(200, 0.005275)(250,0.00434834)
    };
    \addlegendentry{$k_{\textnormal{GCV}}$}
\addplot[
    color=green,
    mark=star,
    ]
    coordinates {
    (50, 0.00934164)(80, 0.00708591)(100, 0.0062164)(160, 0.00527903)(200, 0.00512319)(250,0.00420769)
    };
    \addlegendentry{$k^{*}$}
\end{axis}
\end{tikzpicture}
\caption{}
\label{smooth}
\end{subfigure}
\qquad
\begin{subfigure}[b]{0.45\textwidth}
\begin{tikzpicture}[scale=0.7]
\tikzstyle{every node}=[font=\normalsize]
\begin{axis}[
    title={$k$-NN regressor, $\sigma = 0.15$, sinus r.f.},
    xlabel={Sample size},
    ylabel={Average loss},
    xmin=50, xmax=250,
    ymin=0.0025, ymax=0.0095,
    xtick={50,100,150,200,250},
    ytick={0.003,0.005, 0.008},
    legend pos=north east,
    ymajorgrids=true,
    xmajorgrids=true,
    grid style=dashed,
]
\addplot[
    color=blue,
    mark=star,
    ]
    coordinates {
    (50,0.00856475)(80,0.00527245)(100,0.00520583)(160,0.00414462)(200,0.00360486)(250,0.0029507)
    };
    \addlegendentry{$k^{\tau}$}
   
\addplot[
    color=red,
    mark=+,
    ]
    coordinates {
    (50, 0.00932031)(80, 0.00640974)(100, 0.00562799)(160, 0.00452315)(200, 0.00401955)(250,0.00365806)
    };
    \addlegendentry{$k_{\textnormal{HO}}$}
\addplot[
    color=black,
    mark=x,
    ]
    coordinates {
    (50, 0.00863128)(80, 0.00534549)(100, 0.00504782)(160, 0.00411817)(200, 0.00362582)(250,0.00289884)
    };
    \addlegendentry{$k_{\textnormal{GCV}}$}
\addplot[
    color=green,
    mark=star,
    ]
    coordinates {
    (50, 0.00861029)(80, 0.00477626)(100, 0.00499522)(160, 0.00389686)(200, 0.00349439)(250,0.00275263)
    };
    \addlegendentry{$k^{*}$}
\end{axis}
\end{tikzpicture}
\caption{}
\label{sinus}
\end{subfigure}
\caption{$k$-NN estimator (\ref{estimator}) with two noised regression functions: smooth $f_1^*(x) = 1.5 \cdot \left[ \lVert x - 0.5 \rVert / \sqrt{3} - 0.5 \right]$ for panel (a) and "sinus" $f_2^*(x) = 1.5 \cdot \textnormal{sin}(\lVert x \rVert / \sqrt{3})$ for panel (b), with uniform covariates $x_j \overset{i.i.d.}{\sim} \mathbb{U}[0, 1]^3$. Each curve corresponds to the $L_2(\mathbb{P}_n)$ squared norm error for the stopping rules (\ref{k_tau}), (\ref{k_star}), (\ref{k_ho}), (\ref{k_1_out}), averaged over $1000$ independent trials, versus the sample size $n = \{50, 80, 100, 160, 200, 250 \}$.}
 \label{fig:comp_knn_art}
\end{figure}




\textbf{Results of the simulation experiments.}

Figure \ref{fig:comp_knn_art} displays the resulting (averaged over $1000$ repetitions) $L_2(\mathbb{P}_n)$-error of $k^{\tau}$ (\ref{k_tau}), $k^*$ (\ref{k_star}), $k_{\text{HO}}$ (\ref{k_ho}), and $k_{\text{GCV}}$ (\ref{k_1_out}) versus the sample size $n$. In particular, Figure \ref{smooth} shows the results for the "smooth" regression function, whereas Figure \ref{sinus} provides the results for the "sinus" regression function. 

First, from all the graphs, (almost) all curves do not increase as the sample size $n$ grows. One can notice that in both graphs, the prediction error of the Hold-out strategy is the worst among the model selection criteria.

In more detail, Figure \ref{smooth} indicates that the best performance is achieved by $k^*$ (non-computable in practice bias-variance trade-off).
Besides that, the minimum discrepancy principle rule $k^{\tau}$ is uniformly better than $k_{\textnormal{HO}}$ and has the same performance as the one of $k_{\textnormal{GCV}}$. Moreover, the gap between $k^{\tau}$ and $k^*$ is getting smaller for the sample sizes $n \geq 160$. This behavior supports the theoretical part of the present paper because $k^{\tau}$ should serve as an estimator of $k^*$. Since $k^*$ is the well-known bias-variance trade-off, the minimum discrepancy principle stopping rule seems to be a meaningful model selection method.

In Figure \ref{sinus}, the best performance is achieved again by $k^*$ -- a non-computable in practice stopping rule. As for the data-driven model selection methods, the stopping rules $k^{\tau}$ and $k_{\text{GCV}}$ (an asymptotically optimal model selection strategy) perform almost equivalently. 


\subsection{Real data} \label{real_data}

Here, we tested the performance (prediction error and runtime) of the early stopping rule $k^{\tau}$ (\ref{k_tau}) for choosing the hyperparameter in the $k$-NN estimator on four different data sets mostly taken from the UCI repository (\cite{Dua:2019}). 


\textbf{Data sets description}


The housing data set "Boston Housing Prices" concerns the task of predicting housing values in areas of Boston (USA), the input points are $13$-dimensional. 

The "Diabetes" data set consists of 10 columns that measure different patient's characteristics (age, sex, body mass index, etc), the output is a quantitative measure of disease progression one year after the baseline. 

The "Power Plants" data set contains 9568 data points collected from a Combined Cycle Power Plant over 6 years (2006-2011), when the plant was set to work with the full load. 

"California Houses Prices" data set (\cite{pace1997sparse}) contains information from the 1990's California census. The input variables are "total bedrooms", "total rooms", etc. The output variable is the median house value for households within a block (measured in US Dollars). 

Notice that for "California Houses Prices" and "Power Plants" data sets we take the first 3000 data points in order to speed up the calculations. 

\vspace{0.2cm}

\textbf{Description of the simulation design}

Assume that we are given one of the data sets described above. Let us rescale each variable of this data set $\widetilde{x} \in \mathbb{R}^n$ such that all the components $\widetilde{x}_i, \ i = 1, \ldots, n$, belong to $[0, 1]$:
\begin{equation*}
    \widetilde{x}_i = \frac{\widetilde{x}_i - \min(\widetilde{x})}{\max(\widetilde{x}) - \min(\widetilde{x})}, \ i = 1, \ldots, n,
\end{equation*}
where $\min(\widetilde{x})$ and $\max(\widetilde{x})$ denote the minimum and the maximum component of the vector $\widetilde{x}$.

After that, we split the data set into two parts: one is denoted $S_{\textnormal{train}} = \{ x_{\textnormal{train}}, y_{\textnormal{train}} \}$ (70 \% of the whole data) and is made for training and model selection (model selection rules $k^{\tau}$, $k_{\text{GCV}}$, $k_{\text{5FCV}}$, and $k_{\text{AIC}}$), the other one (30 \% of the whole data) is denoted $S_{\textnormal{test}} = \{ x_{\textnormal{test}}, y_{\textnormal{test}} \}$ and is made for making prediction on it. We denote $n_{\textnormal{train}}$ and $n_{\textnormal{test}}$ as the sample sizes of $S_{\textnormal{train}}$ and $S_{\textnormal{test}}$, respectively. Then, our experiments' design is divided into four parts.

\vspace{0.15cm}

At the beginning, we create a grid of sub-sample size for each data set:
\begin{equation} \label{subsamples}
    n_s \in \Big\{\floor*{n_{\textnormal{train}}/5}, \floor*{n_{\textnormal{train}}/4}, \floor*{n_{\textnormal{train}}/3}, \\ \floor*{n_{\textnormal{train}}/2}, n_{\textnormal{train}}\Big\},
\end{equation}
and a grid of the maximum number of neighbors $k_{\textnormal{max}} = 3 \floor*{\log(n_s)}$, where $n_{\textnormal{train}} = \ceil*{0.7n}$ and $n$ is the sample size of the whole data. 

Further, for each data set and sub-sample size from Eq. (\ref{subsamples}), we estimate the noise variance $\sigma^2$ from the regression model (\ref{nonparametric_regression_model}). There is large amount of work on the estimation of $\sigma^2$ in nonparametric regression (\cite{rice1984bandwidth, hall1990variance}). In our simulated experiments, we take the estimator from \cite[Eq. (5.86)]{wasserman2006all}, which is a consistent estimator of $\sigma^2$ under an assumption that $f^*$ is "sufficiently smooth". It satisfies our simulation experiments' purposes.
\begin{equation} \label{var_est}
    \widehat{\sigma}^2 \coloneqq \frac{\lVert (I_{n_s} - A_k)y_{s} \rVert^2}{n_s (1 - 1/k)} \ \ \ \textnormal{ with } \ \ \ k = 2, 
\end{equation}
where $y_s$ corresponds to the vector of responses from the chosen sub-samples.
After that, we compute our stopping rule $k^{\tau}$ and other model selection strategies from Section \ref{description_early_stopping_rules}. To do that, for each data set and each integer $n_s$ from Eq. (\ref{subsamples}), we randomly sample $n_s$ data points from $S_{\textnormal{train}}$, compute the $k$-NN estimator (\ref{estimator}) and the empirical risk (\ref{empirical_risk}) for $k_{\textnormal{max}}$, and at each step of the iteration process we reduce the value of $k$ by one. Remark that one does not have to calculate the neighborhood matrix $A_k$ for each $k \in \{1, \ldots, k_{\textnormal{max}}\}$, since it is sufficient to do only for $k_{\textnormal{max}}$ (cf. Eq. (\ref{iteration_for_matrix})). This process is repeated until the empirical risk crosses the threshold $\widehat{\sigma}^2$. Fig. \ref{fig:process_learning} provides two illustrations of the minimum discrepancy strategy $k^{\tau}$ for two data sets: "Diabetes" and "Boston Housing Prices".

After that, the AIC criterion (\ref{k_aic}), $5$--fold cross-validation (\ref{k_vfcv}), and the generalized cross-validation $k_{\textnormal{GCV}}$ are calculated. Let us describe how we do that in two steps. We start by defining a grid of values for $k: \{ 1, 2, \ldots, k_{\textnormal{max}} \}$. Further, one should compute $k_{\textnormal{AIC}}$, $k_{\textnormal{GCV}}$, and $k_{\textnormal{5FCV}}$ from Eq. (\ref{k_aic}), Eq. (\ref{k_1_out}), and Eq. (\ref{k_vfcv}) over the mentioned grid.

In the final part, given $k^{\tau}$, $k_{\text{AIC}}$, $k_{\text{5FCV}}$, and $k_{\text{GCV}}$, the goal is to make a prediction on the test data set $x_{\textnormal{test}}$. This can be done as follows. Assume that $x_0 \in x_{\textnormal{test}}$, then the prediction of the $k$-NN estimator on this point can be defined as
\begin{equation*}
    f^k(x_0) = a_k(x_0)^\top y_{s},
\end{equation*}
where $a_k(x_0) = [a_k(x_0, x_1), \ldots, a_k(x_0, x_{n_s})]^\top$ and $x_{s} = [x_1^{\top}, \ldots, x_{n_{s}}^\top]^{\top}$, with $a_k(x_0, x_i) = 1/k$ if $x_i, \ i \in \{1, \ldots, n_{\textnormal{train}} \}$, is a nearest neighbor of $x_0$, otherwise $0$. Further, one can choose $k$ to be equal $k^{\tau}$, $k_{\text{AIC}}$, $k_{\text{5FCV}}$ or $k_{\text{GCV}}$ that are already computed. Combining all the steps together, one is able to calculate the least-squares prediction error $\lVert f^k - y_{\textnormal{test}} \rVert$.

For each sub-sample size $n_s$ from Eq. (\ref{subsamples}) and data set, the procedure has to be performed $25$ times (via new sub-samples from the data set). 

\vspace{0.2cm}

\textbf{Results of the simulation experiments.}

Figures \ref{small_datasets_results} and \ref{large_datasets_results} display the averaged (over 25 repetitions) runtime (in seconds) and the prediction error of the model selection rules $k^{\tau}$ (\ref{k_tau}), $k_{\textnormal{AIC}}$ (\ref{k_aic}), 5-fold cross-validation (\ref{k_vfcv}), and generalized cross-validation (\ref{k_1_out}) for "Boston Housing Prices", "Diabetes" (in Figure \ref{small_datasets_results}), and "California Houses Prices", "Power Plants" data sets (in Figure \ref{large_datasets_results}).

Figures \ref{runtime_boston}, \ref{runtime_diabetes} indicate that the minimum discrepancy principle rule $k^{\tau}$ has the smallest runtime among the model selection criteria. At the same time, Figure \ref{prediction_error_boston} shows that the prediction error of $k^{\tau}$ has better performance than that of $k_{\textnormal{AIC}}$ and $k_{\textnormal{5FCV}}$, and is similar to that of $k_{\textnormal{GCV}}$. Figure \ref{pred_error_diabetes} indicates that the performance of the minimum discrepancy stopping rule $k^{\tau}$ is better than that of $k_{\textnormal{AIC}}$, $k_{\textnormal{GCV}}$, and $k_{\textnormal{5FCV}}$ for the sub-sample size $n_s \geq 100$.

Let us turn to the results for the "California Houses Prices" and "Power Plants" data sets. Figures \ref{prediction_error_california}, \ref{prediction_error_power} display the prediction performance of the model selection rules: for the "California Houses Prices" data set, the prediction performance of $k^{\tau}$ is comparable to the performance of $k_{\textnormal{GCV}}$ and is uniformly better than the $5$FCV rule $k_{\textnormal{5FCV}}$ and $k_{\textnormal{AIC}}$; for the "Power Plants" data set, the prediction error of the minimum discrepancy principle is similar to that of $k_{\textnormal{AIC}}, k_{\textnormal{GCV}}$, and $k_{\textnormal{5FCV}}$ for the sub-sample sizes $n_s \leq 1000$, and is little worse for $n_s = 2100$. This deterioration of the performance could be due to the estimation of the variance $\sigma^2$ from Eq. (\ref{var_est}). Figures \ref{runtime_california}, \ref{runtime_power} show the runtime of the stopping rules: one can conclude that the computational time of the minimum discrepancy rule $k^{\tau}$ is less than the computational time of the generalized cross-validation, AIC criterion, and $5$ fold cross-validation. 

The overall conclusion from the simulation experiments is that the prediction error of the MDP stopping rule $k^{\tau}$ is often better than for standard model selection strategies, such as the AIC criterion or $5$--fold cross-validation, while its computational time is lower.   

\begin{figure}
\begin{subfigure}[b]{0.45\textwidth}
\begin{tikzpicture}[scale=0.7]
\tikzstyle{every node}=[font=\normalsize]
\begin{axis}[
    title={Diabetes data (sub-sample)},
    xlabel={Number of neighbors},
    ylabel={Value},
    xmin=0, xmax=10,
    ymin=0, ymax=6114.66446281,
    xtick={1,3,5,7, 10},
    ytick={1500, 4500},
    legend pos=south east,
    ymajorgrids=true,
    xmajorgrids=true,
    grid style=dashed,
]

\addplot[
    color=blue,
    mark=square,
    ]
    coordinates {
    (1,0)(2,0)(3,3664.18333333)(4,4731.1125)(5,5358.508)(6,5680.28472222)(7,5452.15408163)(8,5851.84765625)(9,5631.3308642)(10,6114.66446281)
    };
    \addlegendentry{Empirical risk}
   
\addplot[
    color=red,
    mark=+,
    ]
    coordinates {
    (1,4215.1)(2,4215.1)(3,4215.1)(4,4215.1)(5,4215.1)(6,4215.1)(7,4215.1)(8,4215.1)(9,4215.1)(10,4215.1)
    };
    \addlegendentry{Threshold}
\addplot +[mark=none, color=orange] coordinates {(3, 0)(3, 6114.66446281)};
\addlegendentry{$k^{\tau}$}
\end{axis}
\end{tikzpicture}
\caption{}
\label{a}
\end{subfigure}
\qquad
\begin{subfigure}[b]{0.45\textwidth}
\begin{tikzpicture}[scale=0.7]
\tikzstyle{every node}=[font=\normalsize]
\begin{axis}[
    title={Boston Housing data (sub-sample)},
    xlabel={Number of neighbors},
    ylabel={Value},
    xmin=0, xmax=10,
    ymin=0, ymax=68.4,
    xtick={1,3,5,10},
    ytick={20, 45, 65},
    legend pos=south east,
    ymajorgrids=true,
    xmajorgrids=true,
    grid style=dashed,
]

\addplot[
    color=blue,
    mark=square,
    ]
    coordinates {
    (1,0)(2,0)(3,0)(4,49.69747283)(5,63.26987826)(6,60.0285628)(7,60.71023957)(8,67.93981658)(9,67.28213634)(10,68.37776957)
    };
    \addlegendentry{Empirical risk}
   
\addplot[
    color=red,
    mark=+,
    ]
    coordinates {
    (1, 53.79369565)(2, 53.79369565)(3, 53.79369565)(4, 53.79369565)(5, 53.79369565)(6,53.79369565)(7,53.79369565)(8,53.79369565)(9,53.79369565)(10,53.79369565)
    };
    \addlegendentry{Threshold}
\addplot +[mark=none, color=orange] coordinates {(4, 0)(4, 6114.66446281)};
\addlegendentry{$k^{\tau}$}
\end{axis}
\end{tikzpicture}
\caption{}
\label{b}
\end{subfigure}
\caption{Stopping the learning process based on the rule (\ref{k_tau}) applied to two data sets: a) "Diabetes" and b) "Boston Housing Prices". "Threshold" horizontal line corresponds to the estimated variance from Eq. (\ref{var_est}).}
\label{fig:process_learning}
\end{figure}




\begin{figure}
\begin{subfigure}[b]{0.45\textwidth}
\begin{tikzpicture}[scale=0.65]
\tikzstyle{every node}=[font=\normalsize]
\begin{axis}[
    title={Boston H.P., runtime},
    xlabel={Sub-sample size},
    ylabel={Average time (sec)},
    xmin=70, xmax=354,
    ymin=0.003, ymax=0.11,
    xtick={70,100,140,200,300,400},
    ytick={0.01,0.1},
    legend pos=north west,
    ymajorgrids=true,
    xmajorgrids=true,
    grid style=dashed,
]

\addplot[
    color=blue,
    mark=square,
    ]
    coordinates {
    (70,0.00855612)(88,0.01127354)(118,0.01240211)(177,0.01705173)(354,0.05850043)
    };
    \addlegendentry{$k^{\tau}$}
   
\addplot[
    color=red,
    mark=+,
    ]
    coordinates {
    (70, 0.01082651)(88, 0.01258839)(118, 0.01454256)(177, 0.01753586)(354, 0.06210802)
    };
    \addlegendentry{$k_{\textnormal{AIC}}$}
\addplot[
    color=black,
    mark=x,
    ]
    coordinates {
    (70, 0.00830978)(88, 0.02092349)(118, 0.01145007)(177, 0.01990505)(354, 0.06201828)
    };
    \addlegendentry{$k_{\textnormal{GCV}}$}
\addplot[
    color=violet,
    mark=*,
    ]
    coordinates {
    (70, 0.02412378)(88, 0.0344868)(118, 0.0204067)(177, 0.03076374)(354, 0.09219963)
    };
    \addlegendentry{$k_{\textnormal{5FCV}}$}

\end{axis}
\end{tikzpicture}
\caption{}
\label{runtime_boston}
\end{subfigure}
\qquad
\begin{subfigure}[b]{0.45\textwidth}
\begin{tikzpicture}[scale=0.65]
\tikzstyle{every node}=[font=\normalsize]
\begin{axis}[
    title={Boston H.P., prediction error},
    xlabel={Sub-sample size},
    ylabel={Average loss},
    xmin=70, xmax=354,
    ymin=50, ymax=85,
    xtick={70,100,140,200,300,400},
    ytick={60, 70, 80},
    legend pos=north east,
    ymajorgrids=true,
    xmajorgrids=true,
    grid style=dashed,
]

\addplot[
    color=blue,
    mark=square,
    ]
    coordinates {
    (70,74.2853426)(88,73.28846512)(118,67.86918275)(177,62.26424591)(354,50.59923418)
    };
    \addlegendentry{$k^{\tau}$}
   
\addplot[
    color=red,
    mark=+,
    ]
    coordinates {
    (70, 75.8850821)(88, 73.72479664)(118, 67.90502476)(177, 62.45029884)(354, 53.17254931)
    };
    \addlegendentry{$k_{\textnormal{AIC}}$}
\addplot[
    color=black,
    mark=x,
    ]
    coordinates {
    (70, 75.6347345)(88, 72.5862216)(118, 68.21207572)(177, 63.64873349)(354, 50.59923418)
    };
    \addlegendentry{$k_{\textnormal{GCV}}$}
\addplot[
    color=violet,
    mark=*,
    ]
    coordinates {
    (70, 75.20374274)(88, 72.83081082)(118, 68.20018387)(177, 63.40840655)(354, 52.9666841)
    };
    \addlegendentry{$k_{\textnormal{5FCV}}$}
\end{axis}
\end{tikzpicture}
\caption{}
\label{prediction_error_boston}
\end{subfigure}
\\
\begin{subfigure}[b]{0.45\textwidth}
\begin{tikzpicture}[scale=0.65]
\tikzstyle{every node}=[font=\normalsize]
\begin{axis}[
    title={Diabetes, runtime},
    xlabel={Sub-sample size},
    ylabel={Average time (sec)},
    xmin=60, xmax=309,
    ymin=0.003, ymax=0.09,
    xtick={61,100,140,200,300},
    ytick={0.01,0.05},
    legend pos=north west,
    ymajorgrids=true,
    xmajorgrids=true,
    grid style=dashed,
]

\addplot[
    color=blue,
    mark=square,
    ]
    coordinates {
    (61,0.00528295)(77,0.006069)(103,0.00993328)(154,0.01101583)(309,0.04544517)
    };
    \addlegendentry{$k^{\tau}$}
   
\addplot[
    color=red,
    mark=+,
    ]
    coordinates {
    (61, 0.0069945)(77, 0.00938141)(103, 0.00800133)(154, 0.01413276)(309, 0.05715765)
    };
    \addlegendentry{$k_{\textnormal{AIC}}$}
\addplot[
    color=black,
    mark=x,
    ]
    coordinates {
    (61, 0.00631852)(77, 0.0088942)(103, 0.00800728)(154, 0.01413908)(309, 0.0594244)
    };
    \addlegendentry{$k_{\textnormal{GCV}}$}
\addplot[
    color=violet,
    mark=*,
    ]
    coordinates {
    (61, 0.01438442)(77, 0.01810751)(103, 0.01518609)(154, 0.02330532)(309, 0.09792542)
    };
    \addlegendentry{$k_{\textnormal{5FCV}}$}
\end{axis}
\end{tikzpicture}
\caption{}
\label{runtime_diabetes}
\end{subfigure}
\qquad
\begin{subfigure}[b]{0.45\textwidth}
\begin{tikzpicture}[scale=0.65]
\tikzstyle{every node}=[font=\normalsize]
\begin{axis}[
    title={Diabetes, prediction error},
    xlabel={Sub-sample size},
    ylabel={Average loss},
    xmin=61, xmax=309,
    ymin=630, ymax=740,
    xtick={61,100,140,200,300},
    ytick={650, 700, 725},
    legend pos=north east,
    ymajorgrids=true,
    xmajorgrids=true,
    grid style=dashed,
]

\addplot[
    color=blue,
    mark=square,
    ]
    coordinates {
    (61,735.32410221)(77,698.07531817)(103,675.96322628)(154,653.35361124)(309,637.9835621)
    };
    \addlegendentry{$k^{\tau}$}
   
\addplot[
    color=red,
    mark=+,
    ]
    coordinates {
    (61, 729.27694327)(77, 694.41295062)(103, 676.86668009)(154, 658.35153366)(309, 649.31011705)
    };
    \addlegendentry{$k_{\textnormal{AIC}}$}
\addplot[
    color=black,
    mark=x,
    ]
    coordinates {
    (61, 715.16104155)(77, 694.23348373)(103, 679.64844807)(154, 656.60010857)(309, 649.31011705)
    };
    \addlegendentry{$k_{\textnormal{GCV}}$}
\addplot[
    color=violet,
    mark=*,
    ]
    coordinates {
    (61, 730.96038738)(77, 716.52927512)(103, 679.82796817)(154, 674.75239488)(309, 645.99680041)
    };
    \addlegendentry{$k_{\textnormal{5FCV}}$}
\end{axis}
\end{tikzpicture}
\caption{}
\label{pred_error_diabetes}
\end{subfigure}
\caption{Runtime (in seconds) and $L_2(\mathbb{P}_n)$ prediction error versus sub-sample size for different model selection methods: MD principle (\ref{k_tau}), AIC (\ref{k_aic}), GCV (\ref{k_1_out}), and $5$--fold cross-validation (\ref{k_vfcv}), tested on the "Boston Housing Prices" and "Diabetes" data sets. In all cases, each point corresponds to the average of $25$ trials. (a), (c) Runtime verus the sub-sample size $n \in \{70, 88, 118, 177, 354 \}$. (b), (d) Least-squares prediction error $\lVert f^k - y_{\textnormal{test}} \rVert$ versus the sub-sample size $n \in \{70, 88, 118, 177, 354 \}$.}
\label{small_datasets_results}
\end{figure}


\begin{figure}
\begin{subfigure}[b]{0.45\textwidth}
\begin{tikzpicture}[scale=0.65]
\tikzstyle{every node}=[font=\normalsize]
\begin{axis}[
    title={California H.P., runtime},
    xlabel={Sub-sample size},
    ylabel={Average time (sec)},
    xmin=420, xmax=2100,
    ymin=0.02, ymax=1.65,
    xtick={500,700,1000,2000},
    ytick={0.1,0.5, 1},
    legend pos=north west,
    ymajorgrids=true,
    xmajorgrids=true,
    grid style=dashed,
]

\addplot[
    color=blue,
    mark=square,
    ]
    coordinates {
    (420,0.06766566)(525,0.1230032)(700,0.14770074)(1050,0.40171222)(2100,1.35318375)
    };
    \addlegendentry{$k^{\tau}$}
   
\addplot[
    color=red,
    mark=+,
    ]
    coordinates {
    (420,0.0669849)(525,0.11980861)(700,0.21817273)(1050,0.43286912)(2100,1.48403368)
    };
    \addlegendentry{$k_{\textnormal{AIC}}$}
   
\addplot[
    color=black,
    mark=x,
    ]
    coordinates {
    (420,0.06887476)(525,0.12885856)(700,0.17433068)(1050,0.42578732)(2100,1.45269584)
    };
    \addlegendentry{$k_{\textnormal{GCV}}$}
   
\addplot[
    color=violet,
    mark=*,
    ]
    coordinates {
    (420,0.11478419)(525,0.16037658)(700,0.26191102)(1050,0.63462547)(2100,1.92028911)
    };
    \addlegendentry{$k_{\textnormal{5FCV}}$}
\end{axis}
\end{tikzpicture}
\caption{}
\label{runtime_california}
\end{subfigure}
\qquad
\begin{subfigure}[b]{0.45\textwidth}
\begin{tikzpicture}[scale=0.65]
\tikzstyle{every node}=[font=\normalsize]
\begin{axis}[
    title={California H.P., prediction error},
    xlabel={Sub-sample size},
    ylabel={Average loss},
    xmin=420, xmax=2100,
    ymin=11.2, ymax=13.5,
    xtick={500,700,1000,2000},
    ytick={12 ,13 , 13.5},
    legend pos=north east,
    ymajorgrids=true,
    xmajorgrids=true,
    grid style=dashed,
]

\addplot[
    color=blue,
    mark=square,
    ]
    coordinates {
    (420,13.19666715)(525,13.02057484)(700,12.61806828)(1050,12.12776051)(2100,11.43766597)
    };
    \addlegendentry{$k^{\tau}$}
   
\addplot[
    color=red,
    mark=+,
    ]
    coordinates {
    (420,13.17638071)(525,13.1502115)(700,12.76413759)(1050,12.18476347)(2100,11.43766597)
    };
    \addlegendentry{$k_{\textnormal{AIC}}$}
\addplot[
    color=black,
    mark=x,
    ]
    coordinates {
    (420,13.15911354)(525,12.98821598)(700,12.54897291)(1050,12.10185179)(2100,11.46773943)
    };
    \addlegendentry{$k_{\textnormal{GCV}}$}
\addplot[
    color=violet,
    mark=*,
    ]
    coordinates {
    (420,13.3404094)(525,13.27833772)(700,12.66432926)(1050,12.27531468)(2100,11.57245704)
    };
    \addlegendentry{$k_{\textnormal{5FCV}}$}
\end{axis}
\end{tikzpicture}
\caption{}
\label{prediction_error_california}
\end{subfigure}
\\
\begin{subfigure}[b]{0.45\textwidth}
\begin{tikzpicture}[scale=0.65]
\tikzstyle{every node}=[font=\normalsize]
\begin{axis}[
    title={Power Plants, runtime},
    xlabel={Sub-sample size},
    ylabel={Average time (sec)},
    xmin=350, xmax=2100,
    ymin=0.03, ymax=1.8,
    xtick={420,700,1000,2000},
    ytick={0.1, 1, 2},
    legend pos=north west,
    ymajorgrids=true,
    xmajorgrids=true,
    grid style=dashed,
]

\addplot[
    color=blue,
    mark=square,
    ]
    coordinates {
    (420,0.05278417)(525,0.09394525)(700,0.14946687)(1050,0.32476645)(2100,1.31787278)
    };
    \addlegendentry{$k^{\tau}$}
   
\addplot[
    color=red,
    mark=+,
    ]
    coordinates {
    (420, 0.05410801)(525, 0.11433973)(700, 0.18229741)(1050, 0.38104916)(2100, 1.55990878)
    };
    \addlegendentry{$k_{\textnormal{AIC}}$}
\addplot[
    color=black,
    mark=x,
    ]
    coordinates {
    (420, 0.05847129)(525, 0.09902878)(700, 0.16085543)(1050, 0.36054396)(2100, 1.41506519)
    };
    \addlegendentry{$k_{\textnormal{GCV}}$}
\addplot[
    color=violet,
    mark=*,
    ]
    coordinates {
    (420, 0.08910254)(525, 0.17428734)(700, 0.21586569)(1050, 0.49328215)(2100, 1.84269897)
    };
    \addlegendentry{$k_{\textnormal{5FCV}}$}
\end{axis}
\end{tikzpicture}
\caption{}
\label{runtime_power}
\end{subfigure}
\qquad
\begin{subfigure}[b]{0.45\textwidth}
\begin{tikzpicture}[scale=0.65]
\tikzstyle{every node}=[font=\normalsize]
\begin{axis}[
    title={Power Plants, prediction error},
    xlabel={Sub-sample size},
    ylabel={Average loss},
    xmin=350, xmax=2100,
    ymin=115, ymax=136,
    xtick={420,700,1000,2000},
    ytick={120, 125, 130},
    legend pos=north east,
    ymajorgrids=true,
    xmajorgrids=true,
    grid style=dashed,
]

\addplot[
    color=blue,
    mark=square,
    ]
    coordinates {
    (420,134.9333523)(525,131.49142363)(700,128.50263383)(1050,124.02441083)(2100,117.00074146)
    };
    \addlegendentry{$k^{\tau}$}
   
\addplot[
    color=red,
    mark=+,
    ]
    coordinates {
    (420, 136.59589921)(525, 131.38593144)(700, 129.14665047)(1050, 124.20131175)(2100, 115.65787105)
    };
    \addlegendentry{$k_{\textnormal{AIC}}$}
\addplot[
    color=black,
    mark=x,
    ]
    coordinates {
    (420, 135.36142759)(525, 131.05128992)(700, 128.72729084)(1050, 124.05305432)(2100, 115.65787105)
    };
    \addlegendentry{$k_{\textnormal{GCV}}$}
\addplot[
    color=violet,
    mark=*,
    ]
    coordinates {
    (420, 135.40563633)(525, 131.06600903)(700, 128.2952843)(1050, 124.04019418)(2100, 116.08168959)
    };
    \addlegendentry{$k_{\textnormal{5FCV}}$}
\end{axis}
\end{tikzpicture}
\caption{}
\label{prediction_error_power}
\end{subfigure}
\caption{Runtime (in seconds) and $L_2(\mathbb{P}_n)$ prediction error versus sub-sample size for different model selection methods: MDP (\ref{k_tau}), AIC (\ref{k_aic}), GCV (\ref{k_1_out}), and $5$--fold cross-validation (\ref{k_vfcv}), tested on the "California Houses Prices" and "Power Plants" data set. In all cases, each point corresponds to the average of $25$ trials. (a), (c) Runtime verus the sub-sample size $n \in \{420, 525, 700, 1050, 2100 \}$. (b), (d) Least-squares prediction error $\lVert f^k - y_{\textnormal{test}} \rVert$ versus the sub-sample size $n \in \{420, 525, 700, 1050, 2100 \}$.}
\label{large_datasets_results}
\end{figure}




\section{Conclusion} \label{sec:conc}
In the present paper, we tackled the problem of choosing the tuning parameter $k$ in the $k$-NN regression estimator. A strategy based on early stopping and the minimum discrepancy principle was proposed. In Section \ref{sec:optimality_result}, it was shown that the minimum discrepancy stopping rule $k^{\tau}$ (\ref{k_tau}) provides a minimax-optimal estimator, in particular, over the class of Lipschitz functions on a bounded domain. Besides that, this theoretical result was confirmed empirically on artificial and real data sets: the stopping rule has better statistical performance than other model selection rules, such as Hold-out or $5$--fold cross-validation, while reducing the computational time of the model selection procedure. The main inconvenience of the proposed strategy is that one has to estimate the variance $\sigma^2$ of the regression model (as it is the case for the AIC or Mallow's $C_p$ criteria), thus a plug-in estimator is needed. We have constructed such an estimator for simulated experiments with artificial and real data in Section \ref{real_data}.

As for perspectives of this work, we are interested in the theoretical performance of the Nadaraya-Watson regressor (\cite{nadaraya1964estimating, wasserman2006all}). Apparently, it should be close to the one of the k-NN regression estimator since these two non-parametric estimators are closely related (see the monographs \cite{gyorfi2006distribution, tsybakov2008introduction}). The main difficulty should come from the fact that, if $h$ is the bandwidth parameter and $A_h$ is the smoothing matrix of the Nadaraya-Watson estimator, then $\textnormal{tr}(A_h^{\top}A_h) \neq \textnormal{tr}(A_h)$. This fact implies that the expectation of the empirical risk minus the noise variance will not be equal to the difference between the bias and variance terms (see Eq. (\ref{exp_empirical_risk})). Therefore, there should be another concentration result that deals with this problem. Besides that, we should emphasize that the early stopping rules in this work were estimating the famous
bias-variance trade-off (\cite[Chapter 7]{hastie2009elements}). However, recently, (\cite{belkin2019reconciling, belkin2020two}) the bias-variance balancing paradigm was rethought by discovering some settings (exact fit to the data) for which a phenomenon of the "double descent" of the risk curve appeared. It would be interesting to understand if early stopping can work for these settings. The reader can look at a recent paper (\cite{dwivedi2020revisiting}) and references therein for another reexamination of the paradigm.









\printbibliography

\newpage


2021-05-06T02:15:58.000
2008
2008.08718
en
https://arxiv.org/abs/2008.08718
We present a novel data-driven strategy to choose the hyperparameter $k$ in the $k$-NN regression estimator. We treat the problem of choosing the hyperparameter as an iterative procedure (over $k$) and propose using an easily implemented in practice strategy based on the idea of early stopping and the minimum discrepancy principle. This model selection strategy is proven to be minimax-optimal, under the fixed-design assumption on covariates, over some smoothness function classes, for instance, the Lipschitz functions class on a bounded domain. The novel method often improves statistical performance on artificial and real-world data sets in comparison to other model selection strategies, such as the Hold-out method and 5-fold cross-validation. The novelty of the strategy comes from reducing the computational time of the model selection procedure while preserving the statistical (minimax) optimality of the resulting estimator. More precisely, given a sample of size $n$, assuming that the nearest neighbors are already precomputed, if one should choose $k$ among $\left\{ 1, \ldots, n \right\}$, the strategy reduces the computational time of the generalized cross-validation or Akaike's AIC criteria from $\mathcal{O}\left( n^3 \right)$ to $\mathcal{O}\left( n^2 (n - k) \right)$, where $k$ is the proposed (minimum discrepancy principle) value of the nearest neighbors. Code for the simulations is provided atthis https URL.
Machine Learning (stat.ML); Machine Learning (cs.LG); Statistics Theory (math.ST)
Minimum discrepancy principle strategy for choosing $k$ in $k$-NN regression
Qwen/Qwen-72B
1. YES
2. YES