Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
105,372 characters · 19 sections · 96 citation commands
Tuning Parameter Selection in Econometrics
\address[D. Chetverikov]{Department of Economics, UCLA, Bunche Hall, 8283, 315 Portola Plaza, Los Angeles, CA 90095, USA.} \email{[email removed]}
Tuning parameter selection problem appears in many econometric settings. For example, one can think of selecting the number of terms in series estimation, selecting the bandwidth value in local estimation, selecting the penalty parameter in penalized estimation, selecting the number of factors in a factor model, selecting the number of components in a mixture model, and selecting the number of neurons and layers in a neural network. For all these settings, there exist multiple methods to select the corresponding tuning parameters. In this paper, I review some of them. For concreteness, I focus on the series estimation of nonparametric mean regression models and on $\ell_1$-penalized estimation of high-dimensional linear regression models.
In the context of the series estimation, I consider the methods of Mallows, Stein, Lepski, cross-validation, penalization, and aggregation. The methods of Mallows, Stein, and cross-validation are based on the idea of unbiased risk estimation, which consists of selecting the tuning parameter, i.e. the number of series terms in our case, by minimizing an unbiased estimator of the risk of the series estimator. As I discuss, the advantage of the Mallows method is that it works not only for series estimators but in fact for any linear estimators. The Stein method requires the assumption of Gaussian noise but works also for non-linear estimators. Cross-validation is an almost universal method that can be applied in many settings well beyond the nonparametric mean regression model but may not be fully efficient ($V$-fold cross-validation) or may be hard to compute (leave-one-out cross-validation). The Lepski method is based on the idea of starting with a small number of series terms and increasing it up to a point where further increasing the number of terms does not yield a significant reduction of the bias of the estimator, where the significance is measured relative to the variance of the estimator. The penalization method is based on the idea of preventing estimation overfitting by penalizing estimators with too many series terms. Finally, instead of selecting the number of series terms, the aggregation method proposes using a weighted average of estimators with different numbers of series terms, with the idea that averaging the estimators may decrease their variance by averaging out the noise in individual estimators. The methods of Mallows, Stein, cross-validation, penalization, and aggregation are global in the sense that they are built to guarantee good average-over-the-domain performance of the series estimator, i.e. they come with estimation guarantees in the prediction and $\ell_2$ norms; see the next section for the definitions. In contrast, the Lepski method can be tuned to yield good performance at a particular point of interest in the domain. The latter is useful, for example, in regression discontinuity designs.
In the context of $\ell_1$-penalized estimation, I mainly focus on the Lasso estimator of the high-dimensional linear mean regression model. I consider methods to choose the penalty parameter for the Lasso estimator based on the theory of self-normalized moderate deviations, bootstrap, Stein's unbiased risk estimation, and cross-validation. The first two methods come with solid theoretical guarantees and the last two methods are intuitively attractive but require more theoretical work to fully understand their properties. I discuss existing results underlying each of these methods. In addition, I cover not only i.i.d. data setting but also clustered and panel data settings. Moreover, I explore extensions to quantile and generalized linear models.
It is surely impossible to describe all existing work on the problem of tuning parameter selection, and my review is in fact highly selective. In particular, I make no effort to trace the original contributions and, in many cases, give only few references which an interested reader can use as a starting point. Moreover, I omit some existing methods. For example, I do not discuss likelihood-based methods, such as AIC, BIC, and other Bayes-based methods. Also, I do not discuss plugin-type methods, where the tuning parameters are selected by a direct estimation of optimal population analogues, as those methods often lack coherent optimality properties. In addition, I do not cover methods suitable for time series data. Moreover, I focus exclusively on the tuning parameter selection problem from the estimation point of view and do not touch the problem from the inference point of view. I refer an interested reader to other available reviews, e.g. AC10, DTY18, WW20, and BBL23, which are all highly complementary to our discussion here.
In this section, I discuss methods to choose tuning parameters in nonparametric settings. For concreteness, I introduce these methods in the context of the series estimator of the nonparametric mean regression models, but note that most methods to be discussed here can be adapted to other estimation problems as well. I consider the methods of Mallows, Stein, Lepski, cross-validation, penalization, and aggregation. At the end of this section, I also carry out a small Monte Carlo simulation study comparing finite-sample performance of these methods.
I consider the nonparametric mean regression model
where $Y\in\mathbb{R}$ is a dependent variable, $X\in\mathcal X$ is a vector of covariates, $e \in \mathbb{R}$ is a zero-mean noise, and $f\colon \mathcal X\to \mathbb{R}$ is a function of interest. I assume that our task is to estimate $f$ using a random sample $(X_1,Y_1),\dots,(X_n,Y_n)$ from the distribution of the pair $(X,Y)$. To describe the series estimator of the function $f$, for all $k\geq 1$, let $p^k = (p_1,\dots,p_k)^{\top}$ be a vector of $k$ basis functions mapping $\mathcal X$ into $\mathbb{R}$ and chosen so that a suitable linear combination of these functions can approximate $f$. Here, all $k$ functions in the vector $p^k$ are allowed to depend on $k$ but I omit this dependence for simplicity of notation. Assuming that $k$ is large enough, the approximation is expected to be sufficiently good in the sense that there exists a vector $\beta_k\in\mathbb{R}^k$ such that
where the residual $r_k$ is small. The series estimator is then based on the idea of OLS estimation of the vector of coefficients $\beta_k$ pretending that the residual $r_k$ is non-existent:
When $X$ is a scalar random variable, the basis functions could be monomials: $p_j(x)=x^{j-1}$ for $j=1,\dots,k$, or quadractic regression splines: $p_1(x)=1$, $p_2(x)=x$, $p_3(x)=x^2$, $p_j(x)=((x-(j-3)/(k-2))\vee 0)^2$ for $j=4,\dots,k$, if we further assume that the support of $X$ is the $[0,1]$ interval. Note that in the case of monomials, all functions $p_j$ are independent of $k$ but in the case of regression splines, only $p_1$, $p_2$, and $p_3$ are independent of $k$. When $X$ is a random vector, the basis functions could be given by tensor products of univariate monomials or regression splines, e.g. see BCCK15 for details and further examples of popular basis functions.
The series estimator $\widehat f_k$ in (ref) contains a tuning parameter $k$, also known as a smoothing or regularization parameter, and one of the most important questions in the theory of nonparametric estimation is how this parameter can be chosen in practice. Intuitively, when $k$ is too small, the residual $r_k$ in (ref) may not be negligible and the estimator $\widehat f_k$ may not be flexible enough to capture the shape of the function $f$. When $k$ is too large, we have to estimate too many parameters, which leads to the estimator $\widehat f_k$ with large variance. We thus want to choose $k$ in a way to strike a balance between flexibility of the estimator $\widehat f_k$ and its variance. Formally speaking, we would like to construct an estimator $\widehat k$ such that
where $d(\widehat f_{\widehat k},f)$ measures the distance between the estimator $\widehat f_{\widehat k}$ and the true function $f$, and $\mathcal K_n$ is a set of candidate values of the tuning parameter $k$ that is chosen large enough to include a value striking a good balance between flexibility and variance. Here, the set $\mathcal K_n$ is typically assumed to depend on $n$ and to grow together with $n$. For example, when the basis functions are given by monomials, we could consider $\mathcal K_n = \{1,\dots,\bar k_n\}$, where $\bar k_n = [n^{1/3}]$, the largest integer smaller than $n^{1/3}$. When the basis functions are given by quadratic splines, we could consider $\mathcal K_n = \{3,\dots,\bar k_n\}$, where $\bar k_n = [n^{1/2}]$, the largest integer smaller than $n^{1/2}$.\footnote{As discussed in BCCK15, $\bar k_n$ should be chosen so that $\bar k_n^2\log n = o(n)$ in the case of monomials and $\bar k_n\log n = o(n)$ in the case of quadratic splines in order to guarantee existence of the series estimator. In addition, to avoid numerical problems with inverting the matrix $\sum_{i=1}^n p^k(X_i)p^k(X_i)^{\top}$, it is recommended to consider a suitable linear transformation of the functions $p^k$. For example, one can use Legendre polynomials instead of monomials and B-splines instead of regression splines.}
The following four distance measures are typically considered in the literature:
where the supremum in the uniform metric is taken either over all $x$ in $\mathcal X$ or over all $x$ in some user-specified subset $\bar{\mathcal X}$ of $\mathcal X$, and $x_0$ in the pointwise metric is a user-specified value in $\mathcal X$. A result in the $\ell_2$ metric is typically useful in the semi-parametric/double-machine-learning estimation, where the estimator $\widehat f_{\widehat k}$ is used as a first step in a two-step procedure to learn some causal parameters, e.g. see CCDDHNR18. A result in the prediction metric is useful in signal processing, where the tast is to get rid of the noise $e_i$ in observable signals $Y_i$. A result in the uniform metric is useful when we want to understand the shape of the function $f$. A result in the pointwise metric is useful when we are interested in a particular point $x_0$ in the domain of the function $f$. Note that the results of the form (ref) in the $\ell_2$ and prediction metrics are closely related, as it is typically possible to translate one into the other by applying the law of large numbers.
When the estimator $\widehat f_{\widehat k}$ satisfies (ref), it is said to be {\em asymptotically optimal}. However, the result of the form (ref) may not always be achievable. Intuitively, a problem occurs when there exists $k\in\mathcal K$ such that $d(\widehat f_k,f)$ is too small, which is the case when a linear combination of the functions $p_1,\dots,p_k$ yields the function $f$ exactly. Alternatively, the result (ref) may be achievable but the rate of convergence to the constant one may be slow. For these reasons, we are also interested in constructing an estimator $\widehat k$ such that for some constant $C\geq 1$ and a sequence $\varepsilon_n$ converging to zero, we have either
with probability approaching one or
Such inequalities are called {\em oracle inequalities} because they show the performance of the estimator $\widehat f_{\widehat k}$ based on our choice of $k$ relative to the performance of the estimator $\widehat f_{k}$ based on the oracle choice $k = k^O$, where $$ k^O = \arg\min_{k\in\mathcal K_n} d(\widehat f_k,f), $$ The purpose of this section is to discuss estimators $\widehat k$ that lead to results in the forms (ref), (ref), (ref), or in some related forms.
For our discussion in this section, it will be convenient to use matrix notations. I therefore denote $\mathrm{Y} = (Y_1,\dots,Y_n)^{\top}$, $\mathrm{f} = (f(X_1),\dots,f(X_n))^{\top}$, and $\mathrm{e} = (e_1,\dots,e_n)^{\top}$, so that $\mathrm{Y} = \mathrm{f} + \mathrm{e}$. Also, I denote $\widehat{\mathrm{f}}_k = (\widehat f_k(X_1),\dots,\widehat f_k(X_n))^{\top}$, $\mathrm{P}_k = (p^k(X_1),\dots,p^k(X_n))^{\top}$, and $\mathrm{M}_k = \mathrm{P}_k(\mathrm{P}_k^{\top}\mathrm{P}_k)^{-1}\mathrm{P}_k^{\top}$, so that $\widehat{\mathrm{f}}_k = \mathrm{M}_k \mathrm{Y}$. In addition, I let $\mathrm{I}_k$ be a $k\times k$ identity matrix. Moreover, I denote $\sigma^2_i = \mathbb{E}[e_i^2\mid X_i]$ for all $i=1,\dots,n$ and let $\Omega$ be the $n\times n$ diagonal matrix with $(\sigma_1^2,\dots,\sigma_n^2)$ on the diagonal. Finally, for any function $g\colon \mathcal X\to\mathbb{R}$, I denote $\|g\|_{2,n}=(n^{-1}\sum_{i=1}^n g(X_i)^2)^{1/2}$.
M73 proposed a method that in our context of choosing $k$ for the series estimator of model (ref) takes the following form:
To explain the intuition behind this method, observe that since $\widehat{\mathrm{f}}_k = \mathrm{M}_k\mathrm{Y} = \mathrm{M}_k \mathrm{f} + \mathrm{M}_k \mathrm{e}$, we have $$ \mathbb{E}[\|\widehat{\mathrm{f}}_k - \mathrm{f}\|^2] = \mathbb{E}[\| (\mathrm{M}_k - \mathrm{I}_k)\mathrm{f} + \mathrm{M}_k \mathrm{e} \|^2] = \mathbb{E}[\| (\mathrm{M}_k - \mathrm{I}_k)\mathrm{f} \|^2] + \mathbb{E}[\|\mathrm{M}_k \mathrm{e} \|^2]. $$ Also, given that $\mathrm{Y} = \mathrm{f} + \mathrm{e}$ and $\widehat{\mathrm{f}} = \mathrm{M}_k \mathrm{f} + \mathrm{M}_k \mathrm{e}$, we have
Here,
Moreover, $\mathbb{E}[\|\mathrm{e}\|^2] = \mathbb{E}[\mathrm{tr}(\Omega)]$ and $$ \mathbb{E}[\mathrm{e}^{\top}\mathrm{M}_k \mathrm{e}] = \mathbb{E}[\mathrm{tr}(\mathrm{e}^{\top}\mathrm{M}_k \mathrm{e})] = \mathbb{E}[\mathrm{tr}(\mathrm{M}_k \mathrm{e}\mathrm{e}^{\top})] = \mathbb{E}[\mathrm{tr}(\mathrm{M}_k\Omega_k)] $$ by interchanging the trace and the expectation operators and using the law of iterated expectations. Thus,
The left-hand side here is the expectation of the criterion function in (ref) and the right-hand side is often referred to as the risk of the estimator $\widehat f_k$, which we want to minimize. The Mallows method (ref) thus consists of minimizing an unbiased risk estimator. In addition, unless the set $\mathcal K_n$ is too large, standard concentration inequalities, e.g. BLM13, can be used to show that the random variables $$ \frac{1}{n}\sum_{i=1}^n (Y_i - \widehat f_k(X_i))^2 + \frac{2}{n}\mathrm{tr}(\mathrm{M}_k\Omega) - \frac{1}{n}\mathrm{tr}(\Omega)\quad\text{and}\quad \| \widehat f_k - f \|^2_{2,n} $$ converge in probability to their expectations uniformly over $k$. Convergence of the former random variable, together with (ref) and (ref), implies that $\widehat k^M$ nearly minimizes the risk $\mathbb{E}[\| \widehat f_k - f \|^2_{2,n}]$. Convergence of the latter random variable implies that $\widehat k^M$ also nearly minimizes $\| \widehat f_k - f \|^2_{2,n}$. Formally, under appropriate regularity conditions, we can then obtain the following result:
see L87 and A91 for details in the homoskedastic and heteroskedastic cases, respectively. Thus, the Mallows method yields an estimator that is asymptotically optimal in the prediction metric. Moreover, as discussed above, (ref) can be combined with the uniform law of large numbers to obtain asymptotic optimality in the $\ell_2$ metric. Related results were also obtained in Shi81, PT90, K94.
Note, however, that the appropriate regularity conditions here include the requirement that $n\|\widehat f_k - f\|_{2,n}^2\to \infty$ for all $k\in\mathcal K_n$, and this requirement is not satisfied when $f(\cdot) = p^k(\cdot)^{\top}\beta_k$ for some $k\in\mathcal K_n$ and $\beta_k \in \mathbb{R}^k$, as follows from the standard OLS estimation theory. To obtain a more general result for the series estimator $\widehat f_k$ with the number of series terms $k$ selected by the Mallows method, we thus need to proceed through oracle inequalities. We will return to this topic in Section (ref), where we discuss the penalization method.
The Mallows method, as described above, is generally infeasible as it requires the knowledge of the heteroskedasticity matrix $\Omega$. However, it is fairly straightforward to obtain a feasible version of this method. Indeed, observe that the term $n^{-1}\mathrm{tr}(\Omega)$ in the criterion function in (ref) is independent of $k$ and thus can be omitted. In addition, given that $\mathrm{M}_k=\mathrm{P}_k(\mathrm{P}_k^{\top}\mathrm{P}_k)^{-1}\mathrm{P}_k$, we have
which can be consistently estimated by
where $\widehat e_i = Y_i - \widehat f_{\bar k}(X_i)$ for all $i=1,\dots,n$ and some $\bar k\in\mathcal K_n$ that is known to yield a consistent estimator $\widehat f_{\bar k}$ of the function $f$, e.g. $\bar k = \bar k_n = [n^{1/3}]$. A feasible version of the Mallows method is thus given by
Moreover, note that if we have a model with homoskedastic noise, i.e. $\sigma_i^2 = \sigma^2$ for all $i=1,\dots,n$ and some $\sigma^2>0$, the Mallows estimator $\widehat k^M$ in (ref) reduces to
as the expression in (ref) in this case reduces to $\mathrm{tr}((\mathrm{P}_k^{\top}\mathrm{P}_k)^{-1}(\mathrm{P}_k^{\top}\mathrm{P}_k)) = \mathrm{tr}(\mathrm{I}_k) = k$. We will use this formula later to relate the Mallows and penalization methods.
Finally, note that the intuition behind the estimator $\widehat k^M$ in (ref) does not really depend on the properties of the series estimator $\widehat \mathrm{f}_k = \mathrm{M}_k\mathrm{Y}$ except for its linearity in $\mathrm{Y}$, and so the Mallows method can in fact be used to choose the regularization parameter $k$ for other linear in $\mathrm{Y}$ estimators of the function $f$, e.g. ridge estimators. In addition, the Mallows method is often used for aggregating linear estimators, i.e. for constructing estimators of the form $\sum_{k\in\mathcal K_n} w_k\widehat f_k$, where the weights $w_k$ are selected by solving an optimization problem analogous to that in (ref), e.g. see H07 in the context of the series estimator and ZWZZ19 and ZWZL23 in some related contexts.
S81 proposed the following estimator of $k$:
where $D^k$ is the $n\times n$ matrix defined by $D^k_{ij} = \partial\widehat f_k(X_i)/\partial Y_j$ for all $i,j=1,\dots,n$. This method is only justified when the conditional distribution of $e$ given $X$ is Gaussian, which we are going to assume throughout this section. To see the intuition behind this method, observe that if $Z\sim N(\mu,\sigma^2)$ and $g$ is a (weakly) differentiable function, then $\sigma^2\mathbb{E}[g'(Z)] = \mathbb{E}[(Z-\mu)g(Z)]$ via integration by parts. Applying this identity, which is often referred to as Stein's lemma, conditional on $\mathrm X = (X_1,\dots,X_n)$ to the function $Y_i\mapsto \widehat f_k(X_i) - Y_i$ yields $$ \sigma_i^2\mathbb{E}\left[\frac{\partial}{\partial Y_i}(\widehat f_k(X_i) - Y_i)\mid \mathrm X\right] = \mathbb{E}[(Y_i - f(X_i))(\widehat f_k(X_i) - Y_i)\mid \mathrm X],\quad \text{for all }i=1,\dots,n. $$ Thus,
Hence, like the the Mallows method, the Stein method minimizes an unbiased risk estimator. Using the same arguments as those used to analyze the Mallows method, it is thus possible to show that the Stein method yields an asymptotically optimal estimator $\widehat f_{\widehat k^S}$ in the prediction and $\ell_2$ metrics. In the literature, this method is often referred to as the Stein's unbiased risk estimation, or SURE, for short.
The advantage of the Stein method over the Mallows method is that it applies for non-linear in $\mathrm{Y}$ estimators as well. In fact, it only requires that the estimators $\widehat f_k$ are weakly differentiable in $\mathrm{Y}$ and can be used to choose regularization parameters of estimators as complicated as the Lasso estimator, which we are going to discuss in Section (ref) below. In addition, for example, MW00 applied the Stein method to estimation problems with shape restrictions and CST13 applied it to low-rank matrix recovery problems. On the other hand, the disadvantage of the Stein method is that it requires that the conditional distribution of $e$ given $X$ is Gaussian, and it seems generally difficult to relax this condition, unless the estimators $\widehat f_k$ are sufficiently smooth in $\mathrm{Y}$.
It is also interesting to note that in our special case of the series estimator, the Stein method yields $\widehat k^S$ that is identical to $\widehat k^M$ given by the Mallows method. Indeed, given that $$ \widehat f_k(X_i) = p^k(X_i)^{\top}\left(\sum_{j=1}^n p^k(X_j)p^k(X_j)^{\top}\right)^{-1}\sum_{j=1}^n p^k(X_j)Y_j,\quad\text{for all }i=1,\dots,n, $$ we have $$ D_{ii}^k = p^k(X_i)^{\top}\left(\sum_{j=1}^n p^k(X_j)p^k(X_j)^{\top}\right)^{-1}p^k(X_i),\quad \text{for all }i=1,\dots,n, $$ and so $\mathrm{tr}(D^k\Omega) = \mathrm{tr}(\mathrm{M}_k\Omega)$, meaning that the criterion functions in (ref) and (ref) coincide. In particular, this implies that a feasible version of the Stein method can be obtained as discussed in the case of the Mallows method.
L90, L91, L92 proposed a method for choosing regularization parameters for nonparametric estimators that has an advantage of yielding results not only in the prediction and $\ell_2$ metrics but also in the uniform and pointwise metrics. The method attracted substantial attention in the literature and many variants of the method have been developed. Here, we describe a slight modification of the variant proposed in SW19. For concreteness, I focus on results in the pointwise metric.
Suppose that we are interested in estimating $f(x_0)$ for some particular $x_0\in\mathcal X$ and our task is to choose one estimator from a collection of the series estimators $\{\widehat f_k(x_0)\}_{k\in\mathcal K_n}$. Intuitively, as we increase the number of series terms $k$, the (absolute value of the) bias $|\mathbb{E}[\widehat f_k(x_0) - f(x_0)]|$ of the estimator $\widehat f_k(x_0)$ goes down and its variance $\mathrm{var}(\widehat f_k(x_0))$ goes up. The main idea of the Lepski method is to start with a small value of $k$ and increase it up to a point where further increasing it does not significantly decrease the bias of the estimator $\widehat f_k(x_0)$, where the significance is understood relative to the variance of the estimator. To formalize this intuition, for all $k,k'\in\mathcal K_n$, denote $$ b_{k,k'} = \mathbb{E}[\widehat f_{k'}(x_0) - \widehat f_{k}(x_0)]\quad\text{and}\quad p_{k,k'} = \mathrm{var}(\widehat f_{k'}(x_0) - \widehat f_{k}(x_0)). $$ Also, let $\beta>0$ be some number and let $\mathcal K^+_n(k) = \{k'\in\mathcal K_n\colon k'>k\}$. Then for each $k\in\mathcal K_n$, we can consider the problem of testing the null hypothesis $$ \mathrm H_{k}\colon |b_{k,k'}|^2 \leq \beta^2p_{k,k'},\quad\text{for all }k'\in\mathcal K_n^+(k). $$ If the hypothesis $\mathrm H_{k}$ is not rejected, there is no reason to set the number of series terms above $k$. The method thus consists of choosing the number of series terms as the smallest number $k\in\mathcal K_n$ such that $\mathrm H_k$ is not rejected.\footnote{The original Lepski method considered choosing the number of series terms as the smallest number $k\in\mathcal K_n$ such that $\mathrm H_{k'}$ for all $k'\geq k$ are not rejected. The more relaxed choice described here was proposed in SV09.}
It thus remains to discuss appropriate tests for hypotheses $\mathrm H_{k}$. One option is as follows. Denote $T_{k,k'}=|\widehat f_{k'}(x_0) - \widehat f_{k}(x_0)|$ for all $k'\in\mathcal K_n^+(k)$. Then under the null hypothesis $\mathrm H_{k}$, we have for all $k'\in\mathcal K_n^+(k)$ that $|b_{k,k'}|\leq \beta \sqrt{p_{k,k'}}$ and so
Thus, letting $c_{k}(\alpha)$ be the $(1-\alpha)$ quantile of
a test of $\mathrm H_{k}$ with size $\alpha$ consists of rejecting $\mathrm H_{k}$ if $$ \max_{k'\in\mathcal K_n^+(k)} \frac{T_{k,k'}}{\sqrt{p_{k,k'}}} > \beta + c_{k}(\alpha). $$ Let $\widehat k^L$ be the value of $k\in\mathcal K_n$ selected by this method.
To understand properties of the estimator $\widehat f_{\widehat k^L}(x_0)$, let $k^*$ be the smallest value of $k\in\mathcal K_n$ such that the hypothesis $\mathrm H_k$ is true. We can think of $\widehat f_{k^*}(x_0)$ as the oracle estimator, and our main task is to derive a bound on how much the estimator $\widehat f_{\widehat k^L}(x_0)$ deviates from this oracle estimator.\footnote{The name “oracle estimator” comes from the observation that we can think of $\widehat f_{k^*}(x_0)$ as the estimator provided by the oracle, who does not need to test hypotheses $\mathrm H_k$ as it knows whether these hypotheses are true or not.} To do so, observe that by construction, the hypothesis $\mathrm H_{k^*}$ is accepted with probability at least $1-\alpha$, and so the probability that $\widehat k^L\leq k^*$ is also at least $1 - \alpha$. Moreover, letting $\mathcal K_n^-(k^*) = \{k\in\mathcal K_n\colon k<k^*\}$, we have for any $k\in \mathcal K_n^-(k^*)$ that
Hence, letting $c_{k^*}^-(\alpha)$ be the $(1-\alpha)$ quantile of $$ \max_{k\in\mathcal K_n^-(k^*)}\frac{|\widehat f_{k}(x_0) - \mathbb{E}[\widehat f_{k}(x_0)] - \widehat f_{k^*}(x_0) + \mathbb{E}[\widehat f_{k^*}(x_0)]|}{\sqrt{p_{k,k^*}}} $$ and $$ \mathcal K_n^0 = \left\{k\in \mathcal K_n^-(k^*)\colon |b_{k,k^*}|^2 > p_{k,k^*}(c_k(\alpha) + c_{k^*}^-(\alpha))^2\right\}, $$ we have $$ \frac{T_{k,k^*}}{\sqrt{p_{k,k^*}}} > c_{k}(\alpha),\quad\text{for all }k\in\mathcal K_n^0, $$ with probability at least $1-\alpha$. Thus, by the union bound, $\widehat k^L\in \mathcal K_{IN}$ with probability at least $1 - 2\alpha$, where $\mathcal K_{IN} = \{k\in\mathcal K_n\colon k\leq k^*\}\setminus \mathcal K_n^0$. In turn,
with the same probability.
The quantity on the right-hand side of (ref) shows how much the series estimator $\widehat f_{\widehat k^L}(x_0)$ with the number of series terms determined by the Lepski method deviates from the oracle estimator $\widehat f_{k^*}(x_0)$ and can be understood as the adaptation price, as it shows how much we need to pay in terms of the precision of the estimator for not knowing the oracle choice $k^*$. To understand this quantity, observe first that $c_k(\alpha)$ is typically of order $\sqrt{\log|\mathcal K_n|} \leq \sqrt{\log n}$, where $|\mathcal K_n|$ denotes the number of elements in the set $\mathcal K_n$. Also, we typically have that $p_{k,k^*}$ is of order $(k^*-k)/n$. Thus, the size of the adaptation price depends on the set $\mathcal K_{IN}$, which SW19 refer to as the {\em insensitivity region}. If the function $k\mapsto b_{k,k^*}$ is relatively flat around $k=k^*$, this region may be rather large, in which case we expect the adaptation price to be of order $\sqrt{k^*\log n/n}$. Since the oracle estimation error $\widehat f_{k^*}(x_0) - f(x_0)$ is typically of order $\sqrt{k^*/n}$, it follows that the adaptation price in this case is actually larger than the oracle estimation error itself. On the other hand, if the function $k\mapsto b_{k,k^*}$ is relatively steep around $k=k^*$, the insensitivity region may be rather small, in which case the adaptation price may be as small as $\sqrt{\log n/n}$ in order, and may be much smaller than the oracle estimation error.
Next, observe that implementing the Lepski method requires the choice of two parameters, $\alpha$ and $\beta$, and one could potentially argue that the Lepski method is strange in the sense that it replaces the problem of choosing one tuning parameter, $k$, by the problem of choosing two tuning parameters, $\alpha$ and $\beta$. However, these parameters are very different from $k$ as they do not generate the bias-variance trade-off. Instead, these parameters specify the properties of the method. Indeed, the purpose of $\beta$ is to fix the oracle estimator we are potentially interested in, and the choice $\beta = 1$ seems reasonable. Similarly, the purpose of $\alpha$ is to specify the probability for theoretical guarantees, and is closely related to the purpose of the size control parameter in testing problems. The choices $\alpha = 0.05$ and $\alpha = 0.01$ seem reasonable.
The Lepski method as described here yields results in the pointwise metric. However, it is fairly straightforward to extend the method to yield results in other metrics as well, which is its main advantage. Indeed, if we want to have results in the uniform metric, we simply need to replace the statistic $T_{k,k'} = |\widehat f_{k'}(x_0) - \widehat f_k(x_0)|$ by $T_{k,k'} = \sup_{x\in\bar{\mathcal X}}|\widehat f_{k'}(x) - \widehat f_{k}(x)|$, and if we want the results in the $\ell_2$ metric, we should use $T_{k,k'} = (\mathbb{E}_X[(\widehat f_{k'}(X) - \widehat f_k(X))^2])^{1/2}$ or its estimated version, where the expectation is taken with respect to the distribution of $X$ that is independent of the data. Of course, in both case, the critical values for the tests should be adjusted accordingly. Note that other metrics are possible here as well. On the other hand, a disadvantage of the method is that it yields results relative to the oracle number of series terms $k^*$ rather than relative to the more relevant oracle number $k^{O} = \arg\min_{k\in\mathcal K}|\widehat f_k(x_0) - f(x_0)|$ appearing above.
Further, the variant of the method described so far is infeasible as it requires the knowledge of $p_{k,k'}$'s and $c_k(\alpha)$'s. However, these quantities are easy to estimate. Indeed, denoting $\mathrm Q_k = n^{-1}\sum_{i=1}^n p^k(X_i)p^k(X_i)^{\top}$ for all $k\in\mathcal K_n$ and letting $\bar k$ be some value of $k\in\mathcal K_n$ that is known to yield a consistent estimator $\widehat f_{\bar k}$ of the function $f$, e.g. $\bar k = \bar k_n = [n^{1/3}]$, and denoting $\widehat e_i = Y_i - \widehat f_{\bar k}(X_i)$ for all $i=1,\dots,n$, it follows that we can estimate $p_{k,k'}$ by $$ \widehat p_{k,k'} = \frac{1}{n^2}\sum_{i=1}^n \widehat e_i^2 \Big(p^{k'}(x_0)^{\top}\mathrm Q_{k'}^{-1}p^{k'}(X_i) - p^{k}(x_0)^{\top}\mathrm Q_{k}^{-1}p^{k}(X_i)\Big)^2 $$ for all $k,k'\in\mathcal K_n$. Also, $c_k(\alpha)$ can be estimated, for example, by the multiplier bootstrap: let $\{w_i\}_{i=1}^n$ be i.i.d. standard Gaussian random variables and let
be the bootstrap version of the random variable in (ref). Then $c_k(\alpha)$ can be estimated by $\widehat c_k(\alpha)$, the $(1-\alpha)$ quantile of the conditional distribution of the random variable in (ref) given the data. Consistency of (the properly normalized versions of) the estimators $\widehat p_{k,k'}$ and $\widehat c_{\alpha}(k)$ can be established using the results in BCCK15 and CCK13, CCK17, CCKK22.
Finally, I note that the Lepski method is nearly universal in the sense that it can be adjusted to work in other estimation frameworks as long as the estimators in these frameworks feature a bias-variance trade-off and can be ordered in terms of decreasing bias. For example, CC18 and CCKan21 used this method in the context of nonparametric instrumental variable estimation, SWH13 used it in the context of nonparametric quantile regression, and CCK14, among many others, used it in the context of density estimation.
Cross-validation, which is traditionally attributed to the work of A74, S74, and G75, is one of the most universal and popular methods to choose tuning parameters in general and in the context of mean regression models via the series estimator in particular. Below, we will consider its three main versions: validation, $V$-fold cross-validation, and leave-one-out cross-validation.
To describe the validation method, let $(I_1,I_2)$ be a random split of the whole sample $I_0 = \{1,\dots,n\}$ into two subsamples, so that $I_1\cap I_2 = \emptyset$ and $I_1\cup I_2 = I_0$. Here, $I_1$ and $I_2$ are referred to as the training subsample and the validation subsample, respectively. Typically, $I_1$ and $I_2$ are assumed to contain roughly two thirds and one third of the whole sample, respectively, but this will not be important for our discussion. For each $k$ in the set of candidate values $\mathcal K_n$, we build a series estimator of the function $f$ using the training subsample only: $\widehat f_{k,1}(\cdot) = p^k(\cdot)^{\top}\widehat\beta_k$, where $$ \widehat\beta_k = \left(\sum_{i\in I_1}p^k(X_i)p^k(X_i)^{\top}\right)\sum_{i\in I_1}p^k(X_i)Y_i. $$ We then choose $k$ by minimizing the out-of-sample prediction error of these estimators:
The ultimate estimator of the function $f$ is then typically $\widehat f_V(\cdot) = \widehat f_{\widehat k^V,1}(\cdot)$.
Theoretical properties of the validation procedure are studied, among others, in W03 and M07. Intuitively, conditional on the data in the training subsample $\{(X_i,Y_i\}_{i\in I_1}$, the sum $\sum_{i\in I_2}(Y_i - \widehat f_{k,1}(X_i))^2$ consists of i.i.d. terms with mean $$ \mathbb{E}_{X,Y}[(Y - \widehat f_{k,1}(X))^2] = \mathbb{E}_{X,Y}[e^2] + \mathbb{E}_{X,Y}[(\widehat f_{k,1}(X) - f(X))^2], $$ where $\mathbb{E}_{X,Y}[\cdot]$ denotes the expectation with respect to the distribution of the pair $(X,Y)$ holding the data $\{(X_i,Y_i\}_{i\in I_1}$, and thus $\widehat f_{k,1}$, fixed. Thus, $\widehat k^V$ in (ref) approximately minimizes $\mathbb{E}_{X,Y}[(\widehat f_{k,1}(X) - f(X))^2]$ and can be shown to lead to a result in the form of (ref) in the $\ell_2$ metric with $\widehat f_k$ replaced by $\widehat f_{k,1}$.
Next, although the validation procedure described above is very intuitive, an issue with this procedure is that it yields the result for the subsample estimators $\widehat f_{k,1}$ instead of the full-sample estimators $\widehat f_k$, which are inefficient because they do not use the validation subsample in estimation. To fix this issue, we can consider $V$-fold cross-validation. To describe it, let $(I_1,\dots,I_V)$ be a random partition of the whole sample $I_0 = \{1,\dots,n\}$ into $V$ subsamples of roughly the same size. For each $v=1,\dots,V$, we build a series estimator of the function $f$ using all the data excluding observations in subsample $v$: $\widehat f_{k,v}(\cdot) = p^k(\cdot)^{\top}\widehat\beta_{k,v}$, where $$ \widehat\beta_{k,v} = \left(\sum_{i\in I_0\setminus I_v} p^k(X_i)p^k(X_i)^{\top}\right)^{-1}\sum_{i\in I_0\setminus I_v}p^k(X_i)Y_i. $$ We then choose the tuning parameter $k$ by minimizing the sum of out-of-sample prediction errors of these estimators: $$ \widehat k^{VCV} = \arg\min_{k\in\mathcal K}\sum_{v=1}^V\sum_{i\in I_v} (Y_i - \widehat f_{k,v}(X_i))^2. $$ The ultimate estimator of the function $f$ can be either the average estimator $$ \widehat f_{AVCV}(\cdot) = \frac{1}{V}\sum_{v=1}^V \widehat f_{\widehat k^{VCV},v}(\cdot) $$ or the full-sample estimator $\widehat f_{FVCV}(\cdot) = \widehat f_{\widehat k^{VCV}}(\cdot)$, with the estimators $\widehat f_k(\cdot)$ being described in (ref). The analysis of the average estimator can be found in LM12, and the full-sample estimator can be analyzed using the arguments in CL21, who studied the quantile version of the estimator $\widehat f_{FVCV}$.
The estimators $\widehat f_{AVCV}$ and $\widehat f_{FVCV}$ are more efficient then $\widehat f^V$ but are still not asymptotically optimal. To obtain an asymptotically optimal estimator, we should assume that the number of folds $V$ gets large together with the sample size $n$. In particular, the arguments in CL21 can be used to show that as $V$ gets large, the full sample estimator $\widehat f_{FVCV}$ becomes asymptotically optimal in the sense that it satisfies (ref) in the $\ell_2$ metric. In the extreme case, when $V = n$, we obtain the leave-one-out cross-validation, which was studied in L87 in the context of the mean regression (ref) and C14 in the context of density estimation, among others. The problem with the leave-one-out cross-validation, however, is that it is computationally difficult, as it requires calculating $n$ estimators for each $k$. On the other hand, in the case of the series estimator, computation turns out to be not so difficult because of the existence of an explicit expression: when $V=n$, the expression for $\widehat k^{VCV}$ above reduces to $$ \widehat k^{LOO} = \arg\min_{k\in\mathcal K}\frac{1}{n}\sum_{i=1}^n \frac{(Y_i - \widehat f_k(X_i))^2}{(1 - m_k)^2}, $$ where $m_k$ is the $k$th diagonal element of the matrix $\mathrm{M}_k$; see equation (1.5) in L87.
More generally, cross-validation procedures can be applied in many settings, e.g. see BBL23, but they always give results in the $\ell_2$ metric and can not be applied if we are interested in the uniform or pointwise metrics. See AC10 for a deep review of the theory underlying cross-validation procedures. See also L20 for more recent results.
The penalization method is based on the idea that the series estimator with too many terms $k$ suffers from overfitting the data, i.e. it yields a small value of the least squares criterion function but has large variance. To avoid choosing too many series terms, the penalization method adds a penalty to the least squares criterion function that is proportional to $k$:
where $\lambda>0$ is the penalty parameter.
The penalization method thus generalizes the Mallows method under homoskedasticity, which corresponds to setting $\lambda = 2\sigma^2$ in (ref); see the expression in (ref). In fact, the Mallows method can be seen as a version of the penalization method where the penalty parameter $\lambda$ is chosen so that the resulting criterion function yields an unbiased risk estimator for the series estimator. The Mallows method yields the series estimator $\widehat f_{\widehat k^M}$ that is asymptotically optimal, and so setting $\lambda = 2\sigma^2$ under homoskedasticity indeed seems reasonable. The benefit of having a flexible penalty parameter, however, becomes apparent when we generalize the problem of selecting the number of terms $k$ for the series estimator and consider the general problem of {\em model selection}. To describe it, consider a set of models $\mathcal M_n$ such that each model $m\in\mathcal M_n$ corresponds to some set of basis functions $\{p_{m,1},\dots,p_{m,k(m)}\}$. Then for each model $m\in\mathcal M_n$, we have a corresponding series estimator $\widehat f_m(\cdot)$ defined by analogy with (ref) where we replace $p^k$ by $p^m = (p_{m,1},\dots,p_{m,k(m)})^{\top}$. We are then interested in constructing an estimator $\widehat m\in\mathcal M_n$ such that $d(\widehat f_m,f)$ is as small as possible in general and satisfies analogs of (ref), (ref), or (ref) in particular. The model selection problem generalizes the problem of selecting the number of terms $k$ for the series estimator in the sense that the latter corresponds to the model selection problem where for each $k\in\mathcal K_n$, we have only one model $m\in\mathcal M_n$ such that $k(m) = k$, and the models are assumed to be ordered in terms of their ability to approximate the true function $f$. The model selection problem also generalizes the {\em variable selection} problem, where we are given a set of basis functions (or variables) $\{p_1,\dots,p_{N_n}\}$ and each model $m$ corresponds to a subset of this set. We say that we have an {\em ordered variable selection} problem if the models are nested, e.g. $\mathcal M_n = \{1,\dots,N_n\}$ and each model $m$ corresponds to the basis functions $p_1,\dots,p_m$, and we say that we have a {\em complete variable selection} problem if all subsets of the set $\{p_1,\dots,p_{N_n}\}$ are included in the set of models. Clearly, the ordered variable selection problem is closely related to the problem of selecting the number of terms for the series estimator. In fact, two problems coincide if the basis functions $p_j$ in the vector $p^k$ are independent of $k$, which is the case when the basis functions are given by monomials, for example.
It turns out that the penalization method can be extended to work with the model selection problem as well:
where $\lambda_m$'s are model specific penalty parameters and $|m|= k(m)$ denotes the number of basis functions in the model $m$. We will see below why it makes sense to consider model-specific penalty parameters.
It is obvious that the Mallows, Stein, and cross-validation methods discussed above can all be extended in a straightforward manner to work with the model selection problem as well but the penalization approach is expected to work better if the set of models $\mathcal M_n$ is large. For example, the Mallows method for model selection under homoskedastisticity would set $\lambda_m = 2\sigma^2$ for all $m\in\mathcal M_n$ in (ref), which makes the criterion function an unbiased risk estimator of the estimator $\widehat f_m$. However, as demonstrated in BM07, this choice is inappropriate if we have too many models in the set $\mathcal M_n$. In particular, it is inappropriate if we are dealing with the complete variable selection problem. The issue here is that even though the criterion function yields an unbiased risk estimator for each model $m$, when we are dealing with many models simultaneously, it is likely to happen that at least for some models, the values of the criterion function are much smaller than the risk of the corresponding series estimators, yielding misleading risk estimates. In fact, when dealing with the complete variable selection problem, setting $\lambda_m = 2\sigma^2$ for all $m\in\mathcal M_n$ may be so bad that the resulting estimator $\widehat f_{\widehat m}$ may not be even consistent. Thus, the Mallows, Stein, and cross-validation methods, which are all based on the idea of unbiased risk estimation, may not be appropriate if the set $\mathcal M_n$ is too large in general and if we are dealing with the complete variable selection problem in particular.
To describe the properties of the penalized model selection estimator $\widehat m$ in (ref), let $\{L_m\}_{m\in\mathcal M_n}$ be a sequence of non-negative numbers and let $\theta\in(0,1)$ and $\kappa>2 - \theta$ be some constants. Assuming Gaussian and homoskedastic noise, $e\mid X\sim N(0,\sigma^2)$, BM07 proved the following result: If we set
then the estimator $\widetilde f = \widehat f_{\widehat m}$ satisfies
where $$ \varepsilon_n = \frac{\sigma^2}{n}\left(\frac{(2-\theta)^2}{\kappa + \theta - 2} + \frac{2}{\theta}\right)\sum_{m\in\mathcal M_n}e^{-L_m|m|}. $$ Based on this extremely nuanced result, BM07 proposed the following choice of the penalty parameters:
where $H_m$ is the number of models with $|m|$ basis functions.\footnote{See A19 for a review of methods to estimate $\sigma^2$ in the model selection framework, which can be used to make the choice (ref) feasible.}
Indeed, setting $\theta = 1/2$, $\delta>0$ such that $3\sqrt \delta + 4\delta <1/2$, $\kappa>3/2$ such that $\kappa + 3\sqrt\delta + 4\delta < 2$, and $L_m = (\log H_m)/|m| + \delta$ for all $m\in\mathcal M_n$, we have that all inequalities in (ref) with $\lambda_m$'s given by (ref) are satisfied, and so it follows from (ref) that $$ \mathbb{E}\left[\|\widetilde f - f\|^2_{2,n}\right] \leq 2 \inf_{m\in\mathcal M_n} \left\{\mathbb{E}\left[\|\widehat f_m - f\|^2_{2,n}\right] + \frac{8\sigma^2|m|\log N_n}{n} \right\} + \frac{c\sigma^2}{n} $$ for the complete variable selection and
for the ordered variable selection, where $c>0$ is some constant. These are oracle inequalities of the form (ref) in the prediction norm.
Note also that the choice (ref) yields the Mallows method for the ordered variable selection under homoskedasticity; compare (ref) with $H_m = 1$ for all $m\in\mathcal M_n$ and (ref). Moreover, the result (ref) is sharp enough to recover the asymptotic optimality of the Mallows method in this case under the assumption that $\min_{m\in\mathcal M_n}n\mathbb{E}[\|\widehat f_m - f\|_{2,n}^2]\to \infty$. Indeed, setting $\theta = \theta_n$, $\kappa = \kappa_n = 2 - \theta_n/2$, $\delta = \delta_n$ such that $3\sqrt{\delta_n} + 4\delta_n = \theta_n / 4$, and $L_m = (\log H_m)/|m| + \delta_n$ for all $m\in\mathcal M_n$, it follows that as long as $\theta_n\to0$ sufficiently slowly, we have $\varepsilon_n / \min_{m\in\mathcal M_n}\mathbb{E}[\|\widehat f_m - f\|_{2,n}^2] \to 0$, yielding the asymptotic optimality result.
More generally, the penalization approach is popular in the empirical risk minimization framework, where the function $f$ is assumed to solve the optimization problem $f = \arg\min_{f\in\mathcal F}\mathbb{E}[m(Z,f)]$, and our task is to choose an estimator in the set $\{\widehat f_k = \arg\min_{f\in\mathcal F_k}n^{-1}\sum_{i=1}^n m(Z_i,f)\colon k\in\mathcal K_n\}$, where $\mathcal F_k$'s are subsets of $\mathcal F$. See M07 and K11 for detailed expositions of this framework.
Instead of selecting one estimator $\widehat f_{\widehat k}$ from the collection of estimators $\{\widehat f_k\colon k\in\mathcal K_n\}$, the aggregation method takes a weighted average of these estimators: $\widetilde f(\cdot) = \sum_{k\in\mathcal K_n} w_k \widehat f_k(\cdot)$ for some weights $w_k$. The idea here is that by taking a weighted average of the estimators, we can average out their estimation errors in the same fashion as diversifying portfolio helps to minimize its risk in finance. In this section, I will describe a weighting scheme due to LB06 that yields an estimator $\widetilde f$ with particularly attractive properties. Throughout this section, we will assume that the model is homoskedastic and Gaussian, i.e. $ e \mid X = N(0,\sigma^2)$.
To describe their weighting scheme, recall that as follows from our discussion in Section (ref), an unbiased estimator of the risk $\mathbb{E}[\|\widehat f_k - f\|_{2,n}^2]$ of the series estimator $\widehat f_k$ based on $k$ terms can be taken to be $$ \widehat r_k = \frac{1}{n}\sum_{i=1}^n (Y_i - \widehat f_k(X_i))^2 + \frac{2\sigma^2 k}{n} - \sigma^2, \quad k\in\mathcal K_n. $$ It is thus natural to give higher weights $w_k$ to estimators $\widehat f_k$ with smaller $\widehat r_k$. Indeed, LB06 proposed setting
Note that the resulting estimator $\widetilde f^{LB}(\cdot) = \sum_{k\in\mathcal K_n}w^{LB}_k\widehat f_k(\cdot)$ is highly non-linear in $Y_i$'s but using the Stein method as described in Section (ref), they showed that an unbiased estimator $\widetilde r^{LB}$ of the risk of $\widetilde f^{LB}$ can be taken to be $\widetilde r^{LB} = \sum_{k\in\mathcal K_n} w_k^{LB}\widehat r_k,$ and this result can be used to derive an oracle inequality for $\widetilde f^{LB}$. Indeed, solving (ref) for $\widehat r_k$, we have
Also, denoting $k^* = \arg\min_{k\in\mathcal K_n} \widehat r_k$ and subtracting the expression (ref) for $k = k^*$ from the same expression for general $k$, $$ \widehat r_k - \widehat r_{k^*} = - \frac{4\sigma^2}{n}\left( \log w_k^{LB} - \log w_{k^*}^{LB} \right) \leq - \frac{4\sigma^2}{n}\log w_k^{LB}. $$ Combining this inequality with $\widetilde r^{LB} = \sum_{k\in\mathcal K_n} w_k^{LB}\widehat r_k$, we obtain $$ \widetilde r^{LB} \leq \widehat r_{k^*} - \frac{4\sigma^2}{n}\sum_{k\in\mathcal K_n} w_k^{LB}\log w_k^{LB} \leq \widehat r_{k^*} + \frac{4\sigma^2\log|\mathcal K_n|}{n}. $$ Finally, taking expectations on both sides of this inequality, we obtain the main result from LB06: if $e\mid X\sim N(0,\sigma^2)$, then
This result is remarkable because it yields an oracle inequality of the type (ref) with the constant $C=1$. Compare, for example, this inequality with the oracle inequality (ref) for the estimator $\widehat f_{\widehat k}$ based on the penalization method, where we had the constant $C=2$. In fact, the right-hand of (ref) is typically smaller than the right-hand of (ref). In particular, as long as $|\mathcal K_n|\leq n$, the bound (ref) is preferred as long as $k^* = \arg\min_{k\in\mathcal K_n} \mathbb{E}[\|\widehat f_k - f\|_{2,n}^2]$ is of larger order than $\log n$. Related results can be found, for example, in C97, Y00, Y04, and DT07.
In this section, I compare the performance of the methods described above via simulations. Since the Stein and penalization methods coincide in the case of the series estimator with the Mallows method, I omit them from consideration. I thus focus on the following four methods: Mallows, Lepski, cross-validation, and aggregation. For all methods, I consider the nonparametric mean regression model (ref) with the function $f$ being given by either $f(x) = \exp(\exp(x))$ or $f(x)=\sin(2\pi x)$ for all $x\in\mathbb{R}$. I assume that $X\sim U([0,1])$ and $e=\varepsilon/\sqrt{1+X^2}$, where $\varepsilon\mid X\sim N(0,1)$. Further, I consider random samples of size $n = 500$ and $n=1000$. I study the series estimator using either monomials or quadratic splines as described above. For both sets of basis functions, I set $\mathcal K_n = \{1,\dots,[n^{1/3}]\}$. For the Mallows, Lepski, and aggregation methods, I obtain regression residuals $\widehat e_i$ using $\widehat f_k$ with $k=[n^{1/3}]$. For the Lepski method, I set $x_0 = 0.5$. For the cross-validation method, I focus on the $V$-fold cross-validation with $V=5$. Finally, since the aggregation method as described in the previous section is aiming at estimating homoskesdatic model, I extend it in a somewhat ad-hoc way as follows. I set $$ w_k = \frac{\exp(-n\widehat r_k/(4\widehat\sigma^2))}{\sum_{k'\in\mathcal K_n}\exp(-n\widehat r_{k'}/(4\widehat\sigma^2))},\quad\text{for all }k\in\mathcal K_n, $$ where $\widehat \sigma^2 = n^{-1}\sum_{i=1}^n \widehat e_i^2$ and $$ \widehat r_k = \frac{1}{n}\sum_{i=1}^n(Y_i - \widehat f_k(X_i))^2 + \frac{2\widehat A_k}{n} - \widehat\sigma^2,\quad\text{for all }k\in\mathcal K_n, $$ is an unbiased estimator of the risk $\|\widehat f_k - f\|_{2,n}^2$ of the series estimator $\widehat f_k$ appearing in our discussion of the Mallows method; see in particular (ref) for the definition of $\widehat A_k$. All other specifications coincide with recommendations for individual methods above.
The simulation results are presented in Table 1. Overall, we observe that estimation results improve as we increase the sample size. Also, the aggregation method tends to give the best results. The Lepski method may fail to give good results in the $\ell_2$, prediction, and uniform metrics, see the case of monomials and $f(x) = \sin(2\pi x)$. This is not a bug but a feature of the method, as it is constructed specifically to yield good results in the pointwise metric, and it does perform very well in the pointwise metric. In fact, in the case of monomials and $f(x) = \sin(2\pi x)$, the Lepski method substantially outperforms the other methods in the pointwise metric. The Mallows and cross-validation methods have similar performance.
Table 1: This table shows the average estimation error in the $\ell_2$ metric ($\ell_2$), prediction metric ($\ell_{2,n}$), uniform metric ($\ell_{\infty}$), and the pointwise metric ($\ell_{pw}$) of the series estimator based on either monomials (M) or quadratic splines (S) with the number of series terms determined by the Mallows, Lepski, CV, and Aggregation methods. CV means 5-fold cross-validation, and the three metrics are given by (ref), (ref), and (ref) with the supremum in (ref) being taken over $\mathcal X = [0,1]$ and $x_0$ in (ref) being equal to 0.5.
In this section, I discuss tuning parameter selection in the context of high-dimensional estimation via $\ell_1$-penalization. I discuss methods based on the theory of self-normalized moderate deviations, bootstrap, Stein's unbiased risk estimation, and cross-validation. I mainly focus on the case of i.i.d. data but describe extensions for the cases of clustered and panel data as well. Similarly, I mainly focus on the case of the mean regression but also touch the cases of the quantile regression and, more generally, generalized linear models.
Since the literature on high-dimensional estimation in general and on tuning parameter selection in high-dimensional settings in particular is large, my review here is rather selective. For example, I do not discuss methods such as the square-root Lasso, where the choice of the tuning parameter is relatively simple. An interested reader can find this and other existing methods, for example, in WW20.
Consider the following mean regression model:
where $Y\in\mathbb{R}$ is a dependent variable, $X\in\mathbb R^p$ is a vector of covariates, $\beta\in\mathbb R^p$ is a vector of parameters, and $e\in\mathbb R$ is noise. Assume that we have a random sample $(X_1,Y_1),\dots,(X_n,Y_n)$ from the distribution of the pair $(X,Y)$. Assume also that the model is high-dimensional in the sense that $p$ is large, potentially larger or much larger than $n$, but sparse in the sense that the number of non-zero components of the vector $\beta$, say $s$, is small.
In the context of this model, it is well-known that the OLS estimator is rather poor and that the Lasso estimator is preferred:
where $\lambda> 0$ is a penalty parameter. In this section, I review the basic Lasso theory, as it guides the choices of $\lambda$ to be discussed below.\footnote{The discussion throughout Section (ref) assumes that the random variables in the vector $X$ are normalized so that their second moments are of the same order. If they are not normalized, it may be appropriate to consider the Lasso estimator of the form $\widehat\beta = \arg\min_{b\in\mathbb{R}^p}\{n^{-1}\sum_{i=1}^n (Y_i - X_i^{\top}b)^2 + \lambda \sum_{j=1}^p \kappa_j |b_j|\}$, where $\kappa_j$'s are appropriate loadings, e.g. see BCCH12.}
Let $T = \{j=1,\dots,p\colon \beta_j\neq 0\}$ be the set of indices corresponding to non-zero components of the vector $\beta = (\beta_1,\dots,\beta_p)^{\top}$ and let $T^c = \{1,\dots,p\}\setminus T$ be the set of remaining indices. Also, for any vector $\delta = (\delta_1,\dots,\delta_p)^{\top}\in\mathbb R^p$, let $\delta_T = (\delta_{T1},\dots,\delta_{Tp})^{\top}\in\mathbb R^p$ be the vector defined by $\delta_{Tj} = \delta_j$ if $j\in T$ and $\delta_{Tj} = 0$ if $j\in T^c$ and let $\delta_{T^c}\in\mathbb R^p$ be the vector defined in the same way where the roles of $T$ and $T^c$ are switched. In addition, let $ \|\delta\|_{2,n} = \sqrt{n^{-1}\sum_{i=1}^n (X_i^{\top}\delta)^2} $ denote the prediction norm of the vector $\delta\in\mathbb{R}^p$. Moreover, for any $c>1$, let $\mathcal R_c = \{\delta\in \mathbb{R}^p\colon \|\delta_{T^c}\|_1 \leq c\|\delta_T\|_1\}$ be the so-called {\em restricted set} and let $$ \kappa_c = \inf_{\delta\in \mathcal R_c}\frac{\sqrt s\|\delta\|_{2,n}}{\|\delta_T\|_1} $$ be the so-called {\em compatibility constant}. As we will see below, the compatibility constant plays an important role in the analysis of the Lasso estimator. Importantly, however, for any $c>1$, the compatibility constant $\kappa_{c}$ is bounded below from zero under so-called {\em sparse eigenvalue conditions}, which can be intuitively understood as requiring that there is no in-sample linear dependence between any $s\log n$ variables in the vector $X$; see BRT09 and BC11 for more precise descriptions.
Using the concept of the compatibility constants, we have the following result: For any $c>1$, if
then
where $\bar c = (c+1)/(c-1)$; see BRT09 and BC11. A bound that is proportional to the bound for $\|\widehat\beta - \beta\|_{2,n}$ can also be derived, under certain conditions, for $\|\widehat \beta - \beta\|_2$; e.g. see BC10. This result suggests the following general principle for choosing the penalty parameter $\lambda$: Choose $\lambda$ as small as possible subject to the constraint that (ref) holds with probability close to one. For example, under the assumption of homoskedastic and Gaussian noise, i.e. if $e\mid X\sim N(0,\sigma^2)$ with some known $\sigma>0$, setting $\lambda$ to be equal to
for some $c>1$ and $\alpha\in(0,1)$, where $\Phi$ is the cdf of the $N(0,1)$ distribution, will ensure via the union bound that (ref) holds with probability at least $1 - \alpha$. Thus, we can use the choice (ref) with some $\alpha$ close to zero and some $c>1$, which I refer to as the BRT method. In practice, it is recommended to set $\alpha = 0.1/\log(p\vee n)$ and $c = 1.1$, where $p\vee n$ denotes the maximum between $p$ and $n$; e.g. see BCCH12.\footnote{The choice of the constant $c$ is actually rather tricky. Indeed, the formal result (ref)-(ref) requires that the constant $c$ satisfies $c>1$. However, it is well-known via simulations that optimal performance of the Lasso estimator based on the penalty parameter $\lambda$ selected according to the BRT method (ref) is often achieved by values of $c$ satisfying $c<1$. This, on the one hand, motivates the choice $c=1.1$ but, on the other hand, also implies that it would be of great interest to better understand the behavior of the Lasso estimator when the constraint (ref) is not satisfied.} In the next sections, I will describe methods that allow to choose $\lambda$ without imposing the assumption of homoskedastic and Gaussian noise.
To conclude this section, I note that under mild regularity conditions, the random variable on the right-hand side of (ref) converges to zero with the rate $\sqrt{\log p /n}$. Thus, given that the compatibility constant $\kappa_{\bar c}$ is bounded away from zero under sparse eigenvalue conditions, it follows that, with a proper choice of the penalty parameter $\lambda$, under certain regularity conditions, the Lasso estimator satisfies
Note also that a similar result can be obtained for the approximately sparse model, where the vector $\beta$ does not have to be sparse but can be well approximated by a sparse vector with $s$ non-zero components.
I now discuss how to choose the penalty parameter $\lambda$ for the Lasso estimator following the general principle described in the previous section but without assuming homoskedastic and Gaussian noise. In this section, I focus on the method based on the theory of self-normalized moderate deviations. In the next section, I will describe the method based on the bootstrap theory.
The method based on the theory of self-normalized moderate deviations was proposed in BCCH12. This theory states that even if the random variables $X_{ij}e_i$ have only three finite moments, the ratio of the distribution function of the random variables $$ T_j = \frac{\sum_{i=1}^n X_{ij} e_i}{\sqrt{\sum_{i=1}^n X_{ij}^2 e_i^2}} $$ and the distribution function of $N(0,1)$ random variables converges to one not only pointwise, which would follow from the classic Central Limit Theorems, but also uniformly over sets $(-\bar x_n,\bar x_n)$ as long as $\bar x_n = o(n^{1/6})$:
e.g. see JSW03 for the original result or Lemma A.1 in BCCHKN18 for exposition.\footnote{The classic Berry-Esseen theorem leads to bounds on $|\mathrm{P}(T_j > x) - P(N(0,1)>x)|$. The result (ref) is much stronger, however, as it provides a bound on $|\mathrm{P}(T_j > x) - P(N(0,1)>x)|/P(N(0,1)>x)$ and $\mathrm{P}(N(0,1)>x)\to 0$ as $x\to\infty$.} Thus, by denoting $$ \bar\lambda_n = \frac{2c}{\sqrt n}\Phi^{-1}\left(1 - \frac{\alpha}{2p}\right)\max_{1\leq j\leq p}\sqrt{\frac{1}{n}\sum_{i=1}^n X_{ij}^2 e_i^2}, $$ it follows via the union bound that
as long as $\log(p/\alpha) = o(n^{1/3})$ since $\Phi^{-1}(1 - \alpha/(2p)) \leq \sqrt{2\log(2p/\alpha)}$ by Proposition 2.5 in D14. Note, however, that setting $\lambda = \bar\lambda_n$ is infeasible since we do not observe $e_i$'s. BCCH12 proposed the following three-step algorithm, which I refer to as the BCCH method, to work around this issue. First, we run Lasso with
to get a preliminary Lasso estimator $\widehat \beta^P$. Second, we calculate regression residuals $\widehat e_i = Y_i - X_i^{\top}\widehat\beta^P$ for all $i=1,\dots,n$. Third, we set $\lambda$ to be equal to
The final Lasso estimator is then $\widehat\beta = \widehat\beta(\widehat\lambda^{BCCH})$. The idea here is that because $\sum_{i=1}^n \mathbb{E}[X_{ij}^2Y_i^2] > \sum_{i=1}^n \mathbb{E}[X_{ij}^2 e_i^2]$ for all $j=1,\dots,p$, which follows from (ref), we have via the uniform law of large numbers that under certain regularity conditions,
with probability approaching one. Thus, the preliminary Lasso estimator $\widehat\beta^P$ satisfies the constraint (ref) with probability at least $1 - \alpha + o(1)$ and, as long as $\alpha=\alpha_n \to 0$, yields consistent estimators $\widehat e_i$ of $e_i$'s. On the other hand, the preliminary Lasso estimator is conservative: because of (ref), its penalty parameter $\lambda$ is larger than necessary. This issue is fixed on the third step, where we replace $Y_i$'s by $\widehat e_i$'s in the penalty choice $\lambda$; compare (ref) and (ref).
Similarly to our discussion of the BRT method, for both the preliminary and the final Lasso estimators, BCCH12 suggested setting $\alpha = 0.1/\log(p\vee n)$ and $c = 1.1$. They proved that under certain regularity conditions, including the sparse eigenvalue conditions and the growth condition $\log p = o(n^{1/3})$, which appeared above because of the requirements of the theory of self-normalized moderate deviations, the final Lasso estimator $\widehat \beta = \widehat\beta(\widehat\lambda^{BCCH})$ satisfies (ref). KPS24 extended this result to the case of a high-dimensional VAR model.
The BCCH method described in the previous section is still somewhat conservative: the third step of the BCCH method eliminates the slackness arising from (ref) but there is still slackness arising from our use of the union bound in (ref). Indeed, as long as the variables in the vector $X$ are correlated, a value of the penalty parameter $\lambda$ that is smaller than that produced by the BCCH method may be sufficient to satisfy constraint (ref) with probability close to one. BC10 illustrated this point via simulations. In turn, smaller values of $\lambda$ will lead to tighter bounds on the estimation error in (ref) and are thus preferable. In this section, I describe a bootstrap method that gives $\lambda$ such that the probability of (ref) is {\em equal} to $1-\alpha + o(1)$.
The method consists of four steps. First and second steps are exactly the same as in the previous section: we first run Lasso with $\lambda$ given in (ref) to obtain the preliminary Lasso estimator $\widehat \beta^P$ and then calculate regression residuals $\widehat e_i = Y_i - X_i^{\top}\widehat\beta^P$ for all $i=1,\dots,n$. On the third step, we estimate the $1-\alpha$ quantile of the random variable on the right-hand side of (ref) via multiplier bootstrap. Specifically, let $w_1,\dots,w_n$ be i.i.d. $N(0,1)$ random variables that are independent of the data and let $\widehat q_{1-\alpha}$ be the $1-\alpha$ quantile of the conditional distribution of
given the data $\{(X_i,Y_i)\}_{i=1}^n$, and thus given $\{(X_i,\widehat e_i)_{i=1}^n\}$. The estimator of the $1-\alpha$ quantile of the random variable on the right-hand side of (ref) is then $\widehat q_{1-\alpha}/\sqrt n$. On the fourth step, we set $\lambda$ to be equal to $\widehat\lambda^{B} = (2c/\sqrt n)\widehat q_{1-\alpha}$, so that the final Lasso estimator is $\widehat \beta = \widehat\beta(\widehat\lambda^B)$, where we again choose $\alpha = 0.1/\log(p\vee n)$ and $c=1.1$.
The idea behind the bootstrap method is that conditional on $\{(X_i,\widehat e_i)_{i=1}^n\}$, the distribution of the vector $n^{-1/2}\sum_{i=1}^n w_i X_i \widehat e_i$ is centered Gaussian with covariance matrix $n^{-1}\sum_{i=1}^n \widehat e_i^2 X_iX_i^{\top}$, and so, using the Gaussian comparison inequalities, as developed in CCK15, it is possible to show that the distribution of the max norm (ref) of this vector can be well approximated by the distribution of the max norm of the centered Gaussian vector with covariance matrix $\mathbb{E}[e^2XX^{\top}]$. In turn, it follows from the high-dimensional Central Limit Theorems as developed in CCK13, CCK17, CCKK22 that the latter distribution also provides a good approximation to the distribution of the max norm of the vector $n^{-1/2}\sum_{i=1}^n X_{i}e_i$, e.g. see BCCHKN18 for an exposition of these results. Formally, one can show that using $\lambda = \widehat\lambda^B = (2c/\sqrt n)\widehat q_{1-\alpha}$ ensures that (ref) holds with probability $1 - \alpha + o(1)$ as long as the growth condition $\log p = o(n^{1/5})$ and other regularity conditions are satisfied, which in turn implies that the Lasso estimator $\widehat\beta = \widehat\beta(\widehat\lambda^B)$ satisfies (ref).
Another approach for choosing the penalty parameter $\lambda$ for the Lasso estimator (ref) is based on the Stein method, as developed in ZHT07 and TT12. The method takes a particularly simple form if we assume that the distribution of the vector $X$ is absolutely continuous with respect to the Lebesgue measure on $\mathbb{R}^p$, as the Lasso estimator is almost surely unique in this case. Assuming also that $e\mid X\sim N(0,\sigma^2)$, their method consists of setting $\lambda$ to be equal to
where we use $\|\widehat\beta(\lambda)\|_0$ to denote the number of non-zero components of the vector $\widehat\beta(\lambda)$. As proven in ZHT07 for the low-dimensional case and in TT12 for the high-dimensional case, the criterion function in (ref) is an unbiased estimator of the risk of the Lasso estimator $\|\widehat\beta(\lambda) - \beta\|_{2,n}^2$, where we condition on $X_1,\dots,X_n$.\footnote{This method can also be referred to as the selection via degrees of freedom, as $\|\widehat\beta(\widetilde\lambda)\|_0$ turns out to be an unbiased estimator of the degrees of freedom of the Lasso estimator $\widehat\beta(\widehat\lambda)$; see ZHT07 and TT12. Also, as proven in ZHT07, the minimum in (ref) is achieved by one of the threshold values of $\lambda$, which are values such that as we increase $\lambda$, crossing them will change the set of non-zero components of the estimator $\widehat\beta(\lambda)$. The latter observation makes the problem of solving (ref) feasible.}
Intuitively, the KKT conditions for the optimization problem in (ref) are $$ \frac{2}{n}\sum_{i=1}^n X_{ij}(Y_i - X_i^{\top}\widehat\beta(\lambda)) = \lambda\gamma_j,\quad\text{for all }j=1,\dots,p, $$ where $\gamma = (\gamma_1,\dots,\gamma_p)^{\top}$ is a vector such that for all $j=1,\dots,p$, we have $\gamma_j = 1$ if $\widehat\beta_j(\lambda) > 0$, $\gamma_j = -1$ if $\widehat\beta_j(\lambda)<0$, and $\gamma_j\in[-1,1]$ if $\widehat\beta_j(\lambda)=0$. Therefore, denoting $\mathrm{Y} = (Y_1,\dots,Y_n)^{\top}$ and $\widehat T = \{j=1,\dots,p\colon \widehat\beta_j(\lambda)\neq 0\}$ and letting $\mathrm{X}_{\widehat T}$ be the matrix consisting of the columns of $\mathrm{X} = (X_1,\dots,X_n)^{\top}$ corresponding to indices in $\widehat T$ and $\gamma_{\widehat T}$ be the vector consisting of the elements of $\gamma$ corresponding to indices in $\widehat T$, it follows that the Lasso estimator of the vector $\mathrm{f} = (X_1^{\top}\beta,\dots,X_n^{\top}\beta)^{\top}$ is $$ \widehat\mathrm{f} = \mathrm{X}\widehat\beta(\lambda) = \mathrm{X}_{\widehat T}(\mathrm{X}_{\widehat T^{\top}}\mathrm{X}_{\widehat T})^{-1}(\mathrm{X}_{\widehat T}\mathrm{Y} - n\lambda \gamma_{\widehat T}/2). $$ Hence, expression (ref) appears from expression (ref) by taking the derivatives of $\widehat\mathrm{f}$ with respect to $\mathrm{Y}$ as discussed in Section (ref) and noting that both $\widehat T$ and $\gamma_{\widehat T}$ are locally independent of $\mathrm{Y}$ almost surely (in the sense that for almost all $\mathrm{Y}$, small perturbations of $\mathrm{Y}$ do not lead to any changes in $\widehat T$ and $\gamma_{\widehat T}$).
The Lasso estimator $\widehat\beta = \widehat\beta(\widehat\lambda^S)$ is potentially attractive but its main disadvantage is that the choice (ref) is only justified under the assumption of homoskedastic and Gaussian noise. In addition, relative to the Mallows method of choosing the number of series terms, the analysis here requires studying concentration properties of the random variables $\|\widehat\beta(\lambda)\|_0$, which is difficult; compare (ref) and (ref). BZ21 derived a bound for the variance of $\|\widehat\beta(\lambda)\|_0$ but, to the best of our knowledge, it is still not clear if the Lasso estimator $\widehat\beta = \widehat\beta(\widehat\lambda^S)$ satisfies (ref).
The penalty parameter $\lambda$ for the Lasso estimator (ref) can also be selected by $V$-fold cross-validation.\footnote{In principle, I could consider here the leave-one-out cross-validation as well. However, even though the Lasso estimator is rather fast to compute, the leave-one-out cross-validation requires computing the Lasso estimator $n$ times, which does typically take too much time.} The algorithm proceeds as follows. Let $\Lambda_n$ be a set of candidate values of $\lambda$. This set is assumed to be large enough to contain a good choice of $\lambda$. In practice, this set can be chosen to consist of a geometric progression starting from some small value, say $1/n$, and ending at some large value, say $n$. Further, let $(I_1,\dots,I_V)$ be a random partition of the whole sample $I_0 = \{1,\dots,n\}$ into $V$ subsamples of roughly the same size. For each $\lambda \in \Lambda_n$ and $v=1,\dots,V$, we build the Lasso estimator of the vector $\beta$ using all the data excluding observations in subsample $v$: $$ \widehat\beta(\lambda,v) = \arg\min_{b\in\mathbb{R}^p}\left\{\frac{1}{|I_0\setminus I_v|}\sum_{i\in I_0\setminus I_v}(Y_i - X_i^{\top}b)^2 + \lambda\|b\|_1\right\}. $$ We then choose the penalty parameter $\lambda$ to be equal to a minimizer of the sum of out-of-sample prediction errors of these estimators:
The final Lasso estimator can be either the average estimator $$ \widehat \beta^{AVCV} = \frac{1}{V}\sum_{v=1}^V \widehat\beta(\widehat\lambda^{VCV},v) $$ or the full-sample estimator $\widehat\beta^{FVCV} = \widehat \beta(\widehat\lambda^{VCV})$, with the estimators $\widehat\beta(\lambda)$ being defined in (ref).
LM12 derived an oracle inequality for the average estimator $\widehat\beta^{AVCV}$. Intuitively, the definition (ref) implies that for any $\lambda\in\Lambda_n$,
Also, by the uniform law of large numbers and independence of $\{(X_i,Y_i)\}_{i\in I_v}$ from $\widehat\beta(\lambda,v)$, under certain regularity conditions,
uniformly over $\lambda\in\Lambda_n$, where $\mathbb{E}_{X,Y}[\cdot]$ denotes expectation with respect to the pair $(X,Y)$ holding the data fixed. Thus, combining (ref) and (ref), noting that $\mathbb{E}_{X,Y}[(Y - X^{\top}\widehat\beta(\lambda,v))^2] = \mathbb{E}[e^2] + E_{X,Y}[(X^{\top}(\beta - \widehat\beta(\lambda,v)))^2]$, and assuming for simplicity that $|I_1|=\dots=|I_V|$, we obtain $$ \sum_{v=1}^V \mathbb{E}_{X,Y}[(X^{\top}(\widehat\beta(\widehat\lambda^{VCV},v) - \beta)^2] \leq \sum_{v=1}^V \mathbb{E}_{X,Y}[(X^{\top}(\widehat\beta(\lambda,v) - \beta)^2] + O_P\left(\sqrt{\frac{\log|\Lambda_n|}{n}}\right). $$ uniformly over $\lambda\in\Lambda_n$. Here, all $V$ terms under the sum on the right-hand side are the same under our i.i.d. data assumption. Thus, by the triangle and Cauchy-Schwarz inequalities,
Moreover, using a more careful argument, LM12 were able to get rid of the square root in the $O_P(\cdot)$ term in expense for having some extra terms in the inequality. Their oracle inequality implies in particular that, under certain regularity conditions, the average estimator $\widehat\beta = \widehat\beta^{AVCV}$ satisfies the $\ell_2$ norm bound in (ref).
The analysis of the full-sample estimator $\widehat\beta^{FVCV}$ is more complicated. The key issue here is that the basic inequality (ref) does not provide any direct information regarding this estimator. Moreover, simulation evidence suggests that the penalty parameter value selected by cross-validation, $\lambda = \widehat\lambda^{VCV}$, typically does not satisfy the constraint (ref), making the basic Lasso theory described in Section (ref) inapplicable. To deal with these issues, CLC21 derived a relationship between the estimation error of the full-sample Lasso estimators $\widehat\beta(\lambda)$ and the corresponding subsample Lasso estimators $\widehat\beta(\lambda,v)$. The relationship turns out to depend on the estimation errors of the subsample Lasso estimators in the $\ell_1$ norm $\|\widehat\beta(\lambda,v) - \beta\|_1$, which can not be bounded using the Lecue and Mitchell arguments directly. Instead, CLC21 derived bounds on $\|\widehat\beta(\lambda,v) - \beta\|_1$ using the degrees of freedom formula for the Lasso estimator derived in ZHT07 and TT12. As a result, they proved that, under certain regularity conditions, the full-sample estimator $\widehat\beta = \widehat\beta^{FVCV}$ satisfies
These bounds, however, are by some $\log p$ factors worse than those in (ref), which are guaranteed for the Lasso estimator with the penalty parameter selected either by the BCCH method or by the bootstrap method. It is therefore not clear whether the bounds (ref) are sharp.
In this section, I point out that one should be careful with methods described above if the data is clustered. To make this point clearly, suppose that $n$ observations $\{1,\dots,n\}$ can be split into $G$ groups such that there are $L$ observations within each group, so that $n=GL$, and all observations within the same group are perfectly correlated, so that $X_{ij}e_i = X_{lj}e_l$ whenever observations $i$ and $l$ belong to the same group. Then, letting $I(g)$ be a function that selects one observation from each group $g=1,\dots,G$, it follows that the constraint (ref) reduces to $$ \lambda \geq 2c\max_{1\leq j\leq p}\left| \frac{1}{G}\sum_{g=1}^G X_{I(g) j}e_{I(g)} \right|. $$ Hence, for example, the BRT rule, which applies when $e\mid X\sim N(0,\sigma^2)$ with known $\sigma$, yields
This expression should be compared with that in (ref), which would be used by the researcher who ignores the clustering in the data. In particular, the correct expression (ref) is $\sqrt L$ times larger, and using the wrong expression (ref) would yield the value of the penalty parameter $\lambda$ that does not necessarily satisfies (ref) with probability at least $1-\alpha$. As a result, the researcher ignoring clustering could obtain misleading estimates that are far away from the true vector $\beta$.
Fortunately, it is possible to fix this issue as long as the number of clusters/groups is large. For concreteness, I focus here on the BCCH method, as studied in C18. To adjust this method so that it works with clustered data, suppose as above that $n$ observations $\{1,\dots,n\}$ can be split into $G$ groups, that observations within each group can be arbitrarily correlated, and that observations across groups are independent. Letting $N(g)$ denote all observations $i$ within each group $g=1,\dots,G$, we then consider the following cluster-robust version of the BCCH method. First, we run Lasso with $$ \lambda = \frac{2c}{\sqrt n}\Phi^{-1}\left(1 - \frac{\alpha}{2p}\right)\max_{1\leq j\leq p}\sqrt{\frac{1}{n}\sum_{g=1}^G\left(\sum_{i\in N(g)} X_{ij} Y_i\right)^2} $$ to get a preliminary Lasso estimator $\widehat \beta^P$. Second, we calculate regression residuals $\widehat e_i = Y_i - X_i^{\top}\widehat\beta^P$ for all $i=1,\dots,n$. Third, we set $\lambda$ to be equal to $$ \widehat\lambda^{BCCH,C} = \frac{2c}{\sqrt n}\Phi^{-1}\left(1 - \frac{\alpha}{2p}\right)\max_{1\leq j\leq p}\sqrt{\frac{1}{n}\sum_{g=1}^G\left(\sum_{i\in N(g)} X_{ij} \widehat e_i\right)^2} $$ The final Lasso estimator is then $\widehat\beta = \widehat\beta(\widehat\lambda^{BCCH,C})$. Using the same ideas as those in BCCH12, it is then straightforward to show that this estimator satisfies (ref) with $n$ replaced by $n/L$, where $L = \max_{1\leq g\leq G}|N(g)|$ denotes the maximal number of observations within groups.
In this section, I consider the Lasso estimator in the context of panel data, where the dependent variable $Y_{it}\in \mathbb{R}$ for unit $i$ at time period $t$ satisfies the model $$ Y_{it} = X_{it}^{\top}\beta + \xi_i + e_{it},\quad i=1,\dots,N, \ t=1,\dots,T, $$ where $X_{it}\in\mathbb{R}^p$ is a vector of covariates, $\beta\in\mathbb R^p$ is a vector of parameters, $\alpha_i\in\mathbb{R}$ is an individual fixed effect that can be correlated with covariates $X_{it}$'s, and $e_{it}\in\mathbb{R}$ is noise that satisfies $\mathbb{E}[e_{it}\mid X_{i1},\dots,X_{iT}] = 0$. I assume that observations are independent across units and that the total number of observable periods $T$ is small.
The Lasso version of the classic fixed-effect estimator in this model takes the following form:
where $\widetilde Y_{it} = Y_{it} - T^{-1}\sum_{s=1}^T Y_{is}$ and $\widetilde X_{it} = X_{it} - T^{-1}\sum_{s=1}^T X_{is}$ for all $i=1,\dots,N$ and $t=1,\dots,T$. Using our discussion in Section (ref), it is then clear that a good choice of the penalty parameter $\lambda$ in (ref) should yield
with probability close to one, where $\widetilde e_{it} = e_{it} - T^{-1}\sum_{s=1}^T e_{is}$ for all $i=1,\dots,N$ and $t=1,\dots,T$; compare (ref) with (ref). In turn, using independence of observations across units and our discussion of clustered data in the previous section, it is straightforward to show that the panel version, for example, of the BCCH method takes the following form. First, we run Lasso with $$ \lambda = \frac{2c}{\sqrt{NT}}\Phi^{-1}\left(1 - \frac{\alpha}{2p}\right)\max_{1\leq j\leq p}\sqrt{\frac{1}{NT}\sum_{i=1}^N\left(\sum_{t=1}^T \widetilde X_{itj} \widetilde Y_{it}\right)^2} $$ to get a preliminary Lasso estimator $\widehat \beta^P$. Second, we calculate regression residuals $\widehat e_{it} = \widetilde Y_{it} - \widetilde X_{it}^{\top}\widehat\beta^P$ for all $i=1,\dots,N$ and $t=1,\dots,T$. Third, we set $\lambda$ to be equal to $$ \widehat\lambda^{BCCH,P} = \frac{2c}{\sqrt{NT}}\Phi^{-1}\left(1 - \frac{\alpha}{2p}\right)\max_{1\leq j\leq p}\sqrt{\frac{1}{NT}\sum_{i=1}^N\left(\sum_{t=1}^T \widetilde X_{itj} \widehat e_{it}\right)^2}. $$ The final Lasso estimator is then $\widehat\beta = \widehat\beta(\widehat\lambda^{BCCH,P})$. Using the same arguments as those in BCCH12, one can show that this estimator satisfies (ref) with $n$ replaced by $N$.
In this section, I briefly discuss estimation of the high-dimensional quantile regression model $$ Y = X^{\top}\beta + e,\quad Q_u(e\mid X) = 0, $$ where $Q_u(e\mid X)$ is the $u$th quantile of the conditional distribution of $e$ given $X$, which is assumed to be continuous, and $u\in(0,1)$ is a fixed quantile index. As before, I assume that $\beta\in\mathbb{R}^p$ is high-dimensional but sparse: the number of its non-zero components is equal to some small $s$. If $(X_1,Y_1),\dots,(X_n,Y_n)$ is a random sample from the distribution of the pair $(X,Y)$, the analog of Lasso for this quantile model is the following $\ell_1$-penalized quantile regression estimator: $$ \widehat \beta = \arg\min_{b\in\mathbb{R}^p} \left\{\frac{1}{n}\sum_{i=1}^n \rho_u(Y_i - X_i^{\top}b) + \lambda\|b\|_1\right\}, $$ where $t\mapsto \rho_u(t) = (u - 1\{t<0\})t$ is Koenker's check function and $1\{\cdot\}$ is the indicator function. BC11 proved the following result: if
then
where $C>0$ is a constant that is independent of $(n, p, s, \lambda)$, $u_i = F(e_i\mid X_i)$ for all $i=1,\dots,n$, and $F(t\mid X) = \mathrm{P}(e\leq t\mid X)$ is the cdf of the conditional distribution function of $e$ given $X$. This result in analogous to that in (ref)-(ref), so that the same general principle for selecting the penalty parameter $\lambda$ applies, but the key point here is that $u_i\mid X_i$ has the $U(0,1)$ distribution, and so the conditional distribution of the random variable on the right-hand side of (ref) given $X_i$'s is pivotal and can be simulated exactly. We can thus choose the penalty parameter $\lambda$ as $cq_{1-\alpha}/\sqrt n$ with $q_{1-\alpha}$ being the $1-\alpha$ quantile of the $\sqrt n$ times the random variable on the right-hand side of (ref), where $\alpha = 0.1/\log(p\vee n)$ and $c=1.1$ as before.
Finally, I consider the generalized linear model
where $m\colon \mathbb{R}\times\mathcal Y\to \mathbb{R}$ is a known function that is convex in its first argument and $\mathcal Y$ denotes the support of $Y$. Although both mean and quantile regression models fit into (ref), I am particularly interested here in other econometric models fitting into (ref), e.g. the logit model $$ \beta = \arg\min_{b\in\mathbb{R}^p}\mathbb{E}[\log(1+\exp(X^{\top}b)) - YX^{\top} b] $$ and the probit model $$ \beta = \arg\min_{b\in\mathbb{R}^p}\mathbb{E}[-Y\log\Phi(X^{\top}b) - (1-Y)\log(1-\Phi(X^{\top}b))]. $$ As before, assuming that $\beta$ is high-dimensional but sparse and that we have a random sample $(X_1,Y_1),\dots,(X_n,Y_n)$ from the distribution of the pair $(X,Y)$, the analogue of Lasso for this framework is the following $\ell_1$-penalized M-estimator:
Applying the same ideas as those in the Lasso theory, we can show that if
then the result of the form (ref) holds. We thus again want to choose the penalty parameter $\lambda$ as small as possible subject to the constraint that (ref) is satisfied with probability close to one. The main issue here, however, is that there may not exist a preliminary estimator allowing us to consistently estimate $\beta$ in the first step, as we did in the BCCH and bootstrap methods for the Lasso estimator. CS21 circumvented this problem and proposed obtaining the preliminary estimator using the subsample versions of the $\ell_1$-penalized estimator (ref) with the penalty parameter $\lambda$ selected by cross-validation, analogous to estimators $\widehat\beta(\widehat\lambda^{VCV},v)$ in Section (ref). The idea here is that even though full sample estimators with the penalty parameter selected by cross-validation are hard to analyze, as discussed in Section (ref), the corresponding subsample estimators yield the estimators of $m_1(X_i^{\top}\beta,Y_i)$ whose consistency is relatively easy to derive. CS21 proved that the $\ell_1$-penalized estimator (ref) with the penalty parameter selected by the multiplier bootstrap analogous to that described in Section (ref) with preliminary estimators of $m_1(X_i^{\top}\beta,Y_i)$ obtained via subsample versions of the estimator (ref) with the penalty parameter selected by cross-validation indeed satisfies a result of the form (ref).
Finally, I note that the result of CS21 can be extended to work with time series data as well: it suffices to replace cross-validation and bootstrap by large-block-small-block versions of both cross-validation and bootstrap. Moreover, this extension is useful not only for the generalized linear models (ref) but for the linear model (ref) as well. Indeed, it follows from the same arguments as those in Section (ref) that the choice of $\lambda$ in (ref) may not yield a consistent Lasso estimator in the case of time series data. In fact, as argued in KPS24, it may be rather difficult to provide a good choice of the penalty parameter $\lambda$ for the Lasso estimator in the case of time series data because the choice should depend on a measure of time-series dependence in the data. This complication can be avoided by using the large-block-small-block version of the method of CS21.