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.
88,890 characters · 22 sections · 0 citation commands
Phase transitions in nonparametric regressions
Estimation of an unknown univariate function $f$ from the nonparametric regression model
has been a central research topic in econometrics, machine learning, numerical analysis and statistics. Many semiparametric estimators involve nonparametric regressions as an intermediate step, and some of the classical examples in economics can be found in several Handbook of Econometrics chapters such as Powell (1994), Chen (2007), and Ichimura and Todd (2007). The typical assumption about $f$ in ((ref)) is that it has derivatives up to a given $(\gamma+1)$th order bounded in absolute values by a common constant everywhere or almost everywhere (a.e.). Given an estimator of $f$, an important object of interest is the convergence rate of the mean integrated squared error (MISE) of this estimator and the minimax optimality property of the MISE rate, which tells one how fast the population mean squared distance between the estimator and $f$ shrinks to zero uniformly when $f$ ranges over a smoothness class, as the sample size $n$ increases. In particular, MISE is a global mean squared error criterion by integrating over all possible values of the covariate and noise with respect to some underlying distribution (see Pagan and Ullah, 1999).
When the smoothness degree $\gamma+1$ is finite and the sample size $n\rightarrow\infty$, existing asymptotic results show that the minimax optimal MISE rate is $\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$ which decreases as $\gamma+1$ increases. In the literature of statistical learning theory, $\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$ is also stated as the nonasymptotic minimax optimal rate and derived without assuming $n\rightarrow\infty$. See Tsybakov (2009) for a comprehensive review of the literature that show the classical asymptotic minimax optimal rate $\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$ under the assumption of $\gamma+1$ being finite and $n$ tending to infinity. See Wainwright (2019) for the nonasymptotic derivation. The rate $\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$ gives arise to the so called “blessing of smoothness” folklore (i.e., the more smoothness one can exploit, the better) in the theoretical community.
In terms of the minimax optimal MISE rates associated with the standard $(\gamma+1)$th degree smoothness classes\footnote{In this paper, a standard $(\gamma+1)$th degree smoothness class refers to one that consists of functions with derivatives belonging to a ball of some constant radius independent of the derivative order and bounded away from zero and from above, with respect to either an $l_{\infty}$ (max) norm or a Hilbert norm.}, we show the following results: (i) if $n\leq\left(\gamma+1\right)^{2\gamma+3}$ (the “small $n$” regime), the minimax optimal MISE rate is $\frac{\log n}{n\log(\log n)}$ and the optimal degree of smoothness to exploit is roughly $\max\left\{ \left\lfloor \frac{\log n}{2\log\left(\log n\right)}\right\rfloor ,\,1\right\} $;\footnote{“$\log$” in this paper is natural logarithm. The more precise characterization of the optimal degree of smoothness to exploit under $n\leq\left(\gamma+1\right)^{2\gamma+3}$ is detailed in Section (ref).} (ii) if $n>\left(\gamma+1\right)^{2\gamma+3}$ (the “large $n$” regime, which clearly includes the case of $\gamma$ being finite and $n$ tending to infinity), the minimax optimal MISE rate is $\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$ and the optimal degree of smoothness to exploit is $\gamma+1$. To our knowledge, this paper is the first in the literature to show the minimax optimal rate in the “small $n$” regime and the sample size (i.e., $\left(\gamma+1\right)^{2\gamma+3}$) at which the minimax optimal rate transitions from $\frac{\log n}{n\log(\log n)}$ to $\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$.
The definition of our minimax optimal rates is based on the existing literature. A rate is said to be minimax optimal in our problem if we can show: (1) the MISE for any estimators (by taking the infimum over all estimators) in the worst case scenario (by taking the supremum over a $(\gamma+1)$th degree smoothness class) is bounded from below, and such a bound is called a minimax lower bound; (2) there exists an estimator such that its MISE in the worst case scenario has an upper bound that matches the lower bound in terms of the rate, the so-called achievability result; to be precise, apart from some universal constant independent of $n$ (and in our interest, also independent of $\gamma$), the upper bound matches the lower bound and the matching part is the minimax optimal rate. The difference in the “constants” between our results and the existing literature is that the constants in our bounds are truly universal constants that are independent of $n$ and $\gamma$ and bounded away from zero and from above, while those in the existing literature cannot be independent of $\gamma$. This difference allows us to provide an asymptotic interpretation for the rate $\frac{\log n}{n\log(\log n)}$ in the “small $n$” regime (in Section (ref)) relative to the classical rate $\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$:
Particularly, we show in this paper that, if the maximum smoothness degree of $f$ is $\gamma+1$, estimators which minimize the sum of squared residuals and are constrained to exploit the optimal degree of smoothness achieve the abovementioned minimax optimal MISE rates. These estimators will be referred to as the constrained nonparametric least squares estimator (CNLS) in the following. CNLS estimators constrained to be in a Sobolev class associated with a Reproducing Kernel Hilbert Space (RKHS) radius have nice closed form expressions via kernel functions and are easy to implement in the regularized form, often referred to as the kernel ridge regression (KRR) estimators (among the most popular nonparametric estimators) in machine learning. There is a rich theory based on RKHS for the asymptotic properties of such estimators under the regime where $\gamma+1$ is finite and $n\rightarrow\infty$ (see, e.g., Sch�lkopf and Smola, 2002; Berlinet and Thomas-Agnan, 2004). These estimators are closely related to smoothing splines methods and Gaussian process regressions (see, Wahba, 1990; Rasmussen and Williams, 2006).
It is worth mentioning the connection between our theoretical results and some numerical findings from the literature. From a practical view-point, Marron and Wand (1992) find the exact MISE of kernel density estimators based on Gaussian kernels can increase with the order of kernel being exploited (the assumed degree of smoothness) when the sample size is moderate. Marron (1994) further shows in simulations that the second order kernel produces a smaller MISE than the fourth order kernel when the sample size is between $70$ and $10000$, and the fourth order kernel is dominantly better than the second order kernel when $n>10000$. Despite that our focus in this paper is on the minimax optimal rates and our achievability results concern global nonparametric procedures such as KRR (instead of kernel density estimators based on Gaussian kernels), Marron and Wand (1992) and our paper share one general message: the classical rate $\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$ is an underestimate of the MISE when $n$ is not large enough. Another message conveyed by our paper is that the optimal degree of smoothness to exploit should depend on the sample size.
In econometrics, there is increasing empirical evidence against exploiting higher degrees of smoothness. Graham, et. al (2010) comment, “As is usual in semiparametric estimation, higher order kernels are required for bias reduction, although the use of such kernels in practice may be ill-advised.” Based on several empirical studies with sample sizes ranging from a couple of thousands to at most thirties of thousands, Gelman and Imbens (2019) recommend researchers to avoid using high order polynomials but use local linear or local quadratic polynomials to estimate the two conditional mean functions of a pretreatment variable in regression discontinuity designs (RDD) analyses.
The implication of our results extends to other applications. Establishing the convergence rate of a semiparametric procedure often requires establishing an MISE rate or its sample analogue concerning a first-step nonparametric regression. Below are a couple of examples:
Notation. Let $\left\lfloor x\right\rfloor $ denote the largest integer smaller than or equal to $x$. For two functions $f(n,\gamma)$ and $g(n,\gamma)$, let us write $f(n,\gamma)\succsim g(n,\gamma)$ if $f(n,\gamma)\geq cg(n,\gamma)$ for a universal constant $c\in(0,\,\infty)$; similarly, we write $f(n,\gamma)\precsim g(n,\gamma)$ if $f(n,\gamma)\leq c^{'}g(n,\gamma)$ for a universal constant $c^{'}\in(0,\,\infty)$; and $f(n,\gamma)\asymp g(n,\gamma)$ if $f(n,\gamma)\succsim g(n,\gamma)$ and $f(n,\gamma)\precsim g(n,\gamma)$. Throughout this paper, we use various $c$ and $C$ letters to denote positive universal constants that are: finite and bounded away from zero (denoted by $\asymp1$) and independent of $n$ and $\gamma$ and the dimension $d$ of the covariates (when $d-$dimensional covariates are of interest); these constants may vary from place to place.
For a $J-$dimensional vector $\theta$, the $l_{q}-$norm $\left|\theta\right|{}_{q}:=\left(\sum_{j=1}^{J}|\theta_{j}|^{q}\right)^{1/q}$ if $1\leq q<\infty$ and $\left|\theta\right|{}_{q}:=\max_{j\in\left\{ 1,...,J\right\} }\left|\theta_{j}\right|$ if $q=\infty$. Let $\mathbb{B}_{q}^{J}\left(R\right):=\left\{ \theta\in\mathbb{R}^{J}:\left|\theta\right|_{q}\leq R\right\} $. For functions on $\left[a,\,b\right]$, the unweighted $L^{2}-$norm $\left|f-g\right|_{2}:=\sqrt{\int_{a}^{b}\left[f\left(x\right)-g\left(x\right)\right]^{2}dx}$, and the weighted $L^{2}\left(\mathbb{P}\right)-$norm $\left|f-g\right|_{2,\mathbb{P}}:=\sqrt{\int_{a}^{b}\left[f\left(x\right)-g\left(x\right)\right]^{2}\mathbb{P}(dx)}$.
For functions on $\left[a,\,b\right]^{d}$, the supremum norm $\left|f-g\right|_{\infty}:=\sup_{x\in\left[a,\,b\right]^{d}}\left|f\left(x\right)-g\left(x\right)\right|$.
Finally, the $\mathcal{L}^{2}(\mathbb{P}_{n})-$norm of the vector $f:=\left\{ f(x_{i})\right\} _{i=1}^{n}$, denoted by $\left|f\right|_{n}$, is $\left[\frac{1}{n}\sum_{i=1}^{n}\left(f(x_{i})\right)^{2}\right]^{\frac{1}{2}}$.
Any function $f$ in a H�lder class on $\left[-1,\,1\right]$ can be written as
where $z$ is some intermediate value between $x$ and $0$, and we follow the convention $0^{0}=0!=1$.
Consequently, we have the following relationships:
As discussed in Wainwright (2019, Chapter 12), any function $f$ in the Sobolev space on $\left[0,\,1\right]$ has the expansion \[ f\left(x\right)=\sum_{k=0}^{\gamma}f^{(k)}(0)\frac{x^{k}}{k!}+\int_{0}^{1}f^{(\gamma+1)}(t)\frac{\left(x-t\right)_{+}^{\gamma}}{\gamma!}dt\quad\textrm{with }(a)_{+}=a\vee0, \] and one (RK)HS norm associated with a Sobolev space takes the form
Therefore, the Sobolev space can be decomposed into a polynomial subspace and a Sobolev subspace imposed with the restrictions that $f^{(k)}(0)=0$ for all $k=0,...,\gamma$ and $f^{(\gamma+1)}$ belongs to the space $\mathcal{L}^{2}\left[0,\,1\right]$ (see Wahba, 1990, Chapter 1; Wainwright, 2019, Chapter 12, Examples 12.17 and 12.29). A Sobolev subclass is a special case of the ellipsoid subclass; in particular, a Sobolev subclass can be expressed in the following way (Wainwright 2019, Chapters 5, 12 and 15).
When considering $\mathcal{H}_{\gamma+1}$ in ((ref)), we will assume that $\mu_{m}=\left(cm\right)^{-2\left(\gamma+1\right)}$ for a positive constant $c\asymp1$ independent of $\gamma$ and $R_{\gamma+1}$. The decay rate of the eigenvalues follows the standard assumption for $(\gamma+1)-$degree smooth functions in the literature (see, e.g., Steinwart and Christmann, 2008; Wainwright, 2019) and $R_{\gamma+1}\asymp1$ in ((ref)) gives the standard ellipsoid subclass. Moreover, ((ref)) is equipped with the inner product $\left\langle h,\,g\right\rangle _{\mathcal{H}}=\sum_{m=1}^{\infty}\frac{\left\langle h,\,\phi_{m}\right\rangle \left\langle g,\,\phi_{m}\right\rangle }{\mu_{m}}$ where $\left\langle \cdot,\cdot\right\rangle $ is the inner product in $\mathcal{L}^{2}\left[0,\,1\right]$.
The following is a standard textbook result that summarizes the relationships between covering and packing numbers:
Given this sandwich result, a lower bound on the packing number gives a lower bound on the covering number, and vice versa; similarly, an upper bound on the covering number gives an upper bound on the packing number, and vice versa.
Metric entropy is an important concept in approximation theory and discrete geometry. In mathematical statistics and machine learning theory, metric entropy is a fundamental building block. Combined with the Fano's inequality from information theory (see, Cover and Thomas, 2005), it allows one to derive the minimax lower bounds for the MISE rates; combined with the notion of “local complexity” in empirical processes theory, it allows one to derive upper bounds on the MISE rates.\footnote{For the standard smoothness classes, one may use methods based on “ranks” (for the polynomial subclass) and “eigenvalues” (for the nonparametric subclass) to derive upper bounds on the MISE. However, “ranks” are not very useful for deriving the minimax lower bounds in general. Even for upper bounds in the case of a H�lder class, once we allow $R_{k}$ to depend on the order of derivative, $k$, the “rank”--based argument is hard to generalize as it does not account for the impact of $R_{k}$.}
In the existing literature, the minimax lower bound and the upper bound (for achievability) for MISE is typically stated as $\mathtt{constant}\cdot\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$ and $\mathtt{constant'}\cdot\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$, respectively, where the specific forms of $\mathtt{constant}$ and $\mathtt{constant}'$ are not derived. Then the matching part “$\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$” is claimed as the minimax optimal rate (e.g., Wainwright, 2019). Deriving constants for global criteria such as MISE and demonstrating achievability (upper bounds) in the context of global nonparametric procedures (the focus of our paper) is known to be difficult, which is why the existing literature does not make an attempt in deriving the constants.
Relative to the existing literature, our results take one step further by revealing more explicit dependence on $n$, $\gamma$ and $\left\{ R_{k}\right\} _{k=0}^{\gamma+1}$. The difference in the “constants” between our results and the existing literature is that the constants in our bounds are truly universal constants that are independent of $n$ and $\gamma$ and bounded away from zero and from above, while those in the existing literature cannot be independent of $\gamma$, as explained in the next section. This difference allows us to provide an asymptotic interpretation for the rate $\frac{\log n}{n\log(\log n)}$ in the “small $n$” regime (in Section (ref)) relative to the classical rate $\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$, as discussed in Section (ref).\footnote{Because of the complexity of our problems, we make no attempt to derive the explicit universal constants that are independent of $n$ and $\gamma$. The following process is a standard one used in this literature for proving minimax optimal rates in MISE and explains why such a derivation is not carried out: First, results on metric entropy are derived. Then, these results are applied with other technical lemmas to establish minimax optimal rates for the MISE. Throughout, the universal constants become complex rather quickly and it becomes physically impossible to track them from lines to lines. Based on the existing mathematical tools from approximation theory, information theory, minimax lower bounds theory, and empirical processes theory (for upper bounds), deriving meaningful universal constants (letting alone sharp universal constants) in the particular context of this paper is infeasible and remains an open problem until a new mathematical foundation can be introduced.}
The fundamental contribution of this paper lies in a set of metric entropy bounds. Some of them are original, and some of them improve and/or extend the ones in the literature to allow general (possibly $k-$dependent) $R_{k}$s. Besides their applications in the empirical processes theory and the theory of minimax lower bounds, our bounds for the covering and packing numbers are of independent interest for researchers working in approximation theory and discrete geometry.
Here we start with a short discussion of why the phase transition phenomenon has been overlooked in the literature while leaving the details to the subsequent sections. For $\mathcal{U}_{\gamma+1}$ under the assumption that $R_{k}=\overline{C}\asymp1$ (which we will refer to as the standard $\mathcal{U}_{\gamma+1}$), Kolmogorov and Tikhomirov (1959) show that
for $\delta-$approximation accuracy. By ((ref)), ((ref)) implies
The minimax optimal MISE rates (as functions of $n$ and $\gamma$) are related to choices of $\delta$ (see, e.g., Yang and Barron, 1999; Wainwright, 2019). Virtually every paper including recent textbooks on nonasymptotic statistics (such as Wainwright, 2019, eq. 5.17 on p. 129) takes $\log(\delta-\textrm{covering/packing number})\asymp\delta^{\frac{-1}{\gamma+1}}$. It would take some diligence for one to recognize that:
Here is a list of our discoveries:
All these issues are addressed in this paper. Particularly, in deriving the metric entropy bounds for the polynomial subclass $\mathcal{U}_{\gamma+1,1}$, we develop our own arguments; in deriving the metric entropy bounds for the ellipsoid subclass, we base our arguments on the existing literature but use an improved truncation strategy that gives our resulting bounds the sharp dependence on $\gamma$ and $R_{\gamma+1}$. In contrast to those in the existing literature, our entropy bounds enable us to reveal the phase transition phenomenon and obtain the sharp MISE rates in the respective “small $n$” and “large $n$” regimes, associated with the standard smoothness classes as well as non-standard smoothness classes motivated in Zhu and Mirzaei (2021).
The results and proofs of the paper are organized as follows.
In this section, we revisit the minimax optimal MISE rates associated with some commonly seen smoothness classes in the literature.
When we write $\mathcal{U}_{\gamma+1}$ or $\mathcal{S}_{\gamma+1}$ in this section, $\mathcal{U}_{\gamma+1}$ refers to the standard H�lder class and $\mathcal{S}_{\gamma+1}$ refers to the standard Sobolev class. The regression model ((ref)) is subject to the following conditions.
The proofs for the minimax optimal rates in this section as well as Appendix (ref) are based on our results in Section (ref) as well as techniques from empirical processes, machine learning theory, and information theory, collected in Appendix (ref).
Assumption (ref) is the most standard in the literature for minimax lower bounds; see, e.g., Yang and Barron (1999), Raskutti, et. al (2011), and Wainwright (2019). The main reason is, a minimax lower bound derived under a more restrictive model would be a lower bound under a less restrictive model that includes the former as a special case. For example, a homoscedastic model is a special case of a heteroscedastic one, and the independence of $\varepsilon_{i}$ and $X_{i}$ along with $\mathbb{E}\left(\varepsilon_{i}\right)=0$ implies $\mathbb{E}\left(\varepsilon_{i}\vert X_{i}\right)=0$. Meanwhile, upper bounds for known estimators under less restrictive model assumptions are typically available or can be established before developing the minimax lower bounds. See the introduction of Chapter 15 in Wainwright (2019) for the insights provided by minimax lower bounds.
Until very recently, Zhao and Yang (2022) find a way to allow the noise terms to be heteroscedastic and have non-Gaussian distributions. Specifically, Zhao and Yang (2022) consider
where, conditional on $X_{i}$, $Y_{i}$ has a density from a location-scale family subjective to some rather mild conditions. This setup allows Zhao and Yang (2022) to derive an upper bound for the Kullback-Leibler divergence between distributions that are not Gaussian (but can be Gaussian). The bound in Zhao and Yang (2022) can be used to relax homoscedasticity without changing the main arguments of our proofs. For more details related to proofs under ((ref)), we refer the interested readers to Appendix (ref) and the remark following Lemma C.1 in Appendix (ref) of our paper. It is worth pointing out that, when $\sigma(\cdot)$ in ((ref)) ranges over the standard smoothness classes with positive members bounded away from zero and bounded from above, the minimax lower bounds have the same rate as those under Assumption (ref).
As for the upper bounds (that demonstrate achievability), our proofs are developed under a weaker assumption: conditional on $\left\{ X_{i}\right\} _{i=1}^{n}$, $\left\{ \varepsilon_{i}\right\} _{i=1}^{n}$ are independent sub-Gaussian random variables with zero mean and sub-Gaussian parameters at most $\overline{\sigma}\asymp1$. This assumption is very standard in the literature of empirical process theory (e.g., van der Vaart and Wellner, 1996; van de Geer, 2000; Zhu, 2017; Wainwright, 2019). However, to establish minimax optimality, it is more sensible to consider a model with the same set of assumptions on $\left\{ \varepsilon_{i}\right\} _{i=1}^{n}$ for the minimax lower bounds and upper bounds. Therefore, to focus on the key points, we will state the results under Assumption (ref) in the main paper and Appendix (ref) to be consistent with the way the well-known results are presented in the literature for minimax optimality.
The proof for Theorem (ref) is given in Section (ref), which relies on the constructions behind the bound $\underline{B}_{2}$ in Lemma (ref), Lemma (ref) and Lemma (ref) in Section (ref).
Recall from Section (ref) that a smoothness class contains two subclasses: a polynomial subclass and an infinite-dimensional nonparametric subclass such that any function $f$ in this subclass has $f^{(k)}\left(0\right)=0$ for all $k=0,...,\gamma$. The proof for Theorem (ref) applies a version of the Fano's inequality, which converts the problem into a multiple classification problem that tries to distinguish among $M$ members in the function classes of interest. The set of $M$ members need be sufficiently large and is naturally connected with a packing set of the function class.
Based on the constructions behind the bound $\underline{B}_{2}$ in Lemma (ref), we can show, the error rate that prevents from distinguishing any pair of members in the packing set of the standard $\gamma$th degree polynomial subclass (in $\mathcal{U}_{\gamma+1}$ and $\mathcal{S}_{\gamma+1}$, respectively) is $\frac{\sigma^{2}\left(\gamma+1\right)}{n}$. Based on Lemmas (ref) and (ref), we can show, the error rate that prevents from distinguishing any pair of members in the packing set of the standard $(\gamma+1)$th degree H�lder subclass and Sobolev subclass is $\left(\frac{\sigma^{2}}{n}\right)^{\frac{2\left(\gamma+1\right)}{2\left(\gamma+1\right)+1}}$. Note that
In deriving the minimax lower bounds, the existing nonasymptotic literature uses the metric entropy $\delta^{\frac{-1}{\gamma+1}}$ (e.g., Wainwright, 2019, Examples 5.11 and 15.23). The entropy $\delta^{\frac{-1}{\gamma+1}}$ ignores the polynomial subclass. For more details, we refer the interested readers to Section (ref). Using $\delta^{\frac{-1}{\gamma+1}}$ to derive the minimax lower bound results in a rate $\left(\frac{\sigma^{2}}{n}\right)^{\frac{2\left(\gamma+1\right)}{2\left(\gamma+1\right)+1}}$ and overlooks the (sharper) rate in Theorem (ref)(ii), the “small $n$” regime.
The following two theorems show the achievability of the rates in Theorem (ref). In particular, we show that estimators constrained to exploit the optimal degree of smoothness in each regime achieve the respective rate in Theorem (ref). We consider the constrained nonparametric least squares estimator (CNLS)
The CNLS is commonly seen in the minimax optimality literature (e.g., Raskutti, et. al 2011) and the Kernel Ridge Regression (KRR) literature. Let $k=\gamma$ in the “large $n$” regime and $k=\gamma^{*}$ in the “small $n$” regime. We consider either $\mathcal{F}=\mathcal{U}_{k+1}$ or
where $\left|\cdot\right|_{\mathcal{H},k+1}$ is the norm defined in ((ref)). Both cases can be of interest but the latter is more widely implemented in practice, as we explain below.
Constraining the estimators to be in $\mathcal{S}_{k+1}$ allows one to implement ((ref)) via kernel functions. The following discussion is based on Wainwright (2019). Let the matrix $\mathbb{K}_{k+1}\in\mathbb{R}^{n\times n}$ consist of entries $\frac{1}{n}\mathcal{K}_{k+1}\left(x_{i},\,x_{i^{'}}\right)$, taking the form
where $\left(a\right)_{+}=a\vee0$. The kernel function $\mathcal{K}_{k+1,1}\left(w,\,z\right)=\sum_{j=0}^{k}\frac{w^{j}}{j!}\frac{z^{j}}{j!}$ generates the $k$th degree polynomial subspace and the kernel function $\mathcal{K}_{k+1,2}\left(w,\,z\right)=\int_{0}^{1}\frac{\left(w-t\right)_{+}^{k}}{k!}\frac{\left(z-t\right)_{+}^{k}}{k!}dt$ for $k>0$ ($w\wedge z$ for $k=0$) generates the $(k+1)$th order Sobolev subspace imposed with the restrictions that $f^{(j)}(0)=0$ for all $j=0,...,k$ and $f^{(k+1)}$ belongs to the space $\mathcal{L}^{2}\left[0,\,1\right]$.
When $\mathcal{F}=\mathcal{S}_{k+1}$ in ((ref)), $\hat{f}$ can be written as
where
In particular, ((ref)) comes from the representation $\left|\check{f}\right|_{\mathcal{H},k+1}^{2}=\pi^{T}\mathbb{K}_{k+1}\pi$ when $\mathcal{F}=\mathcal{S}_{k+1}$ in ((ref)) and $\check{f}$ takes the form $\check{f}\left(\cdot\right)=\frac{1}{\sqrt{n}}\sum_{i^{'}=1}^{n}\pi_{i^{'}}\mathcal{K}_{k+1}\left(\cdot,\,x_{i^{'}}\right)$.
The program above is convex. Hence, by the Lagrangian duality, solving ((ref))--((ref)) is equivalent to solving
for a properly chosen regularization parameter $\lambda>0$ (detailed in the theorems concerning KRR). Consequently, the optimal weight vector $\hat{\pi}$ takes the form
where $Y=\left\{ Y_{i}\right\} _{i=1}^{n}$. The KRR estimators associated with $\mathcal{S}_{k+1}$ are related to the smoothing splines methods and Gaussian process regressions in machine learning. We refer interested readers to Wahba (1990), Sch�lkopf and Smola (2002), as well as Rasmussen and Williams (2006) for more details.
The proof for Theorem (ref) is given in Appendix (ref).
If the smoothness degree $\gamma+1$ is known, there is a conjecture in the literature that cross-validation will yield a regularization parameter with the optimal order $\left(\frac{1}{n}\right)^{\frac{2\left(\gamma+1\right)}{2\left(\gamma+1\right)+1}}$ that $\lambda$ has in Theorem (ref) under the classical setting where the smoothness degree $\gamma+1$ is finite and $n\rightarrow\infty$ (e.g., van de Geer, 2000).\footnote{van de Geer (2000) states “the optimal order” as $\left(\frac{1}{n}\right)^{\frac{\gamma+1}{2\left(\gamma+1\right)+1}}$ as her $\lambda$ corresponds to our $\sqrt{\lambda}$.} To the best of our knowledge, this conjecture has not been proved yet. Our results prompt a much harder question of whether cross validation can also be used to choose a smoothness degree that adapts to the phase transition; that is, in the regime $\frac{n}{\sigma^{2}}\leq\left(\gamma+1\right)^{2\gamma+3}$, choosing a smoothness degree with the order $\gamma^{*}+1$ has in Theorem (ref), and in the regime $\frac{n}{\sigma^{2}}>\left(\gamma+1\right)^{2\gamma+3}$, choosing $\gamma+1$ (also see Section (ref) for a related discussion about rate adaptive estimators).
The takeaway for applied researchers. The conjecture and question above are mostly of theoretical interest. Even if they are proved eventually, using cross-validation to choose a smoothness degree and the corresponding regularization parameter is computationally costly. For applied researchers, the simplest takeaway from our theoretical results is similar to the one given in the applied literature -- exploit no more than two degrees of smoothness (see Section (ref) for more discussions).
Like the Sobolev class, we can show a similar achievability result for the H�lder class.
The proof for Theorem (ref) is given in Appendix (ref).
Table 1 summarizes the results in Section (ref) for easy reference.
When deriving the upper bounds in Theorems (ref) and (ref), we obtain the following bounds on the sample mean squared error (SMSE) as intermediate results.
In this subsection, we consider infinitely differentiable functions with the following form
The standard $\mathcal{U}_{\infty}$ for infinitely smooth functions follows the definition of the standard $\mathcal{U}_{\gamma+1}$ with $\gamma+1=\infty$. The analog of $\mathcal{S}_{\gamma+1}$ in ((ref)) for infinitely smooth functions can be defined as
where $\overline{C}\asymp1$ is a universal constant independent of $n$ and $\gamma$.
We have the following result for analytic functions.
The proof for the lower bound in Corollary (ref) is identical to the proof for part (ii) of Theorem (ref) (Appendix (ref)). The proofs for the upper bounds in Corollary (ref) involve only slight modifications of Appendix (ref) and Appendix (ref). See Appendix (ref).
There is no phase transition in the optimal rate associated with $\mathcal{U}_{\infty}$ and $\mathcal{S}_{\infty}$. As discussed earlier, $\frac{\gamma^{*}+1}{n}\asymp\frac{\log n}{n\log(\log n)}$. Relatedly, see a review article about the estimation of analytic functions by Ibragimov (2001), in particular, Theorem 4.9 where the minimax optimal rate for the MISE associated with analytic functions is $\frac{\log n}{n\log(\log n)}$.\footnote{The author thanks an anonymous referee for pointing out this reference.}
Corollary (ref) highlights the importance of our results on phase transition in this paper. Stating $\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$ as the minimax optimal rate for MISE associated with a standard smoothness class without a condition on $n$ raises the following “paradox”. Recall ((ref)); on the other hand, $\mathcal{U}_{\infty}\subseteq\mathcal{U}_{\gamma+1}$ so ((ref)) is not plausible if $\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$ is stated alone, without a condition on $n$, as the minimax optimal rate in MISE for the standard $\mathcal{U}_{\gamma+1}$. Our results in Section (ref) are crucial for reconciling the aforementioned “paradox”, as the choice of our critical smoothness parameter $\gamma^{*}+1$ guarantees that:
Assuming $\sigma\asymp1$, a better practice would be to state the rate $\left(\frac{\sigma^{2}}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$ under the condition $\frac{n}{\sigma^{2}}>\left(\gamma+1\right)^{2\gamma+3}$ and the rate $\frac{\sigma^{2}\log n}{n\log(\log n)}$ under the condition $\frac{n}{\sigma^{2}}\leq\left(\gamma+1\right)^{2\gamma+3}$, as shown in this paper.
In this section, we present bounds on covering and packing numbers associated with $\mathcal{U}_{\gamma+1,1}$, $\mathcal{U}_{\gamma+1,2}$, and $\mathcal{H}_{\gamma+1}$. Various $c$ and $C$ letters in this section denote positive universal constants that are: finite and bounded away from zero (denoted by $\asymp1$) and independent of $\gamma$ and $\left\{ R_{k}\right\} _{k=0}^{\gamma+1}$ (parameters of the function classes); these constants may vary from place to place.
Table 2 summarizes the results in this section for easy reference.
The proof for Lemma (ref) is given in Appendix (ref).
The lower bounds $\underline{B}_{1}\left(\delta\right)$ and $\underline{B}_{2}$, as well as the upper bound $\overline{B}_{1}\left(\delta\right)$ are original. The (less original) bound $\overline{B}_{2}\left(\delta\right)$ generalizes the upper bound associated with the polynomial subclass in Kolmogorov and Tikhomirov (1959), which takes the form $\left(\gamma+1\right)\log\frac{1}{\delta}$. It is worth pointing out that $\overline{B}_{2}\left(\delta\right)$ holds for all $\delta\in(0,\,1)$ (not just $\delta$ such that $\min_{k\in\left\{ 0,...,\gamma\right\} }\log\frac{4\left(\gamma+1\right)R_{k}}{k!\delta}<0$) but is far from being tight when $\min_{k\in\left\{ 0,...,\gamma\right\} }\log\frac{4\left(\gamma+1\right)R_{k}}{k!\delta}\geq0$. Obviously $\overline{B}_{1}\left(\delta\right)\precsim\overline{B}_{2}\left(\delta\right)$. When it comes to deriving the minimax optimal rates for the MISE under large enough $R_{k}$ (e.g., Theorems (ref) and (ref) in Appendix (ref)), $\overline{B}_{1}\left(\delta\right)$ will be very useful.
When deriving a lower bound for the packing number of the standard $\mathcal{U}_{\gamma+1}$ under the assumption that $R_{k}=\overline{C}\asymp1$, Kolmogorov and Tikhomirov (1959) constructs a set of functions where $\mathcal{U}_{\gamma+1,1}$ is a singleton; therefore, the cardinality of this set only gives a lower bound for the packing number of $\mathcal{U}_{\gamma+1,2}$ and overlooks $\mathcal{U}_{\gamma+1,1}$. This issue is further discussed in Section (ref). Our Lemma (ref) establishes two different lower bounds for the packing number of $\mathcal{U}_{\gamma+1,1}$. In particular, the construction of $\underline{B}_{2}$ in ((ref)) will be useful for deriving the minimax lower bounds for the MISE when $R_{k}$ is relatively small (e.g., Theorem (ref) in Section (ref) and Theorem (ref) in Appendix (ref)), while $\underline{B}_{1}\left(\delta\right)$ in ((ref)) will be useful when $R_{k}$ is relatively large (Theorem (ref) in Appendix (ref)).
To establish $\overline{B}_{1}\left(\delta\right)$, $\underline{B}_{1}\left(\delta\right)$ and $\underline{B}_{2}$, we discard the argument in Kolmogorov and Tikhomirov (1959) and develop our own. The derivation of $\underline{B}_{2}$ is based on a constructive proof. To derive $\overline{B}_{1}\left(\delta\right)$ and $\underline{B}_{1}\left(\delta\right)$, we consider two classes (equivalent to $\mathcal{U}_{\gamma+1,1}$), each involving a $\left(\gamma+1\right)-$dimensional polyhedron. The lower bound $\underline{B}_{1}\left(\delta\right)$ is the more delicate part. In particular, for any $f\in\mathcal{U}_{\gamma+1,1}$, we write $f\left(x\right)=\sum_{k=0}^{\gamma}\tilde{\theta}_{k}\phi_{k}\left(x\right)$, where $\left(\phi_{k}\right)_{k=0}^{\gamma}$ are the Legendre polynomials.
The proof for Lemma (ref) is given in Appendix (ref). Lemma (ref) and Lemma (ref) together are used to prove Theorems (ref) and (ref) in Section (ref), as well as Theorems (ref) and (ref) in Appendix (ref); in addition, Lemma (ref) is also used to prove Theorem (ref)(ii) in Appendix (ref).
Lemma (ref) extends Kolmogorov and Tikhomirov (1959) to allow for general $R_{k}$s. When $R_{k}\leq\overline{C}k!$ for all $k=1,...,\gamma+1$, $R^{*\frac{1}{\gamma+1}}\asymp1$. If $R_{k}\geq\overline{C}k!$ for all $k=1,...,\gamma+1$, $R^{*\frac{1}{\gamma+1}}\succsim1$; for example, taking $R_{k}\geq\overline{C}\left(k!\right)^{2}$ for all $k=1,...,\gamma+1$ yields $R^{*\frac{1}{\gamma+1}}\succsim\gamma$.
Given Lemmas (ref) and (ref), ((ref)) and ((ref)), we have
and \[ \log M_{\infty}\left(\delta,\,\mathcal{U}_{\gamma+1}\right)\succsim\log M_{2}\left(\delta,\,\mathcal{U}_{\gamma+1}\right)\succsim
. \]
Our lower bounds above sharpen the classical result in Kolmogorov and Tikhomirov (1959). In particular, the lower bound for $\mathcal{U}_{\gamma+1}$ in Kolmogorov and Tikhomirov (1959) (derived under the assumption that $R_{k}=\overline{C}\asymp1$) takes the form $\delta^{\frac{-1}{\gamma+1}}$. This result and its proof are inherited later in papers and textbooks including the more recent textbook on nonasymptotic statistics by Wainwright (2019, Example 5.11), where in the derivation of the lower bound, a set of functions are constructed in the way such that their $k$th order derivatives evaluated at zero are zero for all $k=0,...,\gamma$. In other words, $\mathcal{U}_{\gamma+1,1}$ is a singleton in this construction and the cardinality of this set only gives a lower bound for $\mathcal{U}_{\gamma+1,2}$. In particular, the lower bound $\delta^{\frac{-1}{\gamma+1}}$ is not sharp when $\gamma+1$ and $\delta$ are large enough.
The proof for Lemma (ref) is given in Appendix (ref). Lemma (ref) and Lemma (ref) together are used to prove Theorem (ref). Lemma (ref) is also used to prove Theorem (ref)(i) in Appendix (ref).
When $R_{\gamma+1}=1$, Lemma (ref) sharpens the upper bound for $\log N_{2}\left(\delta,\,\mathcal{H}_{\gamma+1}\right)$ in Wainwright (2019) from $\left(\gamma\vee1\right)\delta^{-\frac{1}{\gamma+1}}$ to $\delta^{-\frac{1}{\gamma+1}}$; in particular, the upper and lower bounds in Wainwright (2019) (the last two inequalities on p.131) scale as $\left(\gamma\vee1\right)\delta^{\frac{-1}{\gamma+1}}$ and $\delta^{\frac{-1}{\gamma+1}}$, respectively, while our upper and lower bounds in Lemma (ref) have the same scaling $\delta^{\frac{-1}{\gamma+1}}$. We discover the cause of the gap lies in that the “pivotal” eigenvalue (that balances the “estimation error” and the “approximation error” from truncating for a given resolution $\delta$) in Wainwright (2019) is not optimal. The truncation in Wainwright (2019) is commonly used in the existing literature and seems to originate from Theorem 3 in Mityagin (1961). We close the gap by finding the optimal “pivotal” eigenvalue.
More generally, for the case of $R_{\gamma+1}\precsim\gamma+1$, we consider two different truncations, one giving the upper bound $\delta^{\frac{-1}{\gamma+1}}$ and the other giving the lower bound $\left(R_{\gamma+1}\delta^{-1}\right)^{\frac{1}{\gamma+1}}$. Note that $\left(R_{\gamma+1}\delta^{-1}\right)^{\frac{1}{\gamma+1}}\asymp\delta^{\frac{-1}{\gamma+1}}$ when $R_{\gamma+1}\asymp1$. For the case of $R_{\gamma+1}\succsim\gamma+1$, we use only one truncation to show that both the upper bound and the lower bound scale as $\left(R_{\gamma+1}\delta^{-1}\right)^{\frac{1}{\gamma+1}}$.
In the discussion following Theorem (ref), we have brought up the trade-off between the polynomial subclass and the nonparametric subclass. It is worth mentioning the connection between our theoretical results and the empirical findings from Gelman and Imbens (2019), which implements high order polynomials in regression discontinuity designs (RDD) analyses. In particular, when applying RDD to perform causal inference, two conditional mean functions of a pretreatment variable are estimated from ((ref)). There are several empirical issues of using high order polynomials raised in Gelman and Imbens (2019). The implication of our results is most related to their paper on the issue of mean squared errors (MSE).
Regarding the data sets studied in Gelman and Imbens (2019), the Jacob-Lefgren data (Jacob and Lefgren, 2004), Lee data (Lee, 2008), Matsudaira data (Matsudaira, 2008), the LaLonde data (LaLonde, 1986), and the census data in 1974, 1975 and 1978, the sample sizes used in the implementation of Gelman and Imbens (2019) range from thousands to at most thirties of thousands. Based on their empirical evidence from studying these data sets, Gelman and Imbens (2019) recommend researchers to avoid using high order polynomials but use local linear or local quadratic polynomials. In view of the asymptotic condition $\left(\gamma+1\right)^{2\gamma+3}=o(n)$ for the rate $\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$ to kick in (see ((ref))), let us solve for $\overline{\gamma}+1$ roughly from $n\geq\left(\overline{\gamma}+1\right)^{2\overline{\gamma}+3}$. This heuristic gives a rough degree of smoothness to exploit, $\max\left\{ \left\lfloor \frac{\log n}{2\log\left(\log n\right)}\right\rfloor ,\,1\right\} $, which is quite close to the recommended smoothness degrees in Gelman and Imbens (2019).
We conclude the paper by discussing a few open questions motivated by this work. First, our focus in this paper is on the global criterion MISE while the concern of Gelman and Imbens (2019) is about the MSE of the high order polynomial implementation at a point. The minimax optimality of pointwise MSE and the global MISE (concerning an entire function) would involve different proofs. It is well known that the minimax optimal MISE rate coincides with the minimax optimal pointwise MSE rate in the regime where $\gamma$ is finite and $n\rightarrow\infty$ (see Tsybakov, 2009). We conjecture that the minimax optimal rates would be the same for the MISE and pointwise MSE in our “small $n$” regime, simply because the trade-off between the polynomial subclass and the nonparametric subclass exists whether the interest is the MISE or pointwise MSE.
Second, discussions of related literature in Section (ref) indicate that the classical rate $\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$ is an underestimate of the MISE for local smoothing methods such as kernel density estimators and local polynomials when $n$ is not large enough. For this problem, we could consider ((ref)) where $X_{i}=\frac{i}{n}$ for $i=1,...,n$ and $\left\{ \varepsilon_{i}\right\} _{i=1}^{n}$ satisfies the assumptions in Corollary 2.3 of Tsybakov (2009). This setup is simpler than the one considered in this paper, but serves a good starting point. Like how we establish the results in this paper, we would first show the minimax lower bound under the “small $n$” regime, and then show that the MISE of a local smoothing method has an upper bound that matches the lower bound up to some universal constant independent of $n$ and $\gamma$. The proofs would be different from the ones in this paper. There is some theoretical evidence (although not a proof) suggesting that it would require a large $n$ for higher order local polynomials to become beneficial; for example, Tsybakov (2009) requires the smallest eigenvalue associated with the local polynomials to be bounded away from zero (Assumption LP1) to establish the upper bound $\left(\frac{1}{n}\right)^{\frac{2\gamma+2}{2\gamma+3}}$. This eigenvalue condition in Tsybakov (2009) requires a large enough $n$ and a sufficient condition given in Tsybakov (2009) is that $n\rightarrow\infty$.
Third, our results prompt a hard question of whether a practical estimator can be developed to adapt to the phase transition shown in this paper. Cross validation is one possibility as discussed in Section (ref). In the literature, an alternative construction of adaptive estimators uses the Lepski's method; see, for example, Chen et. al (2021) where the unknown degree of smoothness is assumed to be fixed while $n\rightarrow\infty$. Hence, the construction of smoothness degrees in Chen et. al (2021) does not depend on $n$. Chen raises an interesting point: a better construction of smoothness degrees should depend on $n$, based on our results in this paper (personal communication, October and November 2022). However, it is not clear if this approach will yield an estimator that automatically switches the degree of smoothness to adapt to the phase transition in the rates, which, after all, is a nonlinear phenomenon.