Approximate Sparsity Class and Minimax Estimation

Approximate Sparsity Class and Minimax Estimation

\allowdisplaybreaks

Keywords: Metric entropy; Minimax estimation; Nonparametric estimation; LASSO

JEL classification codes: C1; C13; C14

Introduction

First introduced by cencov62,kromaltarter68,schwartz67,watson69, density estimation using orthogonal series has since been extensively studied. Suppose we observe an i.i.d sample $\{X_i\}_{i=1}^n$ from the distribution of a random variable $X$ on $[0,1]$ with probability density $f_X$. Let $\{\phi_j\}_{j=1}^\infty$ be an orthonormal basis of $L^2([0,1],\mu)$. If $f_X\in L^2([0,1],\mu)$, then it enjoys an expansion $f_X(\cdot) = \sum_{j=1}^\infty\theta_j\phi_j(\cdot)$, where $\theta_j = E\left[\phi_j(X)\right]$ is the j-th Fourier coefficients for all $j\geq 1$. A natural estimator then takes the form

When evaluating the performance of such estimator using MISE criteria, the cutoff value $J$ plays the role of a tuning parameter, which, when chosen properly, balances the variance and bias and hence minimizes MISE. The optimal choices of the cutoff $J$ have been discussed extensively in the literature, see, for example, hall87,hart85,kromaltarter76,watson69. A generalized version of such estimator is based on thresholding

where $\omega_j$'s are the so-called thresholding parameters and typically $\omega_j \in [0,1]$ plays a role of shrinking the estimated coefficients $\hat\theta_j$'s. Note that the previous estimator $\hat{f}_J(x)$ is a special case of the thresholding estimator ((ref)) when we use $\omega_j = \mathbf{1}\{j\leq J\}$. The thresholding estimator of the form ((ref)) was first considered by kromaltarter68, and many thresholding procedures have since been extensively studied, see for example, bunea2010spades,chickencai05,digglehall86,donoho96,efro86,Wahba81 and the references therein.

Various previous works have shown that, when the thresholding parameters $\omega_j$'s are chosen properly, the orthogonal series estimators of the form ((ref)) can achieve minimax rates of convergence over familiar function classes such as Sobolev ellipsoids and Besov spaces (commonly discussed in the context of wavelets thresholding), see for example, bunea2010spades,chickencai05,donoho96,donoho98,efro86,hall86,hardle2012wavelets and the references therein. We note that many function classes considered, such as the Sobolev ellipsoids, are characterized by the restrictions on the Fourier coefficients of the functions in those classes. Those restrictions on the Fourier coefficients in turn help researchers establish the statistical properties of estimators of the form ((ref)).

The type of restriction of Fourier coefficients that we are interested in concerns how fast Fourier coefficients decay. It is motivated, for instance, by the discussion in efro2008,hall86 that for a twice differentiable density $f$ on $[0,1]$ with bounded second derivative, its Fourier expansion $f(\cdot) = \sum_{j=1}^\infty \theta_j\phi_j(\cdot)$ given the cosine orthonormal basis $\{\phi_j\}_{j=1}^\infty$ has the property that the Fourier coefficients $\{\theta_j\}$ decay at rate $j^{-2}$. Similar results concerning Hermite orthonormal basis were discussed in schwartz67. In this paper, we generalize such restrictions. In particular, we allow for that: (i) the non-increasing reordering (in absolute value) of the Fourier coefficients satisfies that the j-th largest Fourier coefficient decays at rate $|\theta_{(j)}|\leq Aj^{-k}$ for some constants $A,k$; (ii) the tail sum of the Fourier coefficients satisfies $\sum_{j=J+1}^\infty \theta_j^2\leq CJ^{-2k+1}$ for some constant $C$. We call such class the “approximate sparsity class”, motivated by the “approximate sparsity condition” from BCCHKN. This class is interesting for the following reasons. First, it generalizes various previously considered classes of functions that are characterized by their Fourier coefficients, including but not limited to the Sobolev ellipsoids and H\"{o}lder classes (e.g. Katznelson04). Second, in practice, researchers may be uncertain about the order (in terms of absolute magnitude) of the true coefficients, and our approximate sparsity class reflects such uncertainty while still maintaining the decaying properties of the re-ordered coefficients that are important for estimation purposes.

While the approximate sparsity class may seem complex, we will show that such class is sandwiched in between two classes with simpler structures, and we will use sandwich arguments to formally establish the $L^2([0,1],\mu)$ $\epsilon$-metric entropy of the approximate sparsity class and of its density subclass. To establish the upper bound on the entropy, we use an existing result from the full approximation set (lorentz66). On the other hand, we prove a lower bound using a volume-type argument inspired by Smolyak60. With these entropy bounds, we can apply the results from yangbarron99 to establish the minimax rate of convergence (in terms of MISE) for both density estimation and nonparametric regression with Gaussian noise. We obtain a minimax rate of order $n^{-(2k-1)/(2k)}$. As mentioned before, one can verify that the Sobolev ellipsoids are subsets of our approximate sparsity class, and as expected, this rate on approximate sparsity class is slower than the well-known minimax rate $n^{-2k/(2k+1)}$ for Sobolev ellipsoids. For a comprehensive review on the minimax estimation and the connections between metric entropy and minimax rates, see for example, tsybakov2008intro,yangbarron97,yangbarron99 and the references therein.

With the obtained minimax rate in mind, we propose an adaptive density estimator based on a data-driven hard thresholding procedure. The main idea is as follows. We first pick a large cutoff $J$, potentially much larger than sample size $n$, and estimate the first $J$ Fourier coefficients by the sample mean $\hat\theta_j = n^{-1}\sum_{i=1}^n\phi_j(X_i)$. However, including all the $J$ terms in the series will inevitably lead to large variance in estimation and hence suboptimal rate. To overcome this issue, in the second step we use a hard-thresholding procedure to select all the $\hat\theta_j$ above a certain data-driven threshold $\lambda$ and penalize the rest to zero. Then we form our estimator

In the final step, we project $\tilde{f}_J$ onto the space of densities to obtain a bona-fide density.

Note that since $\{\phi_j\}_{j=1}^\infty$ is an orthonormal basis, the estimation errors will be characterized entirely by Fourier coefficients $\theta_j =E\left[\phi_j(X_i)\right]$ and estimated $\hat\theta_j = n^{-1}\sum_{i=1}^n \phi_j(X_i)$. This suggests that essentially we are dealing with the problem of estimating and selecting many approximate means. Following the ideas from BCCHKN, we pick the threshold $\lambda$ to be greater than the $(1-\alpha)$-quantile of $\max_{1\leq j\leq J}|\hat\theta_j-\theta_j|$, with $\alpha\downarrow 0$ as $n\to\infty$. Previous literature suggests that such $\lambda$ can be approximated in a data-driven way using results from self-normalized moderate deviation theories or high-dimensional bootstrap, see, for example, BCCH,BCCHKN and the references therein. Although these constructions of $\lambda$ can be used to show the non-asymptotic rate, unfortunately, they are not sufficient for establishing the rate in terms of MISE. Instead, we modify their constructions and propose an alternative data-driven $\lambda$, and we show that such $\lambda$ has the desired property with the help of the Talagrand's inequality (see e.g. Bousquet2003Talagrand). We then show that if the true density $f$ belongs to an approximate sparsity class, our proposed hard-thresholding estimator is nearly minimax up to a log factor. Moreover, the estimator itself does not depend on the assumptions on the parameters of the sparsity class and is therefore adaptive.

The remaining of the paper is organized as follows. In the next section, we introduce notations and formally define the approximate sparsity class. In Section 3 we establish the metric entropy and minimax rates for density estimation and nonparametric regression with Gaussian noise in such classes. In Section 4 we elaborate on the aforementioned adaptive density estimator and we derive its rate of convergence. We then provide a specific example using the cosine basis for twice differentiable densities, in which we verify the assumptions and establish the rate of convergence by applying the results from the main theorem. We conduct simulation studies in Section 5 to illustrate the performance of our estimator, and we conclude in Section 6. The proofs are deferred to the appendix.

Approximate Sparsity Class

Suppose $\Phi:= \{\phi_j\}_{j=1}^\infty$ is an orthonormal basis of $L^2(\mathcal{X},\mu)$ for $\mathcal{X}\subset \mathbf{R}$, and without loss of generality, let $\mathcal{X} = [0,1]$. Then, for any $i\neq j$,

Moreover, for any $f\in L^2([0,1],\mu)$, there is a representation

Here $\mu$ can be either the Lebesgue measure in the context of density estimation or known probability measures on $[0,1]$ in the regression settings.

Restrictions on the Fourier coefficients $\{\theta_j\}_{j=1}^\infty$ lead to several familiar classes such as the Sobolev ellipsoids (e.g. Chapter 1.7.1 in tsybakov2008intro) and full approximation set (e.g. lorentz66). Motivated by recent literature on high dimensional models (e.g. for a comprehensive review, see the handbook chapter by BCCHKN), we introduce a new set of restrictions on the Fourier coefficients, and we call the resulting function class the approximate sparsity class:

First, we note that the re-ordering in the definition refers to the re-ordering of the coefficients by the magnitude of their absolute values. Specifically, after re-ordering, $\theta_{(j)}$ will be the j-th largest element in $\{\theta_j\}_{j=1}^\infty$ by absolute magnitude. We also want to remark that the re-ordering requirement in the definition is equivalent to that $\forall J\geq 1$, the non-increasing re-ordering of $\{\theta_{(j)}\}_{j=1}^J$ satisfies $|\theta_{(j)}|\leq Aj^{-k}$, which is a convenient characterization that we will use to establish various results in this paper. As we discussed in the introduction, it can be shown that the rate of decay of individual series coefficient is closely related to the smoothness of functions. The re-ordering condition relaxes such restriction by imposing less a priori knowledge on which series coefficients are important (as measured by magnitude).

Moreover, the restrictions on the tail sum is a natural one. In particular, if the Fourier coefficients decay at $Aj^{-k}$ without the re-ordering, the tail sum $\sum_{j=J+1}^\infty\theta_j^2$ can be shown to be bounded by $CJ^{-2k+1}$ using integral bound. However, by allowing the re-ordering of the Fourier coefficients, the tail-sum restriction is a necessary one. In particular, the tail-sum restriction imposes structures on the re-ordering, which also helps us bound the estimation bias.

At first sight, the complexity of the approximate sparsity space seems difficult to characterize. We will introduce additional spaces that have simpler structures and “sandwich” the approximate sparsity space.

The class $\mathcal{E}_k(\Phi,A)$ consists of all functions in $L^2([0,1],\mu)$ whose Fourier coefficients decay in an ordered manner in polynomial rates. For example, differentiable functions in $L^2([0,1],\mu)$ can be viewed as elements in such class. On the other hand, the set $\mathcal{A}_k(\Phi,A,C)$ is an example of the full-approximation set discussed in lorentz66. A function $f\in \mathcal{A}_k(\Phi,A,C)$ can be approximated by the partial sums $\sum_{j=1}^J\theta_j\phi_j(\cdot)$, and the tail sum restriction in the definition of $\mathcal{A}_k(\Phi,A,C)$ can be understood as the restriction on the bias from such approximation. We show in the appendix (Lemma (ref)) that for appropriately chosen constant $C$, $\mathcal{E}_k(\Phi,A)\subseteq\Theta_k(\Phi,A,C)\subseteq\mathcal{A}_k(\Phi,A,C)$. In particular, note that if $|\theta_j|<Aj^{-k}$, then

so we can simply take $C = A^2/(2k-1)$ in the definition of $\mathcal{A}_k(\Phi,A,C)$. The structures of $\mathcal{E}_k(\Phi,A)$ and $\mathcal{A}_k(\Phi,A,C)$ are much simpler and they will help us bound the metric entropy of the approximate sparsity class.

Throughout, we will use $M_2(\epsilon,\mathcal{F}):= \log N(\epsilon,\|\cdot\|_{L^2([0,1],\mu)},\mathcal{F})$ to denote the Kolmogorov $\epsilon$-entropy, where $N(\epsilon,\|\cdot\|_{L^2([0,1],\mu)}, \mathcal{F})$ is the cardinality of the largest $\epsilon$-packing set of a set $\mathcal{F}$ under the $L^2([0,1],\mu)$ distance. Next, we formally introduce the definition of minimax rates, borrowing notation from tsybakov2008intro. Let $\{a_n\}$ and $\{b_n>0\}$ be real sequences. We write $a_n\gtrsim b_n$ if $\liminf_{n\to \infty} a_n/b_n >0$, and similarly, we write $a_n\lesssim b_n$ if $\limsup_{n\to \infty} a_n/b_n <\infty$. We use the notation “$\asymp$” as the asymptotic order symbol. That is, we write $a_n\asymp b_n$ if

Entropy and Minimax Rate for Approximate Sparsity Class

As we have discussed in the previous section, for appropriately chosen constant $C$ in the definitions of $\Theta_k(\Phi,A,C)$ and $\mathcal{A}_k(\Phi,A,C)$, we have $\mathcal{E}_k(\Phi,A)\subseteq\Theta_k(\Phi,A,C)\subseteq\mathcal{A}_k(\Phi,A,C)$. Therefore, understanding the metric entropy of $\mathcal{E}_k(\Phi,A)$ and $\mathcal{A}_k(\Phi,A,C)$ can help us control the metric entropy of the approximate sparsity class $\Theta_k(\Phi,A,C)$. In particular, the entropy on $\mathcal{E}_k(\Phi,A)$ and $\mathcal{A}_k(\Phi,A,C)$ will respectively give lower and upper bounds on the entropy of the sandwiched set $\Theta_k(\Phi,A,C)$, as shown in the next lemma.

The proof of this lemma is given in the appendix. As we remarked before, $\mathcal{A}_k(\Phi,A,C)$ is a special case of the full approximation set introduced in lorentz66, and its entropy can be established using Theorem 3 in lorentz66. On the other hand, to the best of our knowledge, the earliest reference of the class $\mathcal{E}_k(\Phi,A)$ can be traced back to Smolyak60 in which the trigonometric basis $\Phi$ is considered \footnote{The class $\mathcal{E}_k(\Phi,A)$ is sometimes referred to as the “hyperrectangles” in the literature, and minimax risks over these hyperrectangles has been studied in the context of Gaussian shift models, see, for example, donoho90 and the references therein.}. Inspired by Smolyak60, we adapt their arguments and extend similar results to $\mathcal{E}_k(\Phi,A)$ for more general orthonormal bases \footnote{Prior to learning the existence of Smolyak60, we pursued an alternative route when establishing the entropy of $\mathcal{E}_k(\Phi,A)$. This alternative proof relies on an isometry between $\mathcal{E}_k(\Phi,A)$ and a generalized Hilbert cube (kloeckner12). kloeckner12 establishes the bounds on the entropy of the Hilbert cubes with which one can infer lower bounds on the entropy of $\mathcal{E}_k(\Phi,A)$. While these lower bounds are sufficient for establishing the minimax rates, we opt to adopt the arguments used by Smolyak60 for the sake of clarity.}. Next, we use Lemma (ref) to establish the following theorem, which can be shown using a sandwich type of argument.

This theorem formally establishes bounds on the $L^2$ metric entropy of our approximate sparsity class $\Theta_k(\Phi,A,C))$with the help of $\mathcal{E}_k(\Phi,A)$ and $\mathcal{A}_k(\Phi,A,C)$. The following corollary establishes a similar entropy result on the density subset of $\Theta_k(\Phi,A,C))$.

We note that many familiar orthonormal bases on $L^2([0,1],\mu)$ contain a constant term. Moreover, the requirement that the functions in $\Theta_k(\Phi,A,C)$ are uniformly bounded, can be satisfied, for example, when the orthonormal basis $\Phi$ is uniformly bounded and $k>1$. Assuming uniform boundedness will allow us to find a set of the densities of the form $(f+M'+1)/(\int f d\mu + M'+1)$ for $f\in\mathcal{E}_k(\Phi,A)$ for some large constant $M'$, which provides a lower bound for the entropy of $\tilde{\Theta}_k(\Phi,A',C')$ for some constants $A'$ and $C'$, which in turn is a lower bound for $\tilde{\Theta}_k(\Phi,A,C)$. Then applying the sandwich argument again gives us the results of the corollary. With the entropy results, we formally establish the minimax rates for nonparametric regression with Gaussian noise and density estimation on the approximate sparsity classes.

Theorem (ref) can be seen as a direct consequence of Theorem (ref) and Corollary (ref) due to the tight connection between the entropy and minimax rates established in yangbarron99 and the references therein. In particular, the minimax rate of estimation on a particular class of functions $\mathcal{F}$ is closely related to the “critical separation” $\epsilon_n$, which is determined via the identity $M_2(\epsilon_n,\mathcal{F}) = n\epsilon_n^2$. Then the results of the theorem should follow from our entropy results. We note that there are some caveats when applying yangbarron99 directly on the density subset $\tilde\Theta_k(\Phi,A,C)$ since it is not convex nor bounded away from zero. However, we show in the proof that this subset is sandwiched in between two convex sets of densities with which we can establish the desired result.

Nearly Minimax Optimal Adaptive Density Estimator

Preliminaries

In this section, we propose a density estimator that is adaptive and achieves the minimax rate of convergence on $\tilde\Theta_k(\Phi,A,C)$ up to a log term. Suppose $X$ is a random variable with probability density $f\in L^2([0,1],\mu)$ and let $\Phi = \{\phi_j\}_{j=1}^\infty$ be an orthonormal basis of $L^2([0,1],\mu)$ with $\mu$ being the Lebesgue measure. Then $f$ has a unique representation: \[ f(\cdot) = \sum_{j=1}^\infty\theta_j\phi_j(\cdot). \] Moreover, $f$ being a probability density and $\Phi = \{\phi_j\}_{j=1}^\infty$ being an orthonormal basis allow us to find expressions of $\theta_j$'s in terms of expectations:

Given an i.i.d sample $\{X_i\}_{i=1}^n$ with $X_i\sim X$, the standard series estimator for the density $f$ builds on identity ((ref)) and takes the form \[ \hat{f}_{J}(x) = \sum_{j=1}^J\hat\theta_j\phi_j(x)\quad\text{where}\quad\hat\theta_j = \frac{1}{n}\sum_{i=1}^n\phi_j(X_i) \] Note that the mean integrated squared error (MISE) of $\hat{f}_n$ equals

The choice of the cutoff $J$ plays an essential role in the trade-off between variance and bias in ((ref)), and yet in order to properly choose $J$ in estimation, one often has to assume some prior knowledge about the smoothness of the unknown $f$.

We propose an alternative estimator that circumvents this issue. In particular, our estimator requires choosing a large cutoff $J$, and then uses LASSO to select the most “relevant” Fourier coefficients among $\{\hat\theta_1,\hat\theta_2,\cdots,\hat\theta_J\}$. Note that many of the familiar bases, such as cosine and Legendre polynomial basis, contain a constant element. Therefore, we always assume the orthonormal basis of choice contains a constant term, which, without loss of generality, we assume to be the first basis term $\phi_1$. Then note that $\hat\theta_1 = n^{-1}\sum_{i=1}^n\phi_1 = \phi_1 = E\left[\phi_1\right] = \theta_1$, that is, there's no estimation error for the first term. Therefore, we should always include $\hat\theta_1$ in our estimator in practice.

For a slight change of notation, let $\theta^J = \{\theta_1,\theta_2,\cdots,\theta_J\}\in\mathbf{R}^{J}$ denote the first $J$ true but unknown series coefficients. Let $\hat\theta^J\in\mathbf{R}^{J}$ be an estimator of $\theta^J$ defined as follows

where each $\hat\theta_j$ is consistent by weak law of large numbers under mild regularity conditions, and $\hat\theta^J$ can be shown to be consistent using Bernstein's inequality and maximal inequality under additional assumptions (e.g. if the orthonormal basis $\Phi$ is bounded or if the density is bounded). Given $\hat\theta^J$, the estimator we consider here is the so-called hard-thresholding estimator:

where $\omega_j$'s are the thresholding parameters depending on the penalty parameter $\lambda$. To be consistent with the notation in ((ref)), we let $\omega_j = 0$ for $j>J$. The hard thresholding estimator we consider here differs from the soft-thresholding estimator that has been considered in some previous literature, for example, bunea2010spades and the references therein. Intuitively, for a properly chosen penalty parameter $\lambda$, we penalize “small” estimates of the series coefficients to $0$ while keeping the rest unchanged. Let $T\subseteq \{1,\cdots,J\}$ denote the set of selected indices, that is,

As a result, the new estimator we get is

As we will show later, the quality of this estimator depends on the penalty parameter and, to a lesser degree, the cutoff $J$.

Moreover, note that $\tilde{f}_J$ is not necessarily a probability density. We need to make sure $\tilde{f}_J$ integrates to $1$ and is nonnegative. The former is easily satisfied if $\mu$ is the Lebesgue measure on $[0,1]$ and the first basis element $\phi_1 = 1$, which, by the definition of orthonormality, implies

This is another reason why we should always include $\hat\theta_1 = \phi_1$ in our estimator.

On the other hand, some forms of post-processing are needed to ensure that the estimator is nonnegative. We will follow the P-Algorithm suggested by Gajek86, which is attractive in both its ease of implementation and the statistical properties of the resulting estimator. Here we review the P-Algorithm and illustrate with our estimator:

Gajek86 shows that the P-Algorithm converges both point-wise and in $\|\cdot\|_{L^2([0,1],\mu)}$ to a density $\hat{f}^* = \max\{0,\tilde{f}_J + c\}$, for some constant $c$. Moreover, $\hat{f}^*$ is the orthogonal projection of $\tilde{f}_J$ onto the space of densities, and $\hat{f}^*$ has at least the same rate of convergence (as measured in MISE) as $\tilde{f}_J$.

Main results

To estimate $\hat{f}^*$ and to establish its statistical properties, we need to consider a data-driven way of choosing the regularization parameter $\lambda$ in ((ref)). In particular, in order to control the penalization error in a uniform manner, we want to pick regularization parameter $\lambda$ in a way such that

There are several ways of achieving this. For example, one can approximate such $\lambda$ using the results from the self-normalized moderate deviation theory or high dimensional bootstrap literature, see for example BCCHKN and the references therein. However, in order to establish the MISE, we need ((ref)) to hold with $\alpha = \alpha_n$ going zero sufficiently fast as sample size $n$ increases. To this end, we will modify the $\lambda$ proposed by BCCHKN based on the moderation deviation theories, and instead we will use the Talagrand's inequality (see Bousquet2003Talagrand) to establish the asymptotic result.

To facilitate the discussion, we borrow the following notations from BCCHKN: Let $\{X_i\}_{i=1}^n$ be an i.i.d random sample in $\mathbf{R}$ and let $\{\phi_j\}_{j=1}^\infty$ be an orthonormal basis of $L^2([0,1],\mu)$. With some abuse of notation, let $\Phi$ denote the CDF of the standard normal variable, and define

and we let $\hat{Z}_{ij}$ be the sample analog of $Z_{ij}$

We introduce the following regularity conditions.

Condition (ref) will be assumed to establish the main results in this section. We want to emphasize that all but (i) in Condition (ref) can be verified, and we note that (i) is stated in the form that is convenient for the proof and the lower bound $1$ in the statement can be replaced by any positive constant. In particular, (ii) is chosen by researchers, and (iii) can be checked for a given orthonormal basis. Note that when the chosen orthonormal basis is uniformly bounded, such as the cosine basis, we have $M_J^2 = M$ for some constant $M$, in which case (iii) is trivially satisfied, and we can potentially choose $J\gg n$. On the other hand, if we have a growing basis, the choices of $J$ can be limited depending on how fast the basis grows. For example, the Legendre basis grows with $M_J^2 = 2J+1$, and we have to choose $J$ such that $J = n^p$ for some $p<1$ ($p$ can be close to $1$ for $n$ large). The requirement (iv) can be viewed as the extra cost we pay when bounding the variance term in the MISE. In the next theorem, we propose a $\lambda$ such that $\lambda\geq (1-\alpha_n)$-quantile of $\|\hat{\theta}^J-\theta^J\|_\infty$ with $\alpha_n$ goes to zero sufficiently fast so that the MISE can be established.

The proof of the theorem is given in the appendix and it can be easily adapted to allow other choices of $\alpha_n$ so that (iv) in Condition (ref) can be satisfied. The expression of $\lambda$ in the theorem bears some similarities to the ones given in BCCH and in BCCHKN (Theorem 2.4) using moderate deviation theories. Compared to their constructions, for $J = n^p$, our $\lambda$ has an extra $\sqrt{\log(n)}$ multiplied (ignoring the constants). Roughly speaking, by the Talagrand's inequality, this extra $\log$ term pushes our $\lambda$ into the exponential tail and allows us to obtain $\alpha_n$ that is sufficiently fast for establishing MISE. We remark that the $\lambda$ based on the moderate deviation theories can still be used to establish the non-asymptotic rate.

With this result, we can establish the second main result of this section. In particular, we are going to assume that the density of the random variable in question belongs to the approximate sparsity space $\tilde\Theta_k(\Phi,A,C)$, and we will show that the post-processed estimator $\hat{f}^*$ defined in (ref) admits the minimax rate of convergence on $\tilde\Theta_k(\Phi,A,C)$ up to a log term.

The proof of the theorem is given in the appendix. To bound the MISE, we first decompose the MISE into roughly the standard “variance” and “bias” components, and then we bound each separately with the help of Theorem (ref). In the assumptions of the theorem, $f\in \tilde\Theta_k(\Phi,A,C)$ imposes restrictions on the Fourier coefficients of $f$, which will play a crucial role in establishing the rates. Moreover, we assume that $\tilde\Theta_k(\Phi,A,C)$ is uniformly bounded, which is also assumed when we establish the minimax rates. Note that for the cases when $k>1$ and the orthonormal basis $\Phi$ is uniformly bounded, one can verify that $\tilde\Theta_k(\Phi,A,C)$ is uniformly bounded. In addition, Condition (ref) is assumed so that we can utilize the results in Theorem (ref), and as we commented before, the requirements in Condition (ref) are easy to verify.

We also want to emphasize the adaptive nature of our estimator. In particular, once the researchers have decided on which orthonormal basis to use, the construction of the estimator $\hat{f}^*$ with data-driven $\lambda$ does not depend on the assumptions on the approximate sparsity class (e.g. how fast the Fourier coefficients decay). Theorem (ref) simply states that our estimator achieves the minimax rate on any approximate sparsity class. This is attractive in practice since the researchers do not have to assume the smoothness of the true data-generating density other than that it belongs to some approximate sparsity class.

Example Using Cosine Basis

In this section, we are going to illustrate our previous results with the differentiable densities and the cosine basis. The standard cosine orthonormal basis $\{\phi_j\}_{j=1}^\infty$ on $L^2([0,1],\mu)$ is defined as follows:

which is uniformly bounded. Suppose $f\in L^2([0,1],\mu)$ is twice differentiable. Then for $\{\phi_j\}_{j=1}^\infty$ being the cosine basis defined above, $f(\cdot) = \sum_{j=1}^\infty\theta_j\phi_j(\cdot)$ and there exists some constant $c$ such that \[ |\theta_j|\leq cj^{-2}\int_0^1|f^{(2)}(x)|~dx,\quad j\geq 1 \] As noted in efro2008 (Section 2.2), unless additional (boundary) assumptions are made on the function $f$, the series coefficients can not decay faster than the rate $j^{-2}$ (regardless of the smoothness of functions). This makes the set of twice differentiable functions with bounded second derivative a special example of the approximate sparsity class. In fact, this is an example of $\mathcal{E}_k(\Phi,A)$, which are themselves a special case of the approximate sparsity class without reordering. The next result follows directly from Theorem (ref).

Note that since the cosine basis is uniformly bounded, we can drop the assumption that the density is bounded and the requirements in Condition (ref) can be verified. Moreover, in view of our results in Section 3, the rate in this corollary is minimax up to a log term. As the approximate sparsity class and $\mathcal{E}_k(\Phi,A)$ class are new (and we commented previously that the familiar Sobolev ellipsoids are subsets of these classes), we have not yet considered the performances of other density estimators on such classes, and whether there are other adaptive rate-optimal estimators on these new classes remains an open question.

Simulations

In this section, we conduct simulation studies in which we compare the simulated MISE of our estimator introduced in Section 4 with an alternative estimator for which the series cutoff has to be properly specified. The true data generating density is constructed to be in the approximate sparsity class.

Design density in approximate sparsity class

For the first simulation, we construct a density whose (cosine) Fourier coefficients satisfy approximate sparsity. Let

where $\phi_j$'s are the cosine basis terms defined in ((ref)). We specify the $\theta_j$'s in the following way:

As we've shown in section 4, since $\theta_1 = 1$, $f_X$ integrates to 1. The constant $A = 2$ is chosen such that $f_X$ in ((ref)) is non-negative and hence a proper probability density.

Simulation procedures

To simulate MISE, we proceed in following steps:

Step 1: We draw $B = 1000$ independent i.i.d. samples $\{X_i\}_{i=1}^{N}$ of size $N$. Here the $N$ we consider will be $N = 5000, 10000, 15000, 20000$.

Step 2(a): For each sample $\{X_i\}_{i=1}^{N}$ from the true density $f_X$, we construct an estimator $\hat{f}^*$ described in section 4:

Step 2(b): For each sample $\{X_i\}_{i=1}^{N}$ from the true density $f_X$, we construct an alternative comparison estimator $\check{f}$ in the following way:

Step 3: For each $1\leq b\leq B$ sample defined in Step 1 and associated estimator from Step 2(a), we estimate the integrated squared error (ISE) using numerical integration, and we use $\widehat{ISE}_b$ to denote the approximated ISE. We then calculate the estimated MISE

Repeat the same process for the estimator from Step 2(b).

Simulation Results

The simulation design described above is coded directly using Python, with code available upon request. Since the design density is not a known density coded in Python packages, we use inverse transform sampling to generate random samples from our design density. We illustrate the performance of the inverse transform sampling in Figure ((ref)), where the solid line is the design density, and the normalized histogram below the solid line is constructed using a random sample of size $10000$ from such inverse transform sampling.

We present the simulation results in Figure ((ref)). Each dot represents the simulated MISE for a given sample size, where the red dots (labeled as “f_star”) are the simulated MISEs for our estimator and the blue dots (labeled as “f_check”) are for the comparison estimator. We make the following observations. First, as expected, for both our estimator and the comparison estimator, the simulated MISE decreases as sample size increases. For the comparison estimator, this reflects the fact that such estimator is using larger series cutoffs as sample size increases, which will eventually include all the nonzero terms of the series coefficients. On the other hand, our estimator requires choosing a large series cutoff ($J=200$) to start with and then uses a data-driven hard-thresholding procedure to select the relevant terms. Second, for each of the sample sizes in our consideration, the simulated MISE of our estimator is smaller than that of the comparison estimator. This is likely due to the special feature of our design density, where the large Fourier coefficients show up in later series terms. Our estimator estimates first $J=200$ series terms and the data-driven hard-thresholding procedure is able to pick up the large series terms to reduce the bias. On the other hand, the comparison estimator has to specify the series cutoff $N^{1/4}$ and may fail to include the large Fourier coefficients that show up in later terms.

These simulation results demonstrate that our estimator performs well (as measured by the simulated MISE) for estimating densities in the approximate sparsity class. Moreover, the results showcase the adaptive nature of our estimator in comparison to an alternative estimator for which the researcher has to make potentially restrictive smoothness assumptions in order to determine the proper series cutoff. Even then, if the sample size is not large enough, such comparison estimator may still miss large Fourier coefficients that show up in later series terms, which can result in larger estimation errors.

Concluding Remarks

In this paper, we have studied a new class of functions, which we call the approximate sparsity class, that is characterized by a new set of restrictions on the Fourier coefficients of those functions for a given orthonormal basis. We have derived the upper and lower bounds on $L^2$ $\epsilon$-metric entropy of the approximate sparsity class and we have established the minimax rates for nonparametric regression with Gaussian noise and for density estimation. We have shown that functions in such classes are natural candidates for the thresholding types of estimators and we proposed an adaptive density estimator that is nearly minimax optimal (up to a log term) over such class. For future research, we hope to study estimators for nonparametric regression for functions in the approximate sparsity class and we hope to generalize the approximate sparsity class to high-dimensional regression settings.

Acknowledgements

The author would like to express his gratitude to Andres Santos and Denis Chetverikov for their generous time and extremely helpful discussions that have lead to substantial improvements to this project. The author would also like to thank Rosa Matzkin, Oscar H. Madrid-Padilla, Mingli Chen, Rodrigo Pinto, and participants at UCLA econometrics proseminars for their helpful suggestions.