Risk Consistency of Cross-Validation with Lasso-Type Procedures

The criterion we focus on for this paper is risk consistency, (alternatively known as persistence). That is, we investigate the difference between the prediction risk of the lasso estimator with tuning parameter estimated by cross-validation and the risk of the best linear oracle predictor (with oracle tuning parameter). Risk consistency of lasso has previously been studied by Greenshtein and Ritov (2004); Bunea, Tsybakov, and Wegkamp (2007); van de Geer (2008) and Bartlett, Mendelson, and Neeman (2012). Their results, in contrast to ours, assume that the tuning parameter is selected in an oracle fashion.

We present two main results which make progressively stronger assumptions on the data generating process and use both forms of the lasso estimator. The first imposes strong conditions on the design matrix, assumes the linear model is true, and that this linear model is sparse. The second is more general and allows the true model to be neither linear nor correctly specified. For this reason, our focus is on risk consistency rather than estimation of a “true” parameter or correct identification of a “true” sparsity pattern. Additionally, well-known results of Shao (1993) imply that cross-validation leads to inconsistent model selection in general, suggesting that results for sparse recovery may not exist. Although prediction is an important goal, one is often interested in variable selection for more interpretable models or follow-up experiments. In light of the negative results in Shao (1993), we are unable to offer theoretical results which promise consistent model selection by cross validation, but simulations in Section 4 suggest that it performs well nevertheless. Both the estimation and sparse recovery settings are frequently studied in the literature assuming the tuning parameter is the oracle and that the data generating model is linear (e.g. Bunea, Tsybakov, and Wegkamp, 2007; Candès and Plan, 2009; Donoho, Elad, and Temlyakov, 2006; Leng, Lin, and Wahba, 2006; Meinshausen and Bühlmann, 2006; Meinshausen and Yu, 2009).

In the first case, when the truth is linear, we examine the Lagrangian formulation in Equation (1). In this scenario, we prove convergence rates which differ only by a $\log$ factor relative to those achievable with the oracle tuning parameter (e.g. Negahban, Ravikumar, Wainwright, and Yu, 2012; Bühlmann and van de Geer, 2011; Bunea, Tsybakov, and Wegkamp, 2007). That is, for an $s^{*}$-sparse linear model with restricted isometry conditions on the covariance of the design, the risk of the cross-validated estimator approaches the oracle risk at a rate of $O(s^{*}\log(p)\log(n)/n)$. Under more general conditions, we follow the approach of Greenshtein and Ritov (2004) and examine the constrained optimization form in Equation (2). Using our methods, we require that the set of candidate predictors, $\mathcal{B}_{t_{n}}$, satisfies $t_{n}=o\left(n^{1/4}/(m_{n}(\log p)^{1/4+1/(2q))}\right)$ where $m_{n}$ is a sequence which approaches infinity and $q$ characterizes the tail behavior of the data. This is essentially as quickly as one could hope relative to Greenshtein and Ritov (2004) under our more general assumptions on the design matrix. We note however that, using empirical process techniques, Bartlett, Mendelson, and Neeman (2012) have been able to improve the rate shown in Greenshtein and Ritov (2004) to $t_{n}=o\left(n^{1/2}/(\log^{3/2}n\log^{3/2}(np))\right)$ for sub-Gaussian design and oracle tuning parameter.

There are several proposed data-dependent techniques for choosing $t$ or $\lambda$. Kato (2009) and Tibshirani and Taylor (2012) investigate estimating the “degrees of freedom” of a lasso solution. With an unbiased estimator of the degrees of freedom, the tuning parameter can be selected by minimizing the empirical risk penalized by a function of this estimator. Note that this approach requires an estimate of the variance, which is nontrivial when $p>n$ (Giraud, Huet, and Verzelen, 2012; Homrighausen and McDonald, 2016). Another risk estimator is the adapted Bayesian information criterion proposed by Wang and Leng (2007) which uses a plug-in estimator based on the second-order Taylor’s expansion of the risk. Arlot and Massart (2009) and Saumard (2011) consider “slope heuristics” as a method for penalty selection with Gaussian noise, paying particular attention to the regressogram estimator in the first case and piecewise polynomials with $p$ fixed in the second. Finally, Sun and Zhang (2012) present an algorithm to jointly estimate the regression coefficients and the noise level. This results in a data-driven value for the tuning parameter which possesses oracle properties under some regularity conditions.

However, many authors (e.g. Efron et al., 2004; Friedman et al., 2010; Greenshtein and Ritov, 2004; Tibshirani, 1996, 2011; Zou et al., 2007 and as discussed by van de Geer and Bühlmann, 2011, Section 2.4.1) recommend selecting $t$ or $\lambda$ in the lasso problem by minimizing a $K$-fold cross-validation estimator of the risk (see Section 2 for the precise definition). Cross-validation is generally well-studied in a number of contexts, especially model selection and risk estimation. In the context of model selection, Arlot and Celisse (2010) give a detailed survey of the literature emphasizing the relationship between the sizes of the validation set and the training set as well as discussing the positive bias of cross-validation as a risk estimator.

Some results supporting the use of cross-validation for statistical methods other than lasso are known. For instance, Stone (1974, 1977) outlines various conditions under which cross-validated methods can result in good predictions. More recently, Dudoit and van der Laan (2005) find finite sample bounds for various cross-validation procedures. These results do not address the lasso nor parameter spaces with increasing dimensions, and furthermore, apply to choosing the best estimator from a finite collection of candidate estimators. Lecué and Mitchell (2012) provide oracle inequalities related to using cross-validation with lasso, however, it treats the problem as aggregating a dictionary of lasso estimators with different tuning parameters, and the results are stated for fixed $p$ rather than the high-dimensional setting investigated here. Most recently, Flynn, Hurvich, and Simonoff (2013) investigate numerous methods for tuning parameter selection in penalized regression, but the theoretical results hold only when $p/n\rightarrow 0$ and not for cross-validation. In particular, the authors state “to our knowledge the asymptotic properties of [$K$]-fold CV have not been fully established in the context of penalized regression” (Flynn, Hurvich, and Simonoff, 2013, p. 1033).

Rather than cross-validation, one may use information criteria such as AIC (Akaike, 1974) or BIC (Schwarz, 1978). These techniques are frequently advocated in the literature (for example Bühlmann and van de Geer, 2011; Wang et al., 2007; Tibshirani, 1996; Fan and Li, 2001), but the classical asymptotic arguments underlying these methods apply only for $p$ fixed and rely on maximum likelihood estimates (or Bayesian posteriors) for all parameters including the noise. The theory in the high-dimensional setting supporting these methods is less complete. Recent work has developed new information criteria with supporting asymptotic results uf $\textrm{rank}(\mathbb{X})=p$ but is still allowed to increase. For example, the criterion developed by Wang et al. (2009) selects the correct model asymptotically even if $p\rightarrow\infty$ as long as $p/n\rightarrow 0$. If $p$ is allowed to increase more quickly than $n$, theoretical results assume $\sigma^{2}$ is known to get around the difficult task of high-dimensional variance estimation (Chen and Chen, 2012; Zhang and Shen, 2010; Kim et al., 2012; Fan and Tang, 2013).

In Section 2, we outline the mathematical setup for the lasso prediction problem and discuss some empirical concerns. Section 3 contains the main result and associated conditions. Section 4 compares different choices of $K$ for cross-validation via simulation, while Section 5 summarizes our contribution and presents some avenues for further research.

Suppose we observe data $Z_{i}^{\top}=(Y_{i},X_{i}^{\top})$ consisting of predictor variables, $X_{i}\in\mathbb{R}^{p_{n}}$, and response variables, $Y_{i}\in\mathbb{R}$, where $Z_{i}\sim\mu_{n}$ are independent and identically distributed for $i=1,2,\ldots,n$ and the distribution $\mu_{n}$ is in some class $\mathcal{F}$ to be specified. Here, we use the notation $p_{n}$ to allow the number of predictor variables to change with $n$. Similarly, we index the distribution $\mu_{n}$ to emphasize its dependence on $n$. For simplicity, in what follows, we omit the subscript $n$ when there is little risk of confusion.

We consider the problem of estimating a linear functional $f(\mathcal{X})=\mathcal{X}^{\top}\beta$ for predicting $\mathcal{Y}$, when $\mathcal{Z}^{\top}=(\mathcal{Y},\mathcal{X}^{\top})\sim\mu_{n}$ is a new, independent random variable from the same distribution as the data and $\beta=(\beta_{1},\ldots,\beta_{p})^{\top}$. For now, we do not assume a linear model, only the usual regression framework where $\mathcal{Y}=f^{*}(\mathcal{X})+\epsilon,$ with $\epsilon$ and $\mathcal{X}$ independent and $f^{*}$ is some unknown function. We will use zero-based indexing for $\mathcal{Z}$ so that $\mathcal{Z}_{0}=\mathcal{Y}$. To measure performance, we use the $L^{2}$-risk of the predictor $\beta$:

Note that this is a conditional expectation: for any estimator $\widehat{\beta}$,

and the expectation is taken only over the new data $\mathcal{Z}$ and not over any observables which may be used to choose $\widehat{\beta}$.

Using the $n$ independent observations $Z_{1},\ldots,Z_{n}$, we can form the response vector $Y:=(Y_{i})_{i=1}^{n}$ and the design matrix $\mathbb{X}:=[X_{1},\ldots,X_{n}]^{\top}$. Then, given a vector $\beta$, we write the squared-error empirical risk function as

Analogously to Equation (5), we write the $K$-fold cross-validation estimator of the risk with respect to the tuning parameter $t$, which we abbreviate to CV-risk or just CV, as

Here, $V_{n}=\{v_{1},\ldots,v_{K}\}$ is a set of validation sets, $\widehat{\beta}_{t}^{(v)}$ is the estimator in Equation (2) with the observations in the validation set $v$ removed, and $|v|$ indicates the cardinality of $v$. Notice in particular that the cross-validation estimator of the risk is a function of $t$ rather than a single predictor $\beta$—it uses $K$ different predictors at a fixed $t$, averaging over their performance on the respective held-out sets. Over the course of the paper, we will freely exchange $\lambda$ for $t$ in this definition.

Lastly, we define the CV-risk minimizing choice of tuning parameter to be

In our setting, we will take $T$ (or $\Lambda$) as an interval subset of the nonnegative real numbers which needs to be defined by the data-analyst. The choice of $T$ is an important part of the performance of $\widehat{\beta}_{\widehat{t}}$ and requires some explanation. First, we provide some insight into the computational load of using CV-risk to find $\widehat{t}$.

In practice, CV-risk is known to be time consuming and somewhat unstable due to the randomness associated with forming $V_{n}$. For a fixed $v\in V_{n}$, suppose $\widehat{\beta}_{\lambda}^{(v)}$ is found for the entire lasso path via the lars (Efron et al., 2004) algorithm, which can be computed in the same computational complexity as a least squares fit. To fix ideas, suppose $n>p$, which means lars has computational complexity $O(np^{2}+p^{3})$. Hence, as the feature matrix $\mathbb{X}$ with $|v|$ rows removed has approximately $n(K-1)/K$ rows, $\widehat{\beta}_{\lambda}^{(v)}$ can be computed for all $\lambda$ in $O((n(K-1)/K)p^{2}+p^{3})$ time. Repeating this $K$ times means the computational complexity for forming $\widehat{R}_{V_{n}}\left(\lambda\right)$ over all $\lambda$ is $O(n(K-1)p^{2}+p^{3})$. If $K$ is a fixed fraction of $n$, CV-risk has computational complexity of order $(np)^{2}$, which is a factor of $n$ slower than a single lasso fit.

Though more expensive on a single processor, CV-risk is readily parallelizable over the $K$ folds and therefore (ignoring communication costs between processors) CV-risk could actually be faster to compute than $\widehat{R}$ (and hence $\widehat{\beta}_{\lambda}$) as

The data analyst must be able to solve the optimization problem in Equation (7). For $\Lambda$, we must choose a lower bound: $\Lambda=[\lambda_{n},\infty)$. This implies we must choose $\lambda_{n}$ as a function of the data. While it is tempting to allow $\lambda_{n}=0$, this results in numerical and practical implementation issues if rank$(\mathbb{X})<p$ and is unnecessary as the theory will show. However, the lower bound will have a nontrivial impact on the quality of the recovery, as choosing a value too large may eliminate the best solutions. We will see in the next section that under some assumptions on the data generating process, we can safely choose a particular $\lambda_{n}>0$ that allows order $\log n$ coefficients to be selected without compromising the quality of the estimator.

In the case of $T$, an upper bound must be selected for any grid-search optimization procedure. As we will impose much weaker conditions on the data generating process, choosing such an upper bound is more complicated. Note that, by Equation (2), $\widehat{\beta}_{t}$ must be in the $\ell_{1}$-ball with radius $t$. This constraint is only binding (Osborne, Presnell, and Turlach, 2000) if

where $\widehat{\beta}_{\infty}=\widehat{\beta}(\mathbb{R}^{p})$ is a least-squares solution and $\mathcal{K}:=\{a:\mathbb{X}a=0\}$. Observe that $\mathcal{K}=\{0\}$ if $n\geq p$ and otherwise $\mathcal{K}$ has dimension $p-n$, which would imply $\widehat{\beta}_{\infty}$ is not unique. Both of these statements assume that the columns of $\mathbb{X}$ contain a linearly independent set of size $\min\{n,p\}$. See Tibshirani (2013) for the more general situation. In either case, if $t\geq t_{0}$, then $\widehat{\beta}_{t}$ is equal to a least-squares solution. Based on this argument and the fact that $0\in\mathcal{K}$, it is tempting to define the upper bound to be $t_{\max}:=||\widehat{\beta}_{\infty}||_{1}$, where $\widehat{\beta}_{\infty}=(\mathbb{X}^{\top}\mathbb{X})^{\dagger}\mathbb{X}^{\top}Y$ is the least-squares solution when $(\cdot)^{\dagger}$ is given by the Moore-Penrose inverse. However, this upper bound has at least two problems. First, although theoretically well-defined, as with setting $\lambda_{n}=0$, it suffers from numerical and practical implementation issues if rank$(\mathbb{X})<p$. The second problem is that it grows much too fast, at least order $\sqrt{n}$, therefore potentially including solutions which will have low bias but very high variance.

Instead, we define an upper bound based on a rudimentary variance estimator $t_{\max}:=\left|\left|Y\right|\right|_{2}^{2}/a_{n},$ where $(a_{n})$ is an increasing sequence of constants. Hence, we consider the optimization interval to be $T=[0,t_{\max}]$. We defer fixing a particular sequence $(a_{n})$ until the next section. As our main results demonstrate, this choice of $T$ is enough to ensure that a risk consistent sequence of tuning parameters can be selected.

If we are willing to impose strong conditions on $\mu_{n}$, as in  Bunea, Tsybakov, and Wegkamp (2007) and Meinshausen (2007), then we can obtain cross-validated lasso estimators which achieve nearly oracle rates.

If we assume that $\mathbb{E}[\mathcal{Y}\ |\ \mathcal{X}]=f^{*}(\mathcal{X})=\mathcal{X}^{\top}\beta^{*}$, then we can write the risk of an estimator $\widehat{\beta}_{\widehat{\lambda}}$ as

where $R(\beta^{*})=\mathbb{E}[(\mathcal{Y}-\mathcal{X}^{\top}\beta^{*})^{2}]=\sigma^{2}$. We write the excess risk as $\mathcal{E}(\widehat{\lambda}):=R(\widehat{\beta}_{\widehat{\lambda}})-\sigma^{2}$ and prove a convergence rate for $\mathcal{E}(\widehat{\lambda})$. In this case, targeting the excess risk is the same as estimating the conditional expectation of $\mathcal{Y}$, but if $f^{*}(\mathcal{X})$ is not linear (as in Section 3.2), the excess risk remains a meaningful way of assessing prediction behavior.

These conditions and the result warrant a few comments. First, note that condition M4 implies that $\Sigma-(1-\nu)\textrm{diag}(\Sigma)$ is positive semi-definite. Second, as $s^{*}$ is fixed, M6 and M7 ensure $\lambda_{\min}\rightarrow 0$ so that $\Lambda$ will eventually allow models with $s^{*}$ covariates. Thus, the procedure is asymptotically consistent for model selection (Meinshausen, 2007). Finally, note that $\mathcal{E}(\widehat{\lambda})$ is random due to the term $R(\widehat{\beta}_{\widehat{\lambda}})$. Here, we emphasize that $R(\widehat{\beta}_{\widehat{\lambda}})$ is a function of the data: the conditional expectation of a new test random variable $\mathcal{Z}$ given the observed data used to choose both $\widehat{\lambda}$ and $\widehat{\beta}_{\lambda}$ as in Equation (4).

While the conditions of Theorem 3 are strong, they are typical of those used to prove persistence of the lasso estimator with oracle tuning parameter (for the case of fixed design, see Negahban et al., 2012). For instance, Bunea et al. (2007) prove an oracle rate for the lasso of $O(s^{*}\log{p}/n)$ with $\lambda_{\min}\propto\sigma\sqrt{\log{p}/n}$. Under similar conditions, our result with cross-validated tuning parameter requires a larger $\lambda_{\min}$ (resulting in smaller models) and a slower convergence rate to the oracle by a factor $\log{n}$. A reasonable choice of $\Lambda$ suggested by Theorem 3 is $\Lambda=[\lambda_{\min},\ \infty)=[(\log{p}\log{n}/n)^{1/2},\ \infty)$.

In order to derive results under more general conditions, we move to the linear oracle estimation framework. For any $t$, define the oracle estimator with respect to $t$ as

Suppose $\widehat{t}$ is an estimator, such as by cross-validation. Then we can decompose the risk of an estimator $\widehat{\beta}_{\widehat{t}}$ as follows:

where $R_{*}$ is the risk of the mean function $f^{*}$. Because the data generating process is not necessarily linear, we study the performance of an estimator $\widehat{\beta}_{\widehat{t}}$ via the excess risk of $\widehat{\beta}_{\widehat{t}}$ relative to $\beta_{t}$, which we define as

Here, $\mathcal{E}(\widehat{t},t)$ depends on the cross-validated tuning parameter $\widehat{t}$ as well as the oracle set through $t$. Focusing on Equation (9), allows for meaningful theory when the approximation error does not necessarily converge to zero as $n$ grows. This is particularly important here, as we do not assume that the conditional expectation of $\mathcal{Y}$ given $\mathcal{X}$ is linear. As $t\rightarrow\infty$, the approximation error decreases and hence we desire to take $t=t_{n}$ for some increasing sequence $(t_{n})$. As discussed in the introduction, Greenshtein and Ritov (2004) show that if $t_{n}=o((n/\log p)^{1/4})$, then $\mathcal{E}(t_{n},t_{n})$ converges in probability to zero. Bartlett, Mendelson, and Neeman (2012) increase the size of this search set allowing $t_{n}=o(n^{1/2}/(\log^{3/2}n\log^{3/2}(np)))$ while still having $\mathcal{E}(t_{n},t_{n})\xrightarrow{P}0.$

Define $b_{n}=\min\{n-c_{n},c_{n}\}$. It follows that $(n^{-1/2}+c_{n}^{-1/2}+(n-c_{n})^{-1/2})\leq 3b_{n}^{-1/2}.$ We state the following corollary to Theorem 4, which outlines the conditions under which we can do asymptotically at least as well at predicting a new observation as the oracle linear model. That is, when the right hand side of Equation (10) goes to zero.

The parameter $b_{n}$ controls the minimum size of the validation versus training sets that comprise cross-validation. It must be true that $b_{n}$ is strictly less than $n$. To get the best guarantee, $b_{n}$ should increase as fast as possible. Hence, our results advocate a cross-validation scheme where the validation and training sets are asymptotically balanced, i.e. $c_{n}\asymp n-c_{n}$. This should be compared with the results in Shao (1993), which state that for model selection consistency, one should have $c_{n}/n\rightarrow 1$. However, Shao (1993) presents results for model selection while we focus on prediction error. Similarly, Dudoit and van der Laan (2005) provide oracle inequalities for cross-validation and also advocate for the validation set to grow as fast as possible. Note that for $K$-fold cross-validation, $c_{n}=\lfloor n/K\rfloor$ so that $b_{n}=O(n)$.

Additionally, the rates $a_{n}$ and $m_{n}$ deserve comment. It is instructive to compare this choice of $t_{\textrm{max}}$ with $\left|\left|Y\right|\right|_{2}^{2}/n$, a standard estimate of the noise variance in high dimensions (e.g. van de Geer and Bühlmann, 2011, p. 104). If $a_{n}=n$, then $\left|\left|Y\right|\right|_{2}^{2}/a_{n}$ is an overestimate of the variance due to the presence of the regression function $f^{*}$. Our results state that we must choose $a_{n}$ to increase slower than $n$, thereby increasing this overestimation and enlarging the potential search set $T$. While $a_{n}$ depends on several parameters, $b_{n}$, $n$, and $p$ are known to the analyst. Also, the choice of $q$ depends on how much approximation error the analyst is willing to make. The $m_{n}$ term is required to slow the growth of $t_{n}$ just slightly. While this shrinks the size of the set $\mathcal{B}_{t_{n}}$ relative to that used by Greenshtein and Ritov (2004), potentially eliminating some better solutions, effectively $m_{n}$ quantifies their requirement $t_{n}=o((n/\log p)^{1/4})$, making explicit the need for $t_{n}$ to grow more slowly than $(n/\log p)^{1/4}$. As such, if we set $b_{n}\asymp n$ and set $q=\infty$, we reaquire the rate shown by Greenshtein and Ritov (2004), where Section 3.2 implies that both $\mathbb{P}_{\mu_{n}}\left(\mathcal{E}(\widehat{t},t_{n})>\delta\right)\rightarrow 0$ and $t_{n}=o((n/\log p)^{1/4})$.

While the entire proof is contained in the supplementary material, we provide a brief sketch here.

We note here that our results generalize to other $M$-estimators which use an $\ell_{1}$-constraint. In particular, relative to the set of coefficients $\beta\in\mathcal{B}_{t_{n}}$ with $t_{n}=o\left(b_{n}^{1/4}/\left(m_{n}(\log p)^{1/4+1/(2q)}\right)\right)$, an empirical estimator with cross-validated tuning parameter has prediction risk which converges to the prediction risk of the oracle.