Sparse Additive Models

Introduction

The nonparametric regression model $Y_{i}=m(X_{i})+\epsilon_{i}$, where $m$ is a general smooth function, relaxes the strong assumptions made by a linear model, but is much more challenging in high dimensions. Hastie and Tibshirani (1999) introduced the class of additive models of the form

This additive combination of univariate functions—one for each covariate $X_{j}$—is less general than joint multivariate nonparametric models, but can be more interpretable and easier to fit; in particular, an additive model can be estimated using a coordinate descent Gauss-Seidel procedure, called backfitting. Unfortunately, additive models only have good statistical and computational behavior when the number of variables $p$ is not large relative to the sample size $n$, so their usefulness is limited in the high dimensional setting.

In this paper we investigate sparse additive models (SpAM), which extend the advantages of sparse linear models to the additive, nonparametric setting. The underlying model is the same as in (1), but we impose a sparsity constraint on the index set $\{j:f_{j}\not\equiv 0\}$ of functions $f_{j}$ that are not identically zero. Lin and Zhang (2006) have proposed COSSO, an extension of lasso to this setting, for the case where the component functions $f_{j}$ belong to a reproducing kernel Hilbert space (RKHS). They penalize the sum of the RKHS norms of the component functions. Yuan (2007) proposed an extension of the non-negative garrote to this setting. As with the parametric non-negative garrote, the success of this method depends on the initial estimates of component functions $f_{j}$.

In Section 3, we formulate an optimization problem in the population setting that induces sparsity. Then we derive a sample version of the solution. The SpAM estimation procedure we introduce allows the use of arbitrary nonparametric smoothing techniques, effectively resulting in a combination of the lasso and backfitting. The algorithm extends to classification problems using generalized additive models. As we explain later, SpAM can also be thought of as a functional version of the grouped lasso (Yuan and Lin, 2006).

The main results of this paper include the formulation of a convex optimization problem for estimating a sparse additive model, an efficient backfitting algorithm for constructing the estimator, and theoretical results that analyze the effectiveness of the estimator in the high dimensional setting. Our theoretical results are of several different types. First, we show that, under suitable choices of the design parameters, the SpAM backfitting algorithm recovers the correct sparsity pattern asymptotically; this is a property we call sparsistence, as a shorthand for “sparsity pattern consistency.” Second, we show that that the estimator is persistent, in the sense of Greenshtein and Ritov (2004), which is a form of risk consistency. Specifically, we show:

Here $S=\{j\,:\,f_{j}\neq 0\}$ is the index set for the nonzero components, $\widehat{S}=\{j\,:\,\widehat{f}_{j}\neq 0\}$ and ${\mathcal{M}}_{n}$ is a class of functions defined by the level of regularization.

In the following section we establish notation and assumptions. In Section 3 we formulate SpAM as an optimization problem and derive a scalable backfitting algorithm. An extension to sparse nonparametric logistic regression is presented in Section 4. Examples showing the use of our sparse backfitting estimator on high dimensional data are included in Section 6. In Section 7.1 we formulate the sparsistency result, when orthogonal function regression is used for smoothing. In Section 7.2 we give the persistence result. Section 8 contains a discussion of the results and possible extensions. Proofs are contained in Section 9.

Notation and Assumptions

We assume that we are given data $(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$ where $X_{i}=(X_{i1},\ldots,X_{ij},\ldots,X_{ip})^{T}\in^{p}$ and

with $\epsilon_{i}\sim N(0,\sigma^{2})$ and

Denote the joint distribution of $(X_{i},Y_{i})$ by $P$. For a function $f$ on $$denote its$L_{2}(P)$ norm by

Let $\{\psi_{jk},k=0,1,\ldots\}$ denote a uniformly bounded, orthonormal basis with respect to Lebesgue measure on $$. Unless stated otherwise, we assume that$f_{j}\in{\cal T}_{j}$ where

for some $0<C<\infty$. We shall take $\nu_{j}=2$ although the extension to other levels of smoothness is straightforward. It is also possible to adapt to $\nu_{j}$ although we do not pursue that direction here.

Let $\Lambda_{\rm min}(A)$ and $\Lambda_{\rm max}(A)$ denote the minimum and maximum eigenvalues of a square matrix $A$. If $v=(v_{1},\ldots,v_{k})^{T}$ is a vector, we use the norms

Sparse Backfitting

The outline of the derivation of our algorithm is as follows. We first formulate a population level optimization problem, and show that the minimizing functions can be obtained by iterating through a series of soft-thresholded univariate conditional expectations. We then plug in smoothed estimates of these univariate conditional expectations, to derive our sparse backfitting algorithm.

where the expectation is taken with respect to $X$ and the noise $\epsilon$. Now consider the following modification of this problem that introduces a scaling parameter for each function, and that imposes additional constraints:

It is convenient to re-express the minimization in the following equivalent form:

The optimization problem in (12) can also be written in the penalized Lagrangian form,

The minimizers $f_{j}\in{\mathcal{H}}_{j}$ of (14) satisfy

At the population level, the $f_{j}$s can be found by a coordinate descent procedure that fixes $(f_{k}:\ k\neq j)$ and fits $f_{j}$ by equation (15), then iterates over $j$.

where ${\mathcal{S}}_{j}$ is a linear smoother, such as a local linear or kernel smoother. Let

This algorithm can be seen as a functional version of the coordinate descent algorithm for solving the lasso. In particular, if we solve the lasso by iteratively minimizing with respect to a single coordinate, each iteration is given by soft thresholding; see Figure 2. Convergence properties of variants of this simple algorithm have been recently treated by Daubechies et al. (2004, 2007). Our sparse backfitting algorithm is a direct generalization of this algorithm, and it reduces to it in case where the smoothers are local linear smoothers with large bandwidths.

Basis Functions. It is useful to express the model in terms of basis functions. Recall that $B_{j}=(\psi_{jk}:\ k=1,2,\ldots)$ is an orthonormal basis for ${\cal T}_{j}$ and that $\sup_{x}|\psi_{jk}(x)|\leq B$ for some $B$. Then

where $\beta_{jk}=\int f_{j}(x_{j})\psi_{jk}(x_{j})dx_{j}$.

where $d=d_{n}$ is a truncation parameter. For the Sobolev space ${\cal T}_{j}$ of order two we have that $\left\|f_{j}-\widetilde{f}_{j}\right\|^{2}=O(1/d^{4})$. Let $S=\{j:\ f_{j}\neq 0\}$. Assuming the sparsity condition $|S|=O(1)$ it follows that $\left\|m-\widetilde{m}\right\|^{2}=O(1/d^{4})$ where $\widetilde{m}=\sum_{j}\widetilde{f}_{j}$. The usual choice is $d\asymp n^{1/5}$ yielding truncation bias $\left\|m-\widetilde{m}\right\|^{2}=O(n^{-4/5})$.

which is the sample version of (14). In Section 7.1 we prove theoretical properties assuming that this particular smoother is being used.

Connection with the Grouped Lasso. The SpAM model can be thought of as a functional version of the grouped lasso (Yuan and Lin, 2006) as we now explain. Consider the following linear regression model with multiple factors,

where $Y$ is an $n\times 1$ response vector, $\epsilon$ is an $n\times 1$ vector of iid mean zero noise, $X_{j}$ is an $n\times d_{j}$ matrix corresponding to the $j$-th factor, and $\beta_{j}$ is the corresponding $d_{j}\times 1$ coefficient vector. Assume for convenience (in this subsection only) that each $X_{j}$ is orthogonal, so that $X_{j}^{T}X_{j}=I_{d_{j}}$, where $I_{d_{j}}$ is the $d_{j}\times d_{j}$ identity matrix. We use $X=(X_{1},\ldots,X_{p_{n}})$ to denote the full design matrix and use $\beta=(\beta_{1}^{T},\dots,\beta_{p_{n}}^{T})^{T}$ to denote the parameter.

The grouped lasso estimator is defined as the solution of the following convex optimization problem:

where $\sqrt{d_{j}}$ scales the $j$th term to compensate for different group sizes.

It is obvious that when $d_{j}=1$ for $j=1,\ldots,p_{n}$, the grouped lasso becomes the standard lasso. From the KKT optimality conditions, a necessary and sufficient condition for $\widehat{\beta}=(\widehat{\beta}_{1}^{T},\ldots,\widehat{\beta}_{p}^{T})^{T}$ to be the grouped lasso solution is

Based on this stationary condition, an iterative blockwise coordinate descent algorithm can be derived; as shown by Yuan and Lin (2006), a solution to (23) satisfies

where $S_{j}=X^{T}_{j}(Y-X\beta_{\backslash j})$, with $\beta_{\backslash j}=(\beta^{T}_{1},\ldots,\beta^{T}_{j-1},\mathbf{0}^{T},\beta^{T}_{j+1},\ldots,\beta^{T}_{p_{n}})$. By iteratively applying (24), the grouped lasso solution can be obtained.

As discussed in the introduction, the COSSO model of Lin and Zhang (2006) replaces the lasso constraint on $\sum_{j}|\beta_{j}|$ with a RKHS constraint. The advantage of our formulation is that it decouples smoothness ($g_{j}\in{\cal T}_{j}$) and sparsity ($\sum_{j}|\beta_{j}|\leq L$). This leads to a simple algorithm that can be carried out with any nonparametric smoother and scales easily to high dimensions.

Sparse Nonparametric Logistic Regression

The SpAM backfitting procedure can be extended to nonparametric logistic regression for classification. The additive logistic model is

where $Y\in\{0,1\}$, and the population log-likelihood is

Recall that in the local scoring algorithm for generalized additive models (Hastie and Tibshirani, 1999) in the logistic case, one runs the backfitting procedure within Newton’s method. Here one iteratively computes the transformed response for the current estimate $f_{0}$

and weights $w(X_{i})=p(X_{i};f_{0})(1-p(X_{i};f_{0})$, and carries out a weighted backfitting of $(Z,X)$ with weights $w$. The weighted smooth is given by

To incorporate the sparsity penalty, we first note that the Lagrangian is given by

In the finite sample case, in terms of the smoothing matrix ${\mathcal{S}}_{j}$, this becomes

If $\|{\mathcal{S}}_{j}(wR_{j})\|<\lambda$, then $f_{j}=0$. Otherwise, this implicit, nonlinear equation for $f_{j}$ cannot be solved explicitly, so we propose to iterate until convergence:

When $\lambda=0$, this yields the standard local scoring update (28). An example of logistic SpAM is given in Section 6.

Choosing the Regularization Parameter

We choose $\lambda$ by minimizing an estimate of the risk. Let $\nu_{j}$ be the effective degrees of freedom for the smoother on the $j^{\rm th}$ variable, that is, $\nu_{j}={\rm trace}({\cal S}_{j})$ where ${\mathcal{S}}_{j}$ is the smoothing matrix for the $j$-th dimension. Also let $\widehat{\sigma}^{2}$ be an estimate of the variance. Define the total effective degrees of freedom as

The first is $C_{p}$ and the second is generalized cross validation but with degrees of freedom defined by $\mathop{\text{df}}(\lambda)$. A proof that these are valid estimates of risk is not currently available; thus, these should be regarded as heuristics.

Based on the results in Wasserman and Roeder (2007) about the lasso, it seems likely that choosing $\lambda$ by risk estimation can lead to overfitting. One can further clean the estimate by testing $H_{0}:f_{j}=0$ for all $j$ such that $\widehat{f}_{j}\neq 0$. For example, the tests in Fan and Jiang (2005) could be used.

Examples

To illustrate the method, we consider a few examples.

Synthetic Data. Our first example is from (Härdle et al., 2004). We generated $n=150$ observations from the following 200-dimensional additive model:

and $f_{j}(x)=0$ for $j\geq 5$ with noise $\epsilon_{i}~{}\sim\mathcal{N}(0,1)$. These data therefore have 196 irrelevant dimensions.

The results of applying SpAM with the plug-in bandwidths are summarized in Figure 3. The top-left plot in Figure 3 shows regularization paths as the parameter $\lambda$ varies; each curve is a plot of $\|\widehat{f}_{j}(\lambda)\|$ versus

for a particular variable $X_{j}$. The estimates are generated efficiently over a sequence of $\lambda$ values by “warm starting” $\widehat{f}_{j}(\lambda_{t})$ at the previous value $\widehat{f}_{j}(\lambda_{t-1})$. The top-center plot shows the $C_{p}$ statistic as a function of $\lambda$. The top-right plot compares the empirical probability of correctly selecting the true four variables as a function of sample size $n$, for $p=128$ and $p=256$. This behavior suggests the same threshold phenomenon that was shown for the lasso by Wainwright (2006).

Boston Housing. The Boston housing data were collected to study house values in the suburbs of Boston. There are 506 observations with 10 covariates. The dataset has been studied by many other authors (Härdle et al., 2004; Lin and Zhang, 2006), with various transformations proposed for different covariates. To explore the sparsistency properties of our method, we add 20 irrelevant variables. Ten of them are randomly drawn from $\text{Uniform}(0,1)$, the remaining ten are a random permutation of the original ten covariates. The model is

The result of applying SpAM to this 30 dimensional dataset is shown in Figure 4. SpAM identifies 6 nonzero components. It correctly zeros out both types of irrelevant variables. From the full solution path, the important variables are seen to be ${\tt rm}$, lstat, ptratio, and crim. The importance of variables nox and b is borderline. These results are basically consistent with those obtained by other authors (Härdle et al., 2004). However, using $C_{p}$ as the selection criterion, the variables ${\tt indux}$, ${\tt age}$, ${\tt dist}$, and ${\tt tax}$ are estimated to be irrelevant, a result not seen in other studies.

SpAM for Spam. Here we consider an email spam classification problem, using the logistic SpAM backfitting algorithm from Section 3. This dataset has been studied by Hastie et al. (2001), using a set of 3,065 emails as a training set, and conducting hypothesis tests to choose significant variables; there are a total of 4,601 observations with $p=57$ attributes, all numeric. The attributes measure the percentage of specific words or characters in the email, the average and maximum run lengths of upper case letters, and the total number of such letters. To demonstrate how SpAM performs with sparse data, we only sample $n=300$ emails as the training set, with the remaining $4301$ data points used as the test set. We also use the test data as the hold-out set to tune the penalization parameter $\lambda$.

The results of a typical run of logistic SpAM are summarized in Figure 5, using plug-in bandwidths. It is interesting to note that even with this relatively small sample size, logistic SpAM recovers a sparsity pattern that is consistent with previous authors’ results. For example, in the best model chosen by logistic SpAM, according to error rate, the 33 selected variables cover 80% of the significant predictors as determined by Hastie et al. (2001).

Functional Sparse Coding. Olshausen and Field (1996) propose a method of obtaining sparse representations of data such as natural images; the motivation comes from trying to understand principles of neural coding. In this example we suggest a nonparametric form of sparse coding.

This optimization problem is not jointly convex in $\beta^{(i)}$ and $X$. However, for fixed $X$, each weight vector $\beta^{(i)}$ is computed by running the lasso. For fixed $\beta^{(i)}$, the optimization is similar to ridge regression, and can be solved efficiently. Thus, an iterative procedure for (approximately) solving this optimization problem is easy to derive.

In the case of sparse coding of natural images, as in Olshausen and Field (1996), the basis vectors $X_{j}$ encode basic edge features. A code with 200 basis vectors, estimated by carrying out the optimization using the lasso and stochastic gradient descent, is shown in Figure 6. The codewords are seen to capture edge features at different scales and spatial orientations.

In the functional version, we no longer assume a linear, parametric fit between the dictionary $X$ and the data $y$. Instead, we model the relationship using an additive model:

Figures 7 and 8 illustrate the reconstruction of different image patches using the sparse linear model compared with the sparse additive model. The codewords $X_{j}$ are those obtained using the Olshausen-Field procedure; these become the design points in the regression estimators. Thus, a codeword for a $16\times 16$ patch corresponds to a vector $X_{j}$ of dimension $256$, with each $X_{ij}$ the gray level for a particular pixel.

It can be seen how the functional fit achieves a more accurate approximation to the image patch, with fewer codewords. Each set of plots shows the original image patch, the reconstruction, and the codewords that were used in the reconstruction, together with their marginal fits to the data. For instance, in Figure 7 it is seen that the sparse linear model uses eight codewords, with a residual sum of squares (RSS) of 0.0561, while the sparse additive model uses seven codewords to achieve a residual sum of squares of 0.0206. It can also be seen that the linear and nonlinear fits use different sets of codewords. Local linear smoothing was used with a Gaussian kernel having fixed bandwidth $h=0.05$ for all patches and all codewords.

These results are obtained using the set of codewords obtained under the sparse linear model. The codewords can also be learned using the sparse additive model; this will be reported in a separate paper.

Theoretical Properties

In the case of linear regression, with $f_{j}(X_{j})=\beta_{j}^{*T}X_{j}$, several authors have shown that, under certain conditions on $n$, $p$, the number of relevant variables $s=|{\text{supp}}(\beta^{*})|$, and the design matrix $X$, the lasso recovers the sparsity pattern asymptotically; that is, the lasso estimator $\widehat{\beta}_{n}$ is sparsistent:

Here, ${\text{supp}}(\beta)=\left\{j:\ \beta_{j}\neq{0}\right\}$. References include Wainwright (2006), Meinshausen and P. Bühlmann (2006), Zou (2005) and Zhao and Yu (2007). We show a similar result for sparse additive models under orthogonal function regression.

In terms of an orthogonal basis $\psi$, we can write

To simplify notation, let $\beta_{j}$ be the $d_{n}$ dimensional vector $\{\beta_{jk},\,k=1,\ldots,d_{n}\}$ and let $\Psi_{j}$ be the $n\times d_{n}$ matrix $\Psi_{j}[i,k]=\psi_{jk}(X_{ij})$. If $A\subset\{1,\ldots,p\}$, we denote by $\Psi_{A}$ the $n\times d|A|$ matrix where for each $j\in A$, $\Psi_{j}$ appears as a submatrix in the natural way.

We now analyze the sparse backfitting Algorithm (1) assuming an orthogonal series smoother is used to estimate the conditional expectation in its Step (2). As noted earlier, an orthogonal series smoother for a predictor $X_{j}$ is the least squares projection onto a truncated set of basis functions $\{\psi_{j1},\ldots,\psi_{jd}\}$. Combined with the soft-thresholding step, the update for $f_{j}$ in Algorithm (1) can thus be seen to solve the following problem,

where $\|v\|_{2}^{2}$ denotes $\sum_{i=1}^{n}v_{i}^{2}$ and $R_{j}=Y-\sum_{l\neq j}\Psi_{l}\beta_{l}$ is the residual for $f_{j}$. The sparse backfitting algorithm can then be seen to solve,

where $R_{n}$ denotes the squared error term and $\Omega$ denotes the regularization term, and each $\beta_{j}$ is a $d_{n}$-dimensional vector. Let $S$ denote the true set of variables $\{j:f_{j}\neq 0\}$, with $s=|S|$, and let $S^{c}$ denote its complement. Let $\widehat{S}_{n}=\{j:\widehat{\beta}_{j}\neq 0\}$ denote the estimated set of variables from the minimizer $\widehat{\beta}_{n}$ of (47). For the results in this section, we will treat the covariates as fixed.

Suppose that the following conditions hold on the design matrix $X$ in the orthogonal basis $\psi$:

Assume that the truncation dimension $d_{n}$ satisfies $d_{n}\rightarrow\infty$ and $d_{n}=o(n)$. Furthermore, suppose the following conditions, which relate the regularization parameter $\lambda_{n}$ to the design parameters $n,p$, the number of relevant variables $s$, and the truncation size $d_{n}$:

The proof is given in an appendix. Note that condition (50) implies that

since $\frac{1}{\sqrt{n}}\|A\|_{\infty}\leq\|A\|\leq\sqrt{m}\,\|A\|_{\infty}$ for an $m\times n$ matrix $A$. This relates the condition to previous $\infty$-norm incoherence conditions that have been used for sparsistency in the linear case (Wainwright, 2006).

For $\nu=2$ we take $d_{n}=n^{1/5}$, which achieves the minimax error rate in the one-dimensional case. The theorem under this design setting, with the simplifying assumption that $s=O(1)$, gives the following

Suppose that $s=O(1)$, and $d_{n}=O(n^{1/5})$. Assume the design conditions (48), (49) and (50). Suppose the penalty $\lambda_{n}$ is chosen to satisfy

For example, under these conditions and $\rho_{n}^{*}\asymp 1$, the dimension $p_{n}$ can be taken as large as $e^{n^{3/5}}$; a suitable choice for the regularization parameter in this case would be $\lambda_{n}={Cn^{-\frac{1}{10}+\delta}\log n}$ for some $\delta>0$.

2 Persistence

The previous assumptions are very strong. They can be weakened at the expense of getting weaker results. In particular, in the section we do not assume that the true regression function is additive. We use arguments like those in Juditsky and Nemirovski (2000) and Greenshtein and Ritov (2004) in the context of linear models. In this section we treat $X$ as random and we use triangular array asymptotics, that is, the joint distribution for the data can change with $n$. Let $(X,Y)$ denote a new pair (independent of the observed data) and define the predictive risk when predicting $Y$ with $v(X)$ by

When $v(x)=\sum_{j}\beta_{j}g_{j}(x_{j})$ we also write the risk as $R(\beta,g)$ where $\beta=(\beta_{1},\ldots,\beta_{p})$ and $g=(g_{1},\ldots,g_{p})$. Following Greenshtein and Ritov (2004) we say that an estimator $\widehat{m}_{n}$ is persistent (risk consistent) relative to a class of functions ${\cal M}_{n}$, if

In this section, we assume that the SpAM estimator $\widehat{m}_{n}$ is chosen to minimize

subject to $\left\|\beta\right\|_{1}\leq L_{n}$ and $g_{j}\in{\cal T}_{j}$. We make no assumptions about the design matrix. Let ${\cal M}_{n}\equiv{\cal M}_{n}(L_{n})$ be defined by

and let $m_{n}^{*}=\mathop{\text{arg\,min}}_{v\in{\cal M}_{n}}R(v)$.

Suppose that $p_{n}\leq e^{n^{\xi}}$ for some $\xi<1$. Then,

and hence , if $L_{n}=o(n^{(1-\xi)/4})$ then SpAM is persistent.

Discussion

An additional direction for future work is to develop procedures for automatic bandwidth selection in each dimension. We have used plug-in bandwidths and truncation dimensions $d_{n}$ in our experiments and theory. It is of particular interest to develop procedures that are adaptive to different levels of smoothness in different dimensions.

Finally, we note that while we have considered basic additive models that allow functions of individual variables, it is natural to consider interactions, as in the functional ANOVA model. One challenge is to formulate suitable incoherence conditions on the functions that enable regularization based procedures or greedy algorithms to recover the correct interaction graph. In the parametric setting, one result in this direction is Wainwright et al. (2007).

Proofs

Consider the minimization of the Lagrangian

Using iterated expectations, the above condition can be rewritten as

Rewriting equation (67) in light of (73), we obtain

Using (70), we thus arrive at the soft thresholding update for $f_{j}$:

Our argument closely follows the approach of Wainwright (2006) in the linear case. In particular, we proceed by a “witness” proof technique, to show the existence of a coefficient-subgradient pair $(\beta,g)$ for which ${\text{supp}}(\beta)={\text{supp}}(\beta^{*})$. To do so, we first set $\widehat{\beta}_{S^{c}}={0}$ and $\widehat{g}_{S}=\partial\Omega(\beta^{*})_{S}$, and we then obtain $\widehat{\beta}_{S}$ and $\widehat{g}_{S^{c}}$ from the stationary conditions in (75). By showing that, with high probability,

this demonstrates that with high probability there exists an optimal solution to the optimization problem in (46) that has the same sparsity pattern as the true model.

Setting $\widehat{\beta}_{S^{c}}={0}$ and

for $j\in S$, the stationary condition for $\beta_{S}$ is

Let $V=Y-\Psi_{S}\beta^{*}_{S}-W$ denote the error due to finite truncation of the orthogonal basis, where $W=(\epsilon_{1},\ldots,\epsilon_{n})^{T}$. Then the stationary condition can be written as

assuming that $\frac{1}{n}\Psi_{S}^{T}\Psi_{S}$ is nonsingular. Recalling our definition

in order to ensure that ${\text{supp}}(\beta_{S}^{*})={\text{supp}}(\beta_{S})=\left\{j\,:\,\|\beta_{j}\|_{\infty}\neq 0\right\}$.

We now proceed to bound the quantities above. First note that for $j\in S$,

and thus $\|g_{j}\|\leq\sqrt{C_{\max}}$. Therefore, since

Now, to bound $\left\|{\textstyle{\frac{1}{n}}}\Psi_{S}^{T}V\right\|_{\infty}$, first note that, as we are working over the Sobolev spaces $\mathcal{S}_{j}$ of order two,

for some constant $B^{\prime}>0$. Therefore,

where $D$ denotes a generic constant. Together then, we have that

By Gaussian comparison results (Ledoux and Talagrand, 1991), we have then that

An application of Markov’s inequality then gives that

which converges to zero under the condition that

which is condition (53) in the statement of the theorem.

We now analyze $\widehat{g}_{S^{c}}$. Recall that we have set $\widehat{\beta}_{S^{c}}=\beta^{*}_{S^{c}}={0}$. The stationary condition for $j\in S^{c}$ is thus given by

it suffices to show that $\max_{j\in S^{c}}\|g_{j}\|\leq\sqrt{C_{\min}}$.

From (104), we see that $\widehat{g}_{j}$ is Gaussian, with mean

By our earlier calculations, we have that

which are conditions (50) and (51) in the statement of the theorem. Then

and in particular $\|\mu_{j}\|\leq\sqrt{C_{\min}}$ for sufficiently large $n$. It therefore suffices to show that

with probability approaching one. To show (115), we again appeal to Gaussian comparison results. Define

for $j\in S^{c}$. Then $Z_{j}$ are zero mean Gaussian random variables, and we need to show that

which converges to zero under the condition that

This is condition (52) in the statement of the theorem. $\;\;\scriptstyle\Box$

for some $K>0$. The bracketing integral is defined to be

From Corollary 19.35 of van der Vaart (1998),

Set $Z\equiv(Z_{0},\ldots,Z_{p})=(Y,X_{1},\ldots,X_{p})$ and note that

where we define $g_{0}(z_{0})=z_{0}$ and $\beta_{0}=-1$. Also define

Hence $\widehat{m}_{n}$ is the minimizer of $\widehat{R}(\beta,g)$ subject to the constraint $\sum_{j}\beta_{j}g_{j}(x_{j})\in{\cal M}_{n}(L_{n})$ and $g_{j}\in{\cal T}_{j}$. For all $(\beta,g)$,

Hence, $J_{[\,]}(C,{\mathcal{M}}_{n})=O(\sqrt{\log p_{n}})$ and it follows from (126) and Markov’s inequality that

and the conclusion follows. $\;\;\scriptstyle\Box$

Acknowledgements

This research was supported in part by NSF grant CCF-0625879 and a Siebel Scholarship to PR. A preliminary version of part of this work appears in Ravikumar et al. (2008).

References