On Cross-validation for Sparse Reduced Rank Regression

TL;DR

This paper investigates cross-validation strategies for sparse reduced rank regression in high-dimensional settings, proposing new criteria and calibration methods to improve model selection and estimation accuracy.

Contribution

It introduces a novel approach to cross-validation that focuses on projection-selection patterns and develops scale-free information criteria for better regularization parameter tuning.

Findings

Proposes a new class of scale-free information criteria.

Introduces a scale-free calibration method for minimax optimal error rate.

Experimental results demonstrate the effectiveness of the proposed methods.

Abstract

In high-dimensional data analysis, regularization methods pursuing sparsity and/or low rank have received a lot of attention recently. To provide a proper amount of shrinkage, it is typical to use a grid search and a model comparison criterion to find the optimal regularization parameters. However, we show that fixing the parameters across all folds may result in an inconsistency issue, and it is more appropriate to cross-validate projection-selection patterns to obtain the best coefficient estimate. Our in-sample error studies in jointly sparse and rank-deficient models lead to a new class of information criteria with four scale-free forms to bypass the estimation of the noise level. By use of an identity, we propose a novel scale-free calibration to help cross-validation achieve the minimax optimal error rate non-asymptotically. Experiments support the efficacy of the proposed methods.

Tables

Table 1. Table 1: Performance comparisons between AIC, BIC, EBIC, PIC, 2 2 2 -fold CV, 10 10 10 -fold CV, 5 5 5 -fold CV and 5 5 5 -fold SCV in the n > p 𝑛 𝑝 n>p experiment with smaller signal strength ( b = 0.1 𝑏 0.1 b=0.1 ). MSEs are scaled for ease of comparison and M- and FA-rates are in percentages.

	$\rho =0.1$					$\rho =0.5$
	MSE	$\widehat{𝑱}$	$\widehat{𝒓}$	M	FA	MSE	$\widehat{𝑱}$	$\widehat{𝒓}$	M	FA
AIC	43	37	6	4	29	44	38	6	7	35
BIC	45	20	3	33	0	42	18	3	41	1
EBIC	104	8	1	74	0	103	7	1	75	0
\hdashlinePIC	26	29	4	11	7	28	27	4	18	9
2-CV	99	50	10	1	68	99	50	10	3	69
10-CV	99	50	10	1	68	99	50	10	3	69
5-CV	96	50	10	1	68	96	50	10	3	69
\hdashline5-SCV	28	29	5	10	9	30	27	5	19	9

Table 2. Table 2: Performance comparisons between AIC, BIC, EBIC, PIC, 2 2 2 -fold CV, 10 10 10 -fold CV, 5 5 5 -fold CV and 5 5 5 -fold SCV in the n > p 𝑛 𝑝 n>p experiment with larger signal strength ( b = 0.5 𝑏 0.5 b=0.5 ). MSEs are scaled for ease of comparison and M- and FA-rates are in percentages.

	$\rho =0.1$					$\rho =0.5$
	MSE	$\widehat{𝑱}$	$\widehat{𝒓}$	M	FA	MSE	$\widehat{𝑱}$	$\widehat{𝒓}$	M	FA
AIC	13	34	6	0	17	12	33	6	0	17
BIC	7	30	5	0	0	7	30	5	0	0
EBIC	7	30	5	0	0	7	30	5	0	0
\hdashlinePIC	7	30	5	0	1	7	30	5	0	0
2-CV	40	50	10	0	67	40	50	10	0	67
10-CV	40	50	10	0	67	40	50	10	0	67
5-CV	40	50	10	0	66	40	50	10	0	66
\hdashline5-SCV	7	30	5	0	2	7	30	5	0	1

Table 3. Table 3: Performance comparisons between AIC, BIC, EBIC, PIC, 2 2 2 -fold CV, 10 10 10 -fold CV, 5 5 5 -fold CV and 5 5 5 -fold SCV in the p > n 𝑝 𝑛 p>n experiment with smaller signal strength ( b = 0.2 𝑏 0.2 b=0.2 ). MSEs are scaled for ease of comparison and M- and FA-rates are in percentages.

	$\rho =0.1$					$\rho =0.5$
	MSE	$\widehat{𝑱}$	$\widehat{𝒓}$	M	FA	MSE	$\widehat{𝑱}$	$\widehat{𝒓}$	M	FA
AIC	93	24	4	44	17	79	24	4	36	16
BIC	42	12	1	59	7	31	10	2	56	5
EBIC	43	7	1	71	3	36	4	1	74	1
\hdashlinePIC	40	10	2	62	5	28	9	2	59	3
2-CV	122	25	4	43	19	97	25	4	35	18
10-CV	130	25	4	43	19	99	25	4	35	18
5-CV	128	25	4	43	19	99	25	4	35	18
\hdashline5-SCV	41	10	2	62	5	29	8	2	61	3

Table 4. Table 4: Performance comparisons between AIC, BIC, EBIC, PIC, 2 2 2 -fold CV, 10 10 10 -fold CV, 5 5 5 -fold CV and 5 5 5 -fold SCV in the p > n 𝑝 𝑛 p>n experiment with larger signal strength ( b = 1 𝑏 1 b=1 ). MSEs are scaled for ease of comparison and M- and FA-rates are in percentages.

	$\rho =0.1$					$\rho =0.5$
	MSE	$\widehat{𝑱}$	$\widehat{𝒓}$	M	FA	MSE	$\widehat{𝑱}$	$\widehat{𝒓}$	M	FA
AIC	47	25	2	40	18	31	25	2	32	16
BIC	43	24	2	40	17	30	24	2	33	15
EBIC	34	21	2	43	12	26	19	2	38	11
\hdashlinePIC	41	23	2	40	15	28	21	2	35	13
2-CV	47	25	2	40	18	34	25	3	31	17
10-CV	48	25	3	40	18	34	25	3	31	17
5-CV	50	25	2	39	18	35	25	4	31	17
\hdashline5-SCV	34	11	2	58	5	22	11	2	53	4

Table 5. Table 5: Prediction error comparison between PIC and SCV in the p > n 𝑝 𝑛 p>n experiment with different signal strengths. Reported numbers are MSE’s with standard errors in parentheses (scaled for ease of comparison).

	$\rho =.1$			$\rho =.5$
	$b=2$	$b=3$	$b=4$	$b=2$	$b=3$	$b=4$
PIC	42 (22)	42 (23)	43 (23)	28 (19)	29 (19)	29 (19)
SCV	34 (16)	33 (17)	33 (17)	21 (10)	21 (10)	21 (10)

Equations

$$\mathbf{X}\mathbf{B}$$

$$=\{\mathbf{X}\mathbf{S}(\mathbf{B})\mathbf{U}(\mathbf{B})\}\cdot \{\mathbf{U}(\mathbf{B}{)}^T\mathbf{S}(\mathbf{B}{)}^T\mathbf{B}\}.$$

$$\underset{\mathbf{B}}{\min }\Vert \mathbf{Y}-\mathbf{X}\mathbf{B}{\Vert }_F^2\text{mbox}s.t.r(\mathbf{B})\le r,\mathbf{B}[{\mathcal{J}}^c,]=0.$$

$$\mathbb{E}[\Vert {\mathbf{Y}}^{\prime }-\mathbf{X}\hat{\mathbf{B}}{\Vert }_F^2]$$

$$=$$

$$\Vert \mathbf{Y}-\mathbf{X}\hat{\mathbf{B}}{\Vert }_F^2+2{\sigma }^2\bar{r}m$$

$$\begin{array}{r}S(J,r)=\{\mathbf{B}\in {\mathbb{R}}^{p\times m}:J(\mathbf{B})\le J,r(\mathbf{B})\le r\}.\end{array}$$

$$P_o(\mathbf{B})=[q\wedge J(\mathbf{B})+m-r(\mathbf{B})]r(\mathbf{B})+J(\mathbf{B})\log (ep/J(\mathbf{B})).$$

$$\underset{\hat{\mathbf{B}}}{\inf }\underset{{\mathbf{B}}^{\ast }\in S(J,r)}{\sup }\mathbb{E}[l(\Vert \mathbf{X}{\mathbf{B}}^{\ast }-\mathbf{X}\hat{\mathbf{B}}{\Vert }_F^2/(c{\sigma }^2{P}_o(J,r)))]\ge c^{\prime }>0,$$

$$\hat{\mathbf{B}}={\mathrm{argmin}}_{{\mathbf{B}}_l}\frac{1}{2}\Vert \mathbf{Y}-\mathbf{X}{\mathbf{B}}_l{\Vert }_F^2+A{\sigma }^2{P}_o({\mathbf{B}}_l),$$

$$\mathbb{E}[\Vert \mathbf{X}\hat{\mathbf{B}}-\mathbf{X}{\mathbf{B}}^{\ast }{\Vert }_F^2]\lesssim \underset{l}{\inf }\mathbb{E}[\Vert \mathbf{X}{\mathbf{B}}_l-\mathbf{X}{\mathbf{B}}^{\ast }{\Vert }_F^2+{\sigma }^2{P}_o({\mathbf{B}}_l)].$$

$$\Vert \mathbf{Y}-\mathbf{X}\mathbf{B}{\Vert }_F^2/[1-\delta (\mathbf{B})],$$

$$\Vert \mathbf{Y}-\mathbf{X}\mathbf{B}{\Vert }_F^2/[1-\delta (\mathbf{B}){]}^2,$$

$$\log \{\Vert \mathbf{Y}-\mathbf{X}\mathbf{B}{\Vert }_F^2\}+\delta (\mathbf{B}),$$

$$\Vert \mathbf{Y}-\mathbf{X}\mathbf{B}{\Vert }_F^2+\delta (\mathbf{B})\Vert \mathbf{Y}-\mathbf{X}\mathbf{B}{\Vert }_F^2.$$

$$\text{mbox}Trn-Err\triangleq \Vert \mathbf{Y}-\mathbf{X}\hat{\mathbf{B}}{\Vert }_F^2,$$

$$\text{mbox}CV-Err\triangleq \sum _{k=1}^K\Vert {\mathbf{Y}}^k-{\mathbf{X}}^k{\hat{\mathbf{B}}}^{-k}{\Vert }_F^2,$$

$$\hat{\mathbf{B}}=\mathbf{S}\mathbf{U}({\mathbf{X}}_{\mathbf{S}\mathbf{U}}^T{\mathbf{X}}_{\mathbf{S}\mathbf{U}}{)}^{-1}{\mathbf{X}}_{\mathbf{S}\mathbf{U}}^T\mathbf{Y},{\hat{\mathbf{B}}}^{-k}=\mathbf{S}\mathbf{U}(({\mathbf{X}}_{\mathbf{S}\mathbf{U}}^{-k}{)}^T{\mathbf{X}}_{\mathbf{S}\mathbf{U}}^{-k}{)}^{-1}({\mathbf{X}}_{\mathbf{S}\mathbf{U}}^{-k}{)}^T{\mathbf{Y}}^{-k}$$

$$\mathbb{E}[\text{mbox}CV-Err]$$

$$\displaystyle D=\$$

$$\begin{array}{rl}U=& \mathbb{E}[\Vert {\mathbf{X}}_{\mathbf{S}\mathbf{U}}({\bar{\mathbf{B}}}_{\mathbf{S}\mathbf{U},{\mathbf{S}}^{\ast }{\mathbf{U}}^{\ast }}-\mathbb{E}{\mathbf{B}}_{\mathbf{S}\mathbf{U},{\mathbf{S}}^{\ast }{\mathbf{U}}^{\ast }}){\Vert }_F^2]\\ +& \mathbb{E}[\Vert {\mathbf{\Sigma }}_{\mathbf{S}\mathbf{U}}^{1/2}({\mathbf{B}}_{\mathbf{S}\mathbf{U},{\mathbf{S}}^{\ast }{\mathbf{U}}^{\ast }}-\mathbb{E}{\mathbf{B}}_{\mathbf{S}\mathbf{U},{\mathbf{S}}^{\ast }{\mathbf{U}}^{\ast }}){\Vert }_F^2].\end{array}$$

$$\mathbb{E}[\text{mbox}CV-Err]=$$

$$\mathbb{E}[\text{mbox}CV-Err]$$

$$+\mathbb{E}[\Vert {\mathbf{X}}_{\mathcal{J}}({\bar{\mathbf{\beta }}}_{\mathcal{J},{\mathcal{J}}^{\ast }}-\mathbb{E}{\mathbf{\beta }}_{\mathcal{J},{\mathcal{J}}^{\ast }}){\Vert }_2^2]+\mathbb{E}[\Vert {\mathbf{\Sigma }}_{\mathcal{J}}^{1/2}({\mathbf{\beta }}_{\mathcal{J},{\mathcal{J}}^{\ast }}-\mathbb{E}{\mathbf{\beta }}_{\mathcal{J},{\mathcal{J}}^{\ast }}){\Vert }_2^2],$$

$$\text{mbox}Trn-err+Am\bar{r}{\sigma }^2.$$

$$\text{mbox}SCV-Err:=\text{mbox}CV-Err+{\alpha }_1(\text{mbox}Trn-Err/(mn))\text{mbox}R+{\alpha }_2(\text{mbox}Trn-Err/(mn))\text{mbox}IF,$$

$$\text{mbox}CV-Err/(1-{\alpha }_1\text{mbox}R/(mn)-{\alpha }_2\text{mbox}IF/(mn))$$

$${\mathcal{B}}^1(r)=\{\mathbf{B}=$$

$$\text{mbox}and{b}_{jk}=0\text{mbox}otherwise\}.$$

$$\log |{\mathcal{B}}^{10}(r)|\ge c_{1}r(J+m-r),$$

$$\rho ({\mathbf{B}}_1,{\mathbf{B}}_2)\ge c_{2}r(J+m-r),\forall {\mathbf{B}}_1,{\mathbf{B}}_2\in {\mathcal{B}}^{10},{\mathbf{B}}_1\ne {\mathbf{B}}_2$$

$$\Vert \mathbf{X}{\mathbf{B}}_1-\mathbf{X}{\mathbf{B}}_2{\Vert }_F^2\ge c_2\underline{\kappa }{\gamma }^2{R}^2(J+m-r)r$$

$$\mathcal{K}(P_0,P_{\mathbf{B}})\le \frac{1}{2{\sigma }^2}\overline{\kappa }{\gamma }^2{R}^2\rho (0,\mathbf{B})\le \frac{{\gamma }^2}{{\sigma }^2}{R}^2(J+m-r)r$$

Full text

1 Background

Modern statistical techniques heavily rely on shrinkage estimation and thus parameter tuning becomes an important task. This work assumes a multivariate regression setup, ${\boldsymbol{Y}}={\boldsymbol{X}}{\boldsymbol{B}}^{*}+{\boldsymbol{E}}$, which is applicable to a variety of real-world applications in machine learning, economics and genetics studies (Stock and Watson,, 2002; Hastie et al.,, 2009; Vounou et al.,, 2012). Here, ${\boldsymbol{Y}}\in\mathbb{R}^{n\times m}$ denotes the response matrix, ${\boldsymbol{X}}\in\mathbb{R}^{n\times p}$ is the design matrix containing $p$ predictors or features, and the noise ${\boldsymbol{E}}$ is assumed to be Gaussian (or sub-Gaussian). The $j$th row of ${\boldsymbol{B}}^{*}$ contains the coefficients associated with the $j$th predictor, and the $k$th column of ${\boldsymbol{B}}^{*}$ corresponds to the coefficients associated with the $k$th response in ${\boldsymbol{Y}}$. Therefore, when a number of features are irrelevant to ${\boldsymbol{Y}}$, ${\boldsymbol{B}}^{*}$ has row sparsity. In the case of a single response, the problem becomes standard variable selection. But variable selection alone may not be sufficient when $m>1$. Suppose that the first $J^{*}$ rows in ${\boldsymbol{B}}^{*}$ are nonzero, i.e., ${\boldsymbol{B}}^{*}=[{\boldsymbol{\beta}}_{1}^{*}\ \cdots\ {\boldsymbol{\beta}}_{J^{*}}^{*}\ \boldsymbol{0}\ \cdots\ \boldsymbol{0}]^{T}$. The number of unknowns, $m\times J^{*}$, can still be large and exceed the total number of observations. In matrix estimation, adding a low-rank constraint is popular and effective (Candès and Plan,, 2011; Rohde and Tsybakov,, 2011). Specifically, suppose further that ${\boldsymbol{B}}^{*}$ has rank $r^{*}<J^{*}$, which means that ${\boldsymbol{B}}^{*}={\boldsymbol{C}}_{1}^{*}\left({\boldsymbol{C}}_{2}^{*}\right)^{T}$ for some ${\boldsymbol{C}}_{1}^{*}\in\mathbb{R}^{p\times r^{*}}$ and ${\boldsymbol{C}}_{2}^{*}\in\mathbb{R}^{m\times r^{*}}$. Then, the model can be rewritten as ${\boldsymbol{Y}}=\left({\boldsymbol{X}}{\boldsymbol{C}}_{1}^{*}\right)\left({\boldsymbol{C}}_{2}^{*}\right)^{T}+{\boldsymbol{E}}$, where ${\boldsymbol{X}}{\boldsymbol{C}}_{1}^{*}$ is composed of $r^{*}$ factors which are linear combinations of $J^{*}$ relevant predictors to explain all response variables. This is an extension of sparse principal component analysis where ${\boldsymbol{X}}$ is an identity matrix (Zou et al.,, 2006; Shen and Huang,, 2008; Witten et al.,, 2009; Johnstone and Lu,, 2009; Ma,, 2013). For more details on jointly rank-deficient and sparse regression, we refer to Bunea et al., (2012), Chen and Huang, (2012), Chen et al., (2012), and She, (2017).

Although joint regularization guarantees effective dimension reduction, it can be challenging to control the amount of shrinkage adaptively in real data. Let $\lambda$ denote the regularization parameter(s) to be tuned. In the literature, a grid search is often used, which calls a sparse learning algorithm for every value of $\lambda$ in a pre-specified grid. So the problem boils down to the design of an appropriate model comparison criterion.

There are two broad classes of comparison criteria: information criteria, appending various penalties on the model complexity to the training error, and cross-validation (CV), which is based on data resampling. These criteria offer data-dependent (or adaptive) parameters, in contrast to some theoretical choices that are derived assuming incoherent designs.

For ordinary variable selection ($m=1$), some examples of information criteria are AIC (Akaike,, 1974), BIC (Schwarz et al.,, 1978), RIC (Foster and George,, 1994) and EBIC (Chen and Chen,, 2008). We refer to Shao, (1997) and Yang, (2005) for asymptotic studies under the classical regime with $n\rightarrow\infty$ and $p$ fixed. In practical data analysis, however, there appears to be no clear conclusion of which criterion to use, and some results even seem contradictory. There is much less theoretical work on rank selection with $p\gg n$ (Anderson,, 1999; Bunea et al.,, 2011), let alone for joint variable selection and rank reduction. Instead of trying to figure out the appropriate information criterion penalty in a specific setting, many practitioners prefer general-purpose $K$-fold cross-validation (Geisser,, 1975). This requires us to split the dataset randomly into $K$ subsets of roughly equal size, then to fit, say, a lasso model (Tibshirani,, 1996) at any given value of $\lambda$ on the data without the $k$-th subset ($1\leq k\leq K$), and to evaluate prediction error on the left-out subset. We refer to Arlot and Celisse, (2010) for a modern survey of cross-validation. It is a common belief that the summarized cross-validation error provides a good index of a model’s goodness of fit. The choice of $K$ is often discretionary but can have a significant effect.

The main contributions of this paper are three-fold. First, we argue that the $K$ trainings in the conventional cross-validation may be associated with inconsistent models, especially when a nonconvex penalty is in use, and we introduce a structural cross-validation (SCV) based on selection-projection patterns. Second, we develop a new class of predictive information criteria (PIC) that can achieve the minimax optimal error rate without any design incoherence or signal strength conditions. Third, we show that $K$-fold cross-validation is not rate-optimal in pursuing sparsity, and we propose a novel scale-free calibration of the CV error to match the optimal rate. Experiments show the superb performance of the proposed methods.

The majority of this work will assume that the model is jointly row-sparse and rank-deficient. The benefits will be seen, for example, in revealing the additive relationship between degrees of freedom and inflation, which may be difficult to perceive in single-response regression. But all of our results and discussions apply to pure variable selection and low-rank modeling.

The following notation and symbols will be used. Given ${\boldsymbol{B}}\in\mathbb{R}^{p\times m}$, we shall use $r({\boldsymbol{B}})$ to denote the rank of ${\boldsymbol{B}}$, $\mathcal{J}({\boldsymbol{B}})=\{1\leq j\leq p:\|{\boldsymbol{B}}[j,]\|_{2}\neq 0\}$ to denote the row support of ${\boldsymbol{B}}$, and $J({\boldsymbol{B}})=|\mathcal{J}({\boldsymbol{B}})|$. We use ${\boldsymbol{B}}[\mathcal{I},\mathcal{J}]$ to denote a submatrix of ${\boldsymbol{B}}$ with rows and columns indexed by $\mathcal{I},\mathcal{J}$ respectively, and abbreviate ${\boldsymbol{B}}[1\!:\!p,\mathcal{J}]$ to ${\boldsymbol{B}}[,\mathcal{J}]$. Let $\mathbb{O}^{m\times n}=\{{\boldsymbol{A}}\in\mathbb{R}^{m\times n}:{\boldsymbol{A}}^{T}{\boldsymbol{A}}={\boldsymbol{I}}\}$ represent the class of orthogonal matrices. The symbol $\gtrsim$ means that the inequality holds up to multiplicative constants (similarly for $\lesssim$). We use $\|{\boldsymbol{A}}\|_{F}$ to denote the Frobenius norm of ${\boldsymbol{A}}$. Finally, ${\boldsymbol{\mathcal{P}}}_{{\boldsymbol{A}}}$ is the orthogonal projection matrix onto the column space of ${\boldsymbol{A}}$, i.e., ${\boldsymbol{\mathcal{P}}}_{{\boldsymbol{A}}}={\boldsymbol{A}}({\boldsymbol{A}}^{T}{\boldsymbol{A}})^{-}{\boldsymbol{A}}^{T}$ where - is the Moore-Penrose inverse; when there is no ambiguity, we also use ${\boldsymbol{\mathcal{P}}}_{{\boldsymbol{A}}}$ to denote the column space of ${\boldsymbol{A}}$.

2 Structural Cross-Validation

When using cross-validation for parameter tuning, one often holds $\lambda$ fixed at a given value across all folds. Some recent theoretical studies partially support this idea. For instance, consider the lasso which minimizes $\|{\boldsymbol{y}}-{\boldsymbol{X}}{\boldsymbol{\beta}}\|_{2}^{2}+\lambda\|{\boldsymbol{\beta}}\|_{1}$ with ${\boldsymbol{X}}\in\mathbb{R}^{n\times p}$, ${\boldsymbol{y}}\in\mathbb{R}^{n}$. Under some incoherent conditions on the design matrix, $\lambda$ has a universal rate choice of $O(\sigma\sqrt{n\log p})$; see, for example,  Bickel et al., (2009). This rate depends only on problem dimensions and makes the lasso achieve nearly optimal prediction error rate with high probability. Similar conclusions exist for singular-value thresholds in low-rank matrix estimation; see, for example, Candès and Plan, (2011) and Bunea et al., (2011). These universal-rate results may lend some support to fixing the regularization parameter(s) in cross-validation, as long as the designs in different trainings are of the same size.

However, the rate choice, as well as the error bounds obtained with high probability, is too crude and restrictive in real-life data analysis. Moreover, the design matrix is required to satisfy, say, the restricted eigenvalue condition (Bickel et al.,, 2009), the restricted isometry property (Candès,, 2008), the sparse Riesz condition (Zhang and Huang,, 2008) or the comparison condition (She,, 2016), all placing stringent requirements on the Gram matrix ${\boldsymbol{X}}^{T}{\boldsymbol{X}}$ and easy to violate in practice. There is yet no theory that is sufficiently fine to take all characteristics of data into consideration to give a precise formula of $\lambda$. The issue is perhaps easier to understand from a Lagrangian perspective. Indeed, the $\lambda$ in lasso is a Lagrangian multiplier of the constrained problem $\min_{\boldsymbol{\beta}}\|{\boldsymbol{y}}-{\boldsymbol{X}}{\boldsymbol{\beta}}\|_{2}^{2}\mbox{ s.t. }\|{\boldsymbol{\beta}}\|_{1}\leq c$. In contrast to the constraint parameter $c$ that is directly defined on ${\boldsymbol{\beta}}$, $\lambda=\lambda({\boldsymbol{X}},{\boldsymbol{y}},c)$, as a dual parameter of $c$, depends on ${\boldsymbol{X}}$ and ${\boldsymbol{y}}$ as well, and this data dependence is not just restricted to problem’s dimensions. Therefore, on different datasets (associated with the same ${\boldsymbol{\beta}}$), the penalty parameter may have to be changed to maintain the same constraint value. Therefore, with $\lambda$ held fixed in CV, the $K$ fitted models may not be consistent and the judgement based on the total CV-error could be spurious. See Figure 1 for an illustration of the issue, where for a given value of $\lambda$ the cardinality of the trained models varies across the folds (for example, at $\lambda=1.2$ the number of non-zeros varies from nine to sixteen). The issue was discussed in She, (2012) and She et al., (2013).

What if we were to cross-validate $c$? This constraint form is exactly what Tibshirani, (1996) used to define the lasso criterion and to call cross-validation on. But the idea puts some limitations on computation, especially when a nonconvex regularization is in use. Solving a constrained problem is often not easy; for example, changing a SCAD penalty to a constrained form results in a much harder problem (Fan and Li,, 2001). Moreover, it is well known that the solution path associated with a nonconvex penalty or constraint is discontinuous, and algorithms may be easily trapped into a local minimum. These matters only exacerbate the training inconsistency issue.

There is a simple idea to address the problem. Instead of finding the optimal regularization parameter, our goal should be to obtain the best coefficient matrix estimate. So we can cross-validate factor loadings instead of regularization parameters to maintain the same model in $K$ trainings and validations. Below we describe how to extract selection-projection patterns in the jointly row-sparse and rank-deficient setup given an arbitrary nonzero estimate ${\boldsymbol{B}}$, where $r=r({\boldsymbol{B}})$, ${\mathcal{J}}=\mathcal{J}({\boldsymbol{B}})$ and $J=J({\boldsymbol{B}})$.

1. a)

Let ${\boldsymbol{S}}({\boldsymbol{B}})={\boldsymbol{I}}[,{\mathcal{J}}]$, the submatrix formed by the columns indexed by $\mathcal{J}$ in a $p\times p$ identity matrix. 2. b)

If $r<J\wedge m$, find ${\boldsymbol{U}}({\boldsymbol{B}})\in\mathbb{O}^{{J}\times r}$ spanning the column space of ${\boldsymbol{B}}[{\mathcal{J}},]$ (which can be obtained by SVD or QR decomposition). Otherwise, set ${\boldsymbol{U}}({\boldsymbol{B}})={\boldsymbol{I}}$, the identity matrix of size $J\times J$. 3. c)

The structural pattern is given by the orthogonal matrix ${\boldsymbol{S}}({\boldsymbol{B}}){\boldsymbol{U}}({\boldsymbol{B}})$.

According to the procedure, when ${\boldsymbol{B}}$ has row sparsity, the index set of nonzero rows, ${\mathcal{J}}$, determines the sparsity pattern. If ${\boldsymbol{B}}[{\mathcal{J}},]$ has further rank deficiency, its range provides the projection pattern. This can be seen by writing ${\boldsymbol{X}}{\boldsymbol{B}}$ as

[TABLE]

By construction, all structural patterns ${\boldsymbol{S}}({\boldsymbol{B}}){\boldsymbol{U}}({\boldsymbol{B}})$ are orthogonal.

In our jointly sparse model, ${\boldsymbol{X}}{\boldsymbol{S}}({\boldsymbol{B}}){\boldsymbol{U}}({\boldsymbol{B}})$ is a new design matrix that has only $r$ factors constructed from $J$ raw predictors. Explicitly extracting the factors removes estimation bias in regularized estimation. Given $\mathcal{J}$ and $r$, let ${\boldsymbol{B}}_{ub}$ denote the unbiased estimate that solves

[TABLE]

Assume ${\boldsymbol{B}}_{ub}\neq\boldsymbol{0}$ and let ${\boldsymbol{S}}{\boldsymbol{U}}$ be the associated structural pattern of ${\boldsymbol{B}}_{ub}$. Then it is not difficult to prove that ${\boldsymbol{B}}_{ub}$ can always be decomposed into ${\boldsymbol{S}}{\boldsymbol{U}}\cdot{\boldsymbol{C}}$, with ${\boldsymbol{C}}$ the ordinary least squares (OLS) estimate from regressing ${\boldsymbol{Y}}$ on ${\boldsymbol{X}}{\boldsymbol{S}}{\boldsymbol{U}}$.

The procedure of cross-validating the structural patterns (or s-patterns, for short) that are associated with each estimate then follows, which is termed structural cross-validation (SCV).

Split the row indices of $({\boldsymbol{Y}},{\boldsymbol{X}})$ into $K$ subsets. Denote the $k$th ($1\leq k\leq K$) one of the data by $({\boldsymbol{Y}}^{k},{\boldsymbol{X}}^{k})$, and the remaining by $({\boldsymbol{Y}}^{-k},{\boldsymbol{X}}^{-k})$. 2. 2.

For each s-pattern ${\boldsymbol{S}}{\boldsymbol{U}}$, compute $\hat{\boldsymbol{C}}^{-k}$, which is the OLS estimate from regressing ${\boldsymbol{Y}}^{-k}$ on ${\boldsymbol{X}}^{-k}{\boldsymbol{S}}{\boldsymbol{U}}$. Let $\hat{\boldsymbol{B}}^{-k}={\boldsymbol{S}}{\boldsymbol{U}}\hat{\boldsymbol{C}}^{-k}$. Evaluate the prediction error of $\hat{\boldsymbol{B}}^{-k}$ on $({\boldsymbol{Y}}^{k},{\boldsymbol{X}}^{k})$. Repeat the process for all $k:1\leq k\leq K$. 3. 3.

Summarize the total validation error for each candidate model.

Clearly, the $K$ models in SCV are comparable because a pattern only acts on the columns of ${\boldsymbol{X}}$ and remains consistent across the (re-sampled) rows. The fitting step is cost-effective—indeed, each training fits only a low-dimensional OLS model restricted to the selected and/or projected dimensions. This amounts to sparse low-rank estimation with bias calibration in view of (2.1). In comparison, ordinary $K$-fold cross-validation runs a more involved learning algorithm for $K$ times. In structural cross-validation, the learning algorithm is called only once on the overall data to generate candidate patterns.

It is worth mentioning that SCV is applicable to any sparsity-inducing penalty, including those that are nonconvex. It also applies to sole variable selection or rank selection and can be extended to non-Gaussian models such as generalized linear models; see She, (2012).

3 In-sample Optimal Complexity Penalty

Another problem is if $K$-fold SCV is rate optimal or if it applies a sufficiently large penalty for overly complex models. Some popular choices are $K=5$, $10$, or $n$. Under fixed $p$ and $n\rightarrow\infty$, some asymptotic conclusions have been obtained (Shao,, 1997; Yang,, 2005), but it is difficult to determine how large $n$ should be to apply the results in practice. Noticing that the problem is similar to picking an appropriate complexity-penalty in information criteria, this section studies a given model’s in-sample error (Hastie et al.,, 2009). In-sample error is a form of prediction error conditioned on the design matrix (Hastie et al.,, 2009). Our goal is to find a non-asymptotic information criterion to achieve the optimal error rate in prediction without placing any signal-to-noise or design incoherence restriction or requiring the true model to be among the candidate models. The predictive learning principle allows the data to inform the best parsimonious model.

Assume ${\boldsymbol{Y}}={\boldsymbol{X}}{\boldsymbol{B}}^{*}+{\boldsymbol{E}}$ with $\,\mathbb{E}[{\boldsymbol{E}}{\boldsymbol{E}}^{T}]=m\sigma^{2}{\boldsymbol{I}}$. Given an arbitrary s-pattern ${\boldsymbol{S}}{\boldsymbol{U}}$ with $\bar{r}=r({\boldsymbol{S}}{\boldsymbol{U}})$, the associated restricted estimate is $\hat{\boldsymbol{B}}={\boldsymbol{S}}{\boldsymbol{U}}({\boldsymbol{X}}_{{\boldsymbol{S}}{\boldsymbol{U}}}^{T}{\boldsymbol{X}}_{{\boldsymbol{S}}{\boldsymbol{U}}})^{-1}{\boldsymbol{X}}_{{\boldsymbol{S}}{\boldsymbol{U}}}^{T}{\boldsymbol{Y}}$, where we write ${\boldsymbol{X}}{\boldsymbol{S}}{\boldsymbol{U}}$ as ${\boldsymbol{X}}_{{\boldsymbol{S}}{\boldsymbol{U}}}$ for simplicity. The in-sample test error for an independent copy ${\boldsymbol{Y}}^{\prime}$ of ${\boldsymbol{Y}}$ can be computed as follows

[TABLE]

It appears that the criterion

[TABLE]

could be used for model comparison. When $m=1$, $\bar{r}=J$, (3.2) is the AIC (but not when $m>1$). Nevertheless, the above derivation only shows the unbiasedness of (3.2) for a given model. Although we are searching in a $p$-dimensional space of all possible models, the model complexity penalty $2\sigma^{2}\bar{r}m$ does not involve $p$. In fact, minimizing (3.2) will not lead to an estimator with the lowest prediction error in sparse scenarios.

To motivate the correct complexity penalty, we study the minimax error rate for the class of models having both row sparsity and low rank. Given $J$ and $r$, define a signal class

[TABLE]

Any model in (3.3) is jointly parsimonious in its rank and row support. Let $l(\cdot)$ be an arbitrary nondecreasing loss function with $l(0)=0,l\not\equiv 0$. Let $q=r({\boldsymbol{X}})$. Define

[TABLE]

We also write $P_{o}(J,r)$ because $P_{o}$ depends on ${\boldsymbol{B}}$ through $J({\boldsymbol{B}})$ and $r({\boldsymbol{B}})$.

Theorem 1**.**

Let ${\boldsymbol{Y}}={\boldsymbol{X}}{\boldsymbol{B}}^{*}+{\boldsymbol{E}}$ where ${\boldsymbol{E}}=[e_{ik}]$ with $e_{ik}$ $\overset{i.i.d.}{\sim}\mathcal{N}(0,\sigma^{2})$, and ${\boldsymbol{B}}^{*}\in S(J,r)$ for some $J,r$ satisfying $1\leq J\leq p/2$, $p\geq 2$, $r(q\wedge J+m-r)\geq 16$, $1\leq r\leq\min\{(J+m)/2,2J\}$. Assume the following restricted condition-number condition: (i) when $J\leq q$, there exist $\underline{\kappa},\overline{\kappa}>0$ such that $\underline{\kappa}\|{\boldsymbol{B}}\|_{F}^{2}\leq\|{\boldsymbol{X}}{\boldsymbol{B}}\|_{F}^{2}\leq\overline{\kappa}\|{\boldsymbol{B}}\|_{F}^{2}$ for any ${\boldsymbol{B}}:J({\boldsymbol{B}})\leq J,r({\boldsymbol{B}})\leq r$ and $\underline{\kappa}/\overline{\kappa}$ is a positive constant; (ii) when $J>q$, $\sigma_{\min}^{2}({\boldsymbol{X}})/\sigma_{\max}^{2}({\boldsymbol{X}})$ is a positive constant, where $\sigma_{\min}({\boldsymbol{X}})$ and $\sigma_{\max}({\boldsymbol{X}})$ are the smallest and the largest nonzero singular values of ${\boldsymbol{X}}$, respectively. Then there exist positive constants $c,c^{\prime}$, depending on $l(\cdot)$ only, such that when $J\leq q$,

[TABLE]

and when $J>q$, $\inf_{\hat{\boldsymbol{B}}}\sup_{{\boldsymbol{B}}^{*}\in S(J,r)}\,\mathbb{E}[l(\|{\boldsymbol{X}}{\boldsymbol{B}}^{*}-{\boldsymbol{X}}\hat{\boldsymbol{B}}\|_{F}^{2}/(c\sigma^{2}(q+m-r)r))]\geq c^{\prime}>0.$

Notice that the restricted condition-number condition takes a more relaxed form when $J$ is large ($J>q$) and allows $p$ to be greater than $n$. Setting $l(t)=t$, the theorem shows that the minimax optimal rate for the risk $\,\mathbb{E}\|{\boldsymbol{X}}{\boldsymbol{B}}^{*}-{\boldsymbol{X}}\hat{\boldsymbol{B}}\|_{F}^{2}$ is $\sigma^{2}P_{o}(J,r)$. Also, with $l(t)=1_{t\geq 1}$, for any estimator $\hat{\boldsymbol{B}}$, there exists a ${\boldsymbol{B}}^{*}\in S(J,r)$ such that $\|{\boldsymbol{X}}{\boldsymbol{B}}^{*}-{\boldsymbol{X}}\hat{\boldsymbol{B}}\|_{F}^{2}\gtrsim\sigma^{2}P_{o}(J,r)$ with positive probability. The minimax lower bound suggests $P_{o}$ be an ideal model-complexity penalty, which leads to a new information criterion for model comparison.

Theorem 2**.**

Assume ${\boldsymbol{Y}}={\boldsymbol{X}}{\boldsymbol{B}}^{*}+{\boldsymbol{E}}$ where $\operatorname{\mbox{vec}\,}({\boldsymbol{E}})$ is sub-Gaussian with mean $\boldsymbol{0}$ and scale bounded by $\sigma$. For any collection of (possibly random) non-zero matrices ${\boldsymbol{B}}_{1}$, ${\boldsymbol{B}}_{2}$, $\ldots$, select the optimal one by

[TABLE]

with the constant $A$ appropriately large. Then $\hat{\boldsymbol{B}}$ satisfies

[TABLE]

We refer to this new information criterion as the Predictive Information Criterion (PIC). (3.6) does not require ${\boldsymbol{B}}_{l}$ ($l\geq 1$) to cover the true model. But when ${\boldsymbol{B}}^{*}$ is among the candidate matrices or there exists some ${\boldsymbol{B}}_{l}$ close to ${\boldsymbol{B}}^{*}$, the risk of $\hat{\boldsymbol{B}}$ is bounded above by $\sigma^{2}P_{o}\left({\boldsymbol{B}}^{*}\right)$ up to multiplicative constants. When $m=1$, the complexity penalty in PIC is of the order $\sigma^{2}J\log(ep/J)$, is similar to (but finer than) the rate in RIC (Foster and George,, 1994). However, we do not assume that $n-p$ must be large as in the derivation of RIC. In fact, PIC is non-asymptotic in nature, and requires no large-$n$ assumption. Moreover, it does not require any incoherence condition which is commonly assumed in the literature. From the bias-variance tradeoff in (3.6), setting ${\boldsymbol{B}}_{l}={\boldsymbol{B}}^{*}$, the noise-free truth, may not yield the most accurate and parsimonious model, especially when the noise level is large. This is well known in wavelet studies (Donoho and Johnstone,, 1994). When the signal is properly large, minimizing PIC also recovers the row support of ${\boldsymbol{B}}^{*}$ with high probability under some regularity conditions; see Theorem 5 in the Appendix.

In (3.4), $[J({\boldsymbol{B}})+m-r({\boldsymbol{B}})]r({\boldsymbol{B}})$ corresponds to the number of free parameters of the model or degrees of freedom (DF), typically less than $mn$, while $J({\boldsymbol{B}})\log(ep/J({\boldsymbol{B}}))$ characterizes the inflation (IF) due to the selection among $p$ predictors. The model-complexity penalty can then be written as $A_{1}\sigma^{2}\mbox{{DF}}+A_{2}\sigma^{2}\mbox{{IF}}$, where $A_{1},A_{2}$ are constants. The additive form is in contrast to the multiplicative form $c(n,p)\times\sigma^{2}\mbox{{DF}}$ that is widely used in the literature, for example, $c(n,p)=\log p$ in BIC.

The unknown noise scale remains an issue in supervised learning, because in large-$p$ settings estimating $\sigma$ could be as challenging as estimating the mean. The following theorem, as a generalization of Theorem 3 in She, (2017), presents four scale-free forms of PIC to bypass the scale estimation under a model sparsity assumption.

Theorem 3**.**

Let ${\boldsymbol{Y}}={\boldsymbol{X}}{\boldsymbol{B}}^{*}+{\boldsymbol{E}}$ where ${\boldsymbol{E}}=[e_{ik}]$ with $e_{ik}$ $\overset{i.i.d.}{\sim}\mathcal{N}(0,\sigma^{2})$. Suppose the true model is parsimonious in the sense that $P_{o}({\boldsymbol{B}}^{*})\allowbreak<mn/A_{0}$ for some constant $A_{0}>0$. Let $\delta({\boldsymbol{B}})=AP_{o}({\boldsymbol{B}})/(mn)$ for some constant $A<A_{0}$, and so $\delta({\boldsymbol{B}}^{*})<1$. Consider the following criteria

[TABLE]

Then, for sufficiently large values of $A_{0},A$, any $\hat{\boldsymbol{B}}$ that minimizes (3.7), (3.8), (3.9), or (3.10) subject to $\delta({\boldsymbol{B}})<1$ satisfies $\|{\boldsymbol{X}}\hat{\boldsymbol{B}}-{\boldsymbol{X}}{\boldsymbol{B}}^{*}\|_{F}^{2}\lesssim\sigma^{2}P_{o}({\boldsymbol{B}}^{*})$ with probability at least $1-Cp^{-c}-C^{\prime}\exp(-c^{\prime}mn)$ for some constants $C,C^{\prime},c,c^{\prime}>0$.

We call (3.7), (3.8), (3.9), (3.10) the fractional form, GCV form, logarithmic form, and plug-in form of PIC, respectively. With the inflation term removed, (3.8) shares similarity with GCV (Wahba,, 1990); the log form is commonly seen when applying AIC or BIC with an unknown $\sigma^{2}$; the penalty in (3.10) can be written as $A\hat{\sigma}^{2}P_{o}({\boldsymbol{B}})$ with $\hat{\sigma}^{2}=\|{\boldsymbol{Y}}-{\boldsymbol{X}}{\boldsymbol{B}}\|_{F}^{2}/(mn),$ which resembles Mallows’ $C_{p}$ (Mallows,, 1973). The sparsity assumption and constraint cannot be dropped, which in turn rule out over-complex models.

In common with most non-asymptotic analyses, we showed the optimal rate but not the optimal numerical constants. The absolute constants can be determined by Monte Carlo experiments—for example, in the fractional form we recommend using $\|{\boldsymbol{Y}}-{\boldsymbol{X}}{\boldsymbol{B}}\|_{F}^{2}/[1-(2\cdot\mbox{DF}+1.8\cdot\mbox{IF})/(mn)]$.

4 Rate Calibration of Cross-validation

This section studies extra-sample error that cross-validation methods try to estimate (Hastie et al.,, 2009, Chapter 7). It turns out that there is a connection between the in-sample error and the extra-sample error, which can be used to guide cross-validation.

To define the extra-sample error, we assume that the row observations in the design are i.i.d. and independent of the noise component. Specifically, assume ${\boldsymbol{Y}}={\boldsymbol{X}}{\boldsymbol{B}}^{*}+{\boldsymbol{E}}$, where ${\boldsymbol{E}}=[e_{ik}]$ satisfies $\,\mathbb{E}[{\boldsymbol{E}}{\boldsymbol{E}}^{T}]=\sigma^{2}m{\boldsymbol{I}}$, ${\boldsymbol{X}}=[\tilde{\boldsymbol{x}}_{1}\,\cdots\,\tilde{\boldsymbol{x}}_{n}]^{T}$ is independent of ${\boldsymbol{E}}$ and has i.i.d. rows with ${\boldsymbol{\Sigma}}$ (positive-definite) as the covariance matrix. The $s$-pattern from ${\boldsymbol{B}}^{*}$ is denoted by ${\boldsymbol{S}}^{*}{\boldsymbol{U}}^{*}$. Let ${\boldsymbol{S}}{\boldsymbol{U}}$ be a given candidate s-pattern with $\bar{r}=r({\boldsymbol{S}}({\boldsymbol{B}}){\boldsymbol{U}}({\boldsymbol{B}}))$. Recall the training error and the structural cross-validation error

[TABLE]

where $\hat{\boldsymbol{B}}$, $\hat{\boldsymbol{B}}^{-{k}}$ are the restricted OLS estimates associated with the given s-pattern, obtained on the overall data and the data without the $k$-th subset, respectively. Concretely,

[TABLE]

with ${\boldsymbol{X}}_{{\boldsymbol{S}}{\boldsymbol{U}}}={\boldsymbol{X}}{\boldsymbol{S}}{\boldsymbol{U}}$ and ${\boldsymbol{X}}_{{\boldsymbol{S}}{\boldsymbol{U}}}^{-k}={\boldsymbol{X}}^{-k}{\boldsymbol{S}}{\boldsymbol{U}}$, where all the restricted OLS problems are assumed to be non-degenerate. Suppose that $n=dK$ for some integer $d$ (and so $d=n/K$). The theorem below gives an identity of the cross-validation error.

Theorem 4**.**

Given any s-pattern ${\boldsymbol{S}}{\boldsymbol{U}}$ with $\bar{r}=r({\boldsymbol{S}}{\boldsymbol{U}})$, the following identity holds

[TABLE]

where

[TABLE]

and

[TABLE]

*Here, ${\boldsymbol{Z}}_{n-d}\in\mathbb{R}^{(n-d)\times\bar{r}}$ is a submatrix composed of (any) $n-d$ rows of ${\boldsymbol{Z}}={\boldsymbol{X}}_{{\boldsymbol{S}}{\boldsymbol{U}}}({\boldsymbol{S}}{\boldsymbol{U}})^{T}\allowbreak{\boldsymbol{\Sigma}}^{-1/2}{\boldsymbol{S}}{\boldsymbol{U}}$, ${\boldsymbol{\Sigma}}_{{\boldsymbol{S}}{\boldsymbol{U}}}=({\boldsymbol{S}}{\boldsymbol{U}})^{T}(n{\boldsymbol{\Sigma}}){\boldsymbol{S}}{\boldsymbol{U}}$, ${\boldsymbol{B}}_{{\boldsymbol{S}}{\boldsymbol{U}},{\boldsymbol{S}}^{*}{\boldsymbol{U}}^{*}}$ is $(({\boldsymbol{X}}_{{\boldsymbol{S}}{\boldsymbol{U}}}^{-k})^{T}{\boldsymbol{X}}_{{\boldsymbol{S}}{\boldsymbol{U}}}^{-k})^{-1}\allowbreak({\boldsymbol{X}}_{{\boldsymbol{S}}{\boldsymbol{U}}}^{-k})^{T}({\boldsymbol{X}}^{-k}{\boldsymbol{\mathcal{P}}}^{\circ}{\boldsymbol{B}}^{*})$ with say $k=1$, $\bar{\boldsymbol{B}}_{{\boldsymbol{S}}{\boldsymbol{U}},{\boldsymbol{S}}^{*}{\boldsymbol{U}}^{*}}=({\boldsymbol{X}}_{{\boldsymbol{S}}{\boldsymbol{U}}}^{T}{\boldsymbol{X}}_{{\boldsymbol{S}}{\boldsymbol{U}}})^{-1}{\boldsymbol{X}}_{{\boldsymbol{S}}{\boldsymbol{U}}}^{T}\allowbreak({\boldsymbol{X}}{\boldsymbol{\mathcal{P}}}^{\circ}{\boldsymbol{B}}^{*})$, and ${\boldsymbol{\mathcal{P}}}^{\circ}$ is the orthogonal projection onto ${\boldsymbol{\mathcal{P}}}_{{\boldsymbol{S}}^{*}{\boldsymbol{U}}^{*}}\setminus({\boldsymbol{\mathcal{P}}}_{{\boldsymbol{S}}^{*}{\boldsymbol{U}}^{*}}\cap{\boldsymbol{\mathcal{P}}}_{{\boldsymbol{S}}{\boldsymbol{U}}})$, i.e., the orthogonal complement of ${\boldsymbol{\mathcal{P}}}_{{\boldsymbol{S}}^{*}{\boldsymbol{U}}^{*}}\cap{\boldsymbol{\mathcal{P}}}_{{\boldsymbol{S}}{\boldsymbol{U}}}$ in the subspace ${\boldsymbol{\mathcal{P}}}_{{\boldsymbol{S}}^{*}{\boldsymbol{U}}^{*}}$. *

According to (4.3), $U$ is always nonnegative. In the classical asymptotic regime where the sample size tends to $\infty$ and $m,p$ are fixed, ${\boldsymbol{B}}_{{\boldsymbol{S}}{\boldsymbol{U}},{\boldsymbol{S}}^{*}{\boldsymbol{U}}^{*}}$, as well as $\bar{\boldsymbol{B}}_{{\boldsymbol{S}}{\boldsymbol{U}},{\boldsymbol{S}}^{*}{\boldsymbol{U}}^{*}}$, will approach $\,\mathbb{E}{{\boldsymbol{B}}_{{\boldsymbol{S}}{\boldsymbol{U}},{\boldsymbol{S}}^{*}{\boldsymbol{U}}^{*}}}$ because of the law of large numbers. In finite samples, however, $U$ incurs extra cost for underfitting models, which could be an advantage over some information criteria. From the definition of ${\boldsymbol{\mathcal{P}}}^{\circ}$, $U$ vanishes for over-fitting models that satisfy ${\boldsymbol{B}}^{*}\in{\boldsymbol{\mathcal{P}}}_{{\boldsymbol{S}}{\boldsymbol{U}}}$, and hence is not active for “large" models.

The term $D$ penalizes the cardinality of the model, and relates to the degrees of freedom. The definition of $D$ is based on ${\boldsymbol{Z}}$ which is just ${\boldsymbol{X}}{{\boldsymbol{S}}{\boldsymbol{U}}}$ decorrelated. Because the rows of ${\boldsymbol{Z}}$ are isotropic (i.e., the covariance of each row vector is ${\boldsymbol{I}}$), it is not difficult to show an upper bound of $D$ using random matrix theory. We give a corollary as an illustration under the assumption that the rows of ${\boldsymbol{X}}$ are i.i.d. $\mathcal{N}(\boldsymbol{0},{\boldsymbol{\Sigma}})$ (such designs are widely used in simulation studies).

Corollary 1**.**

Under the setup of the previous theorem and the Gaussian assumption of the design,

[TABLE]

where $U$ is as defined in Theorem 4. In particular, for ordinary variable selection with $m=1$, the identity for any given support $\mathcal{J}$ reads

[TABLE]

where $J=|\mathcal{J}|$, ${\boldsymbol{\Sigma}}_{\mathcal{J}}=n{\boldsymbol{\Sigma}}[\mathcal{J},\mathcal{J}]$, ${\boldsymbol{X}}_{\mathcal{J}}={\boldsymbol{X}}[\mathcal{J},]$, $\bar{\boldsymbol{\beta}}_{\mathcal{J},\mathcal{J}^{*}}=({\boldsymbol{X}}_{\mathcal{J}}^{T}{\boldsymbol{X}}_{\mathcal{J}})^{-1}{\boldsymbol{X}}_{\mathcal{J}}^{T}{\boldsymbol{X}}_{\mathcal{J}^{*}\cap\mathcal{J}^{c}}\allowbreak{\boldsymbol{\beta}}_{{\mathcal{J}^{*}\cap\mathcal{J}^{c}}}^{*}$ and ${\boldsymbol{\beta}}_{\mathcal{J},\mathcal{J}^{*}}=(({\boldsymbol{X}}_{\mathcal{J}}^{-k})^{T}\allowbreak{\boldsymbol{X}}_{\mathcal{J}}^{-k})^{-1}({\boldsymbol{X}}_{\mathcal{J}}^{-k})^{T}{\boldsymbol{X}}_{\mathcal{J}^{*}\cap\mathcal{J}^{c}}^{-k}{\boldsymbol{\beta}}_{{\mathcal{J}^{*}\cap\mathcal{J}^{c}}}^{*}$.

In (4.4), $D$ is a decreasing function in $K$. When $K=n$, $D$ is essentially $2m\bar{r}\sigma^{2}$; even when $K$ is as small as $2$, $D\leq\{3+4\theta/(1-2\theta)\}m\bar{r}\sigma^{2}$ under $(\bar{r}+1)/n\leq\theta$, and so for $n$ large, the CV-Err of an over-fitting model is no larger than

[TABLE]

This means that the extra-sample error of $K$-fold cross-validation is not significantly different from that in the in-sample error formula (3.1), which is, however, perhaps natural.

Returning to the jointly sparse model (where $\bar{r}=r({\boldsymbol{S}}({\boldsymbol{B}}){\boldsymbol{U}}({\boldsymbol{B}}))=r({\boldsymbol{B}})=r$), combining the identity with PIC gives a choice of $K$ to match the minimax rate in (3.4): $K=\{AP_{o}(J,r)-mr\}/\{AP_{o}(J,r)-2mr-(AP_{o}(J,r)-mr)(r+1)/n\}$. In large-$p$ problems, however, the value is below $2$ and unattainable; consequently, $K$-fold cross-validation ($K\geq 2$) cannot penalize the complexity of the model sufficiently heavily.

One possible fix is to use delete-$d$ cross-validation (Shao,, 1993) which removes $d$ observations in each training and $d$ can be larger than $n/2$. A similar identity to that in Theorem 4 holds for this form of cross-validation, and to match the PIC rate, the training sample size $n-d$ should be of order ${mnr}/(AP_{o}(J,r)-mr)+r+1$ or $O(n\cdot{mr}/({Jr+J\log p}))$. Delete-$d$ cross-validation can be implemented by enumerating all subsets of size $d$, or in a stochastic fashion by randomly splitting the whole dataset many times. But neither is computationally efficient for large $n$ and we shall not pursue further in this work.

By contrast, the commonly used 5-fold CV and 10-fold CV are much less expensive, and these fold choices enjoy small variance (Hastie et al.,, 2009). Although $K$-fold CV ($2\leq K\leq n$) is suboptimal as seen from (4.5) and (3.4), we can make a rate calibration as in the in-sample error case. Let $\mbox{R}=(q\wedge J-r)r$ and recall $\mbox{IF}=J\log(ep/J)$. From Theorems 2 and 4, if $\sigma^{2}$ is known, one can append a bias correction term of the order $\sigma^{2}\mbox{R}+\sigma^{2}\mathrm{IF}$ to the CV-Err. When $\sigma$ is unknown, motivated by the plug-in form of the $\sigma$-free PIC (cf. Theorem 3), we can use the following calibrated structural cross-validation error as the model selection criterion:

[TABLE]

where $\alpha_{1},\alpha_{2}$ are constants. The rate correction is evident from the identity: the calibrated CV error (4.6) has a complexity penalty of the form $\sigma^{2}\mbox{{DF}}+\sigma^{2}\mbox{{IF}}$, complying with the studies in Section 3. According to the constraint $\delta({\boldsymbol{B}})<1$ in Theorem 3, the candidate models that are far too complex (say $\alpha_{1}\mbox{DF}+\alpha_{2}\mbox{IF}>mn$) should be excluded; equivalently, we set their SCV-Err to be $+\infty$. For $K=5$, we recommend $\alpha_{1}=4.6$ and $\alpha_{2}=3.5$ on the basis of Monte Carlo experiments. (The rate correction of SCV applies to any choice of $K$, but the numerical constants may be different.) Of course, other forms are possible—for example, the fractional form, $\mbox{CV-Err}/[1-\alpha_{1}\mbox{R}/(mn+\alpha mr)-\alpha_{2}\mbox{IF}/(mn+\alpha mr)]$, or simply

[TABLE]

with $\alpha_{1}=2$ and $\alpha_{2}=2.4$ also gives satisfactory performance. We apply the plug-in form (4.6) as the SCV error in the regression setting.

5 Experiments

5.1 Simulations

We use the following setup for the simulation studies. The design matrix ${\boldsymbol{X}}$ has i.i.d. rows from $\mathcal{N}\left(\boldsymbol{0},{\boldsymbol{\Sigma}}\right)$ with $\Sigma_{jk}=\rho^{|j-k|}$, $\rho>0$, $1\leq j$, $k\leq p$. The coefficient matrix has the form ${\boldsymbol{B}}^{*}=[b({\boldsymbol{A}}_{0}{\boldsymbol{A}}_{1})^{T}\,\,\boldsymbol{0}]^{T}$ where $b$ is a constant, ${\boldsymbol{A}}_{0}$ is a $J\times r$ matrix and ${\boldsymbol{A}}_{1}$ is an $r\times m$ matrix. Entries in ${\boldsymbol{A}}_{0}$ and ${\boldsymbol{A}}_{1}$ are i.i.d. standard Gaussian. Entries in ${\boldsymbol{B}}^{*}$ past the $J$th row are all zero. The matrix ${\boldsymbol{E}}$ has i.i.d. standard Gaussian entries and the response matrix is ${\boldsymbol{Y}}={\boldsymbol{X}}{\boldsymbol{B}}^{*}+\sigma{\boldsymbol{E}}$.

We consider the following cases of $p>n$ and $n>p$, with two different correlation levels and signal strengths in generating the design ${\boldsymbol{X}}$:

$n>p$: $n=100$, $J=30$, $p=60$, $m=15$, $r=5$, $\sigma=1$, $\rho=0.1,0.5$, $b=0.1,0.5$. 2. 2.

$p>n$: $n=30$, $J=15$, $p=100$, $m=10$, $r=2$, $\sigma=1$, $\rho=0.1,0.5$, $b=0.2,1$.

We ran $200$ simulations for each setting. In common with other studies on high-dimensional parameter tuning (Chen and Chen,, 2008), we called a learning algorithm (cf. Bunea et al., (2012) and Chen and Huang, (2012)) on each synthetic dataset to compute a solution path of candidate estimates. Then the associated projection-selection patterns were extracted according to the procedure in Section 2.

Six model selection criteria were compared: AIC, BIC, EBIC, all taking the logarithmic form due to the unknown noise level, PIC as in (3.7), $K$-fold CV with $K=2,5,10$, and $5$-fold SCV—see (4.6). (We found that the performance of $K$-fold SCV was similar for $K=2,5,10$, and so focused on $5$-fold SCV.) The prediction accuracy was evaluated by the mean squared error (MSE): $\,\mathbb{E}[Tr\{(\hat{\boldsymbol{B}}-{\boldsymbol{B}}^{*})^{T}{\boldsymbol{\Sigma}}(\hat{\boldsymbol{B}}-{\boldsymbol{B}}^{*})\}]/m$. The median of the MSEs over all simulations was then computed to represent the goodness of fit of a model. Selection performance was assessed by the median number of predictors ($\hat{J}$) and the median rank estimate ($\hat{r}$) over all simulations. The rates of left-out noise-free true variables (M for missing) and incorrectly included variables (FA for false alarms) averaged across all runs are also reported. When the signal strength is high, a successful variable selection method minimizes the M- and FA-rates with a preference for a low miss rate since it is undesirable to leave out true features. In contrast, when the noise level is very high, there is no reason to believe that the noise-free simulation truth yields the best predictive model from the observed data. Therefore, in the low SNR situations our ultimate concern is the prediction error. The results for the six model comparison criteria are presented in Tables 1–4 and are summarized below.

The PIC and SCV methods have superior prediction performance across nearly all combinations of signal strength and correlation considered here. EBIC’s prediction in the lower signal-to-noise ratio (SNR) experiments (see, for example, Table 1) is evidence that its comparatively high complexity penalty induces too much regularization for these particular setups. In the larger SNR experiments (cf. Tables 2 and 4), the EBIC’s prediction performance is comparable to PIC and SCV. In the $p>n$ cases, AIC has the highest MSEs among all the information criteria, which is unsurprising because it is well known that AIC underpenalizes for large $p$.

CV often has the highest MSEs among all methods, indicating that a rate correction is absolutely necessary. In all experiments the prediction performance of CV across $K=2,5,10$ is similar, verifying the discussion following Corollary 1. $K=5$ showed a modest improvement compared to $K=2,10$.

Interestingly, the PIC and SCV miss some of the variables specified in the absence of noise when the signal strength is weak (cf. Tables 1 and 3), but this is perhaps natural because a larger amount of regularization may be required to achieve low prediction error. In these small SNR experiments, the AIC tends to miss the fewest of these noise-free variables, but this comes at the cost of higher FA-rates and as stated before it often has poor prediction accuracy. In both weak and strong signal strength situations, the BIC and especially the EBIC tend to have high missing rates compared to the AIC because of their larger penalty terms but, unlike the PIC and SCV, they are rarely able to select a model with parsimony and low prediction error in both weak and strong SNR situations.

In the $n>p$ experiments with relatively stronger signal strength (Table 2), PIC and SCV’s variable selection performance nearly equals or exceeds that of all other methods. In fact, their M- and FA- rates are nearly 0 in Table 2, suggesting that here the true variables are also highly predictive. Indeed, in such a large SNR situation (see Theorem 5) PIC can recover the noise-free row support with high probability.

Compared to the other methods, the median rank values of the PIC and SCV are most consistently equal to their true values. In nearly every experiment the AIC and CV overestimate $r$ while the BIC and especially EBIC underestimate the rank when the signal strength is relatively weak.

The prediction performance of the PIC and SCV was similar in the $n>p$ experiments (cf. Tables 1 and 2), but SCV revealed some advantage in Table 4. We made additional comparisons between these two methods by varying the signal strength in the $p>n$ experiments; see Table 5. SCV produced the lowest standard errors and seemed to be more successful at picking models with the lowest MSE’s, as well as reducing instability.

Overall, in almost all cases either the PIC or SCV has the lowest prediction error, which is the ultimate goal of this work. In the weak signal strength situations, it could be argued the AIC is better at selection because it misses the fewest variables in terms of the noise-free simulation truth. However, as previously discussed, because the low SNR data are heavily contaminated by noise, variable selection as measured by the M and FA-rates may not be that meaningful; a more parsimonious and predictive model may exist that diverges from the zero-noise model used in generating synthetic data.

5.2 Yeast cell cycle data

In an experiment conducted by Spellman et al., (1998), $106$ transcription factors (TFs) (also known as DNA binding proteins) were collected for $800$ yeast genes that regulate RNA levels within the eukaryotic cycle. The cell cycle was measured by taking RNA levels on the $800$ genes at $18$ time points using the $\alpha$ factor arrest method. In this data analysis we use a subset of the original dataset obtained from the R (R Core Team,, 2017) package “spls" (Chun and Keleş,, 2010). The ${\boldsymbol{X}}$ matrix consists of $106$ transcription factors (TF) collected on $542$ genes; the ${\boldsymbol{Y}}$ matrix contains RNA levels measured on the same subset of genes at $18$ time points. For the dataset, there are $21$ experimentally verified TFs related to cell cycle regulation (Wang et al.,, 2007). They can serve as a biological truth and should be consistently picked by a variable selection technique. We centered and scaled both ${\boldsymbol{X}}$ and ${\boldsymbol{Y}}$, performed selective reduced rank regression (She,, 2017) and compared the selection performance of ordinary CV and structural cross-validation (SCV). (The prediction performance of CV and SCV was found to be quite similar over 200 repeated training-test splits with 50% for training and 50% for test.)

We bootstrapped the data $200$ times to measure the stability of the two methods in terms of both rank and variable selection. The bootstrap distributions of the estimated ranks and cardinalities are displayed in Figure 2. SCV’s median $\hat{J}$ and $\hat{r}$ were equal to $86$ and $4$ and CV’s median $\hat{J}$ and $\hat{r}$ were equal to $46$ and $7$, respectively. Clearly, the spread of the optimal ranks and cardinalities selected by CV is much larger than SCV, showing that CV is unstable in terms of model selection. Even though CV’s median $\hat{J}$ is smaller, its large variance suggests that this method picks and leaves out TFs in a more random fashion than does SCV. The number of free parameters is determined jointly by $r$ and $J$; the median of the degrees of freedom was $400$ for SCV and $410$ for CV, suggesting that SCV does tend to select smaller models. In addition, CV showed a five-fold increase in computation time than SCV.

To assess the selection performance in terms of the $21$ experimentally verified TFs, the left panel of Figure 3 shows the percentage of bootstrap replicates in which each TF was selected. Each point corresponds to a TF and the dotted line labels all identical selection frequencies by the two methods. Notably, every TF lies either on or above the line, showing that SCV’s selection frequencies for the confirmed TFs were uniformly larger than CV’s. It is clear that CV often fails to select all of the TFs confirmed to be related to cell cycle regulation. For example, BAS1 is likely to be a significant regulator of the yeast cell cycle (Cokus et al.,, 2006) and was selected by SCV in nearly $75\%$ of the replicates, but CV selected this TF in under $25\%$ of the replicates. The right panel of Figure 3 tabulates the number of confirmed TFs selected at various percentage cut-offs of the bootstrap replicates. The large gap between SCV and CV shows that SCV is far more successful at selecting the verified TFs at any cut-off. For example, roughly $50\%$ of the confirmed TFs are selected by CV in at least $50\%$ of the replicates, but SCV has a nearly $100\%$ success rate at the same cut-off. SCV also identified some TFs that are not part of the confirmed subset. For example, SCV selected SKO1 $194$ times while CV picked this TF less than 50% of the time; Niu et al., (2008) experimentally determined that overexpression of this TF is related to cell cycle progression.

6 Discussion

The conventional approach to cross-validation involves splitting the data into $K$ subsets and calling a learning algorithm $K$ times. This procedure can be expensive and unstable because it may result in $K$ fitted models not directly comparable. As a remedy to this problem, we proposed structural cross-validation, which maintains the same model in data resampling and is computationally efficient. Theoretically, we showed that the optimal complexity rate for joint variable and rank selection is achieved by the predictive information criterion non-asymptotically. Based on an identity built for cross-validation error, we proposed scale-free rate correction for the commonly used $K$-fold cross-validation to match the optimal error rate. To the best of our knowledge, this calibrated cross-validation is novel.

The notions of structural cross-validation and PIC are applicable to either pure variable selection or pure rank selection. For example, in the case of a single response variable, the SCV in (4.6) takes the form ${\mbox{CV-Err}}+\alpha(\mbox{Trn-Err}/n)J\log(ep/J)$ for some positive constant $\alpha$, whereas in pure ridge-type problems the inflation term disappears. The extraction of $s$-patterns extends to generalized lasso problems that pursue sparsity of ${\boldsymbol{T}}{\boldsymbol{\beta}}$ for a thin ${\boldsymbol{T}}$ matrix (She,, 2010; Tibshirani and Taylor,, 2011). For example, given an estimate $\hat{\boldsymbol{\beta}}$ satisfying ${\boldsymbol{T}}[j,]{\boldsymbol{\beta}}\neq 0$ for all $j\in\mathcal{J}$, or ${\boldsymbol{T}}[{\mathcal{J}}^{c},]\hat{\boldsymbol{\beta}}=\boldsymbol{0}$, we can construct a structural pattern ${\boldsymbol{O}}\in\mathbb{O}^{p\times r}$, with $r=p-r({\boldsymbol{T}}[{\mathcal{J}}^{c},])$, that spans the orthogonal complement of the row space of ${\boldsymbol{T}}[{\mathcal{J}}^{c},]$.

This framework pursues prediction accuracy as its ultimate principle and leads to some general theorems without large-$n$ assumption or incoherent design conditions. If the SNR is not too small, adopting the predictive learning perspective also automatically implies faithful support recovery. Indeed, when the noise contamination was relatively small, our simulations showed that the proposed methods had low prediction error and satisfactory selection performance. But our experiments also demonstrated that when the SNR is small, it may not be valid to assess variable selection in reference to the noise-free simulation truth. These studies suggest that the aphorism “all models are wrong but some are useful" (Box,, 1979) seems to apply in the small SNR scenarios.

Appendix A Proof of Theorem 1

We consider three cases as follows.

Case (i): $(J+m-r)r\geq J\log(ep/J)$ and $q\geq J$. Consider a signal subclass

[TABLE]

where $R=\frac{\sigma}{{\overline{\kappa}^{1/2}}}$ and $\gamma>0$ is a small constant to be chosen later. Clearly, $|{\mathcal{B}}^{1}(r)|=2^{(J+m-r/2)r/2}$, $\mathcal{B}^{1}(r)\subset S(J,r)$, and $r({\boldsymbol{B}}_{1}-{\boldsymbol{B}}_{2})\leq r$, $\forall{\boldsymbol{B}}_{1},{\boldsymbol{B}}_{2}\in\mathcal{B}^{1}(r)$. Also, since $r\leq(J+m)/2$, $(J+m-r/2)r/2\geq c(J+m-r)r$ for some universal constant $c$.

Let $\rho({\boldsymbol{B}}_{1},{\boldsymbol{B}}_{2})=\|\operatorname{\mbox{vec}\,}({\boldsymbol{B}}_{1})-\operatorname{\mbox{vec}\,}({\boldsymbol{B}}_{2})\|_{0}$ be the Hamming distance. By the Varshamov-Gilbert bound (cf. Lemma 2.9 in Tsybakov, (2009)), there exists a subset ${\mathcal{B}}^{10}(r)\subset{\mathcal{B}}^{1}(r)$ and $\boldsymbol{0}\in{\mathcal{B}}^{10}(r)$ such that

[TABLE]

and

[TABLE]

for some universal constants $c_{1},c_{2}>0$. Then $\|{\boldsymbol{B}}_{1}-{\boldsymbol{B}}_{2}\|_{F}^{2}=\gamma^{2}R^{2}\rho({\boldsymbol{B}}_{1},{\boldsymbol{B}}_{2})\geq c_{2}\gamma^{2}R^{2}(J+m-r)r$. It follows from the conditional number assumption that

[TABLE]

for any ${\boldsymbol{B}}_{1},{\boldsymbol{B}}_{2}\in\mathcal{B}^{10}$, ${\boldsymbol{B}}_{1}\neq{\boldsymbol{B}}_{2}$, where $\underline{\kappa}/\overline{\kappa}$ is a positive constant.

For Gaussian models, the Kullback-Leibler divergence of $\mathcal{N}(\operatorname{\mbox{vec}\,}({\boldsymbol{X}}{\boldsymbol{B}}_{2}),\sigma^{2}\boldsymbol{I}\otimes{\boldsymbol{I}})$ (denoted by $P_{{\boldsymbol{B}}_{2}}$) from $\mathcal{N}(\operatorname{\mbox{vec}\,}({\boldsymbol{X}}{\boldsymbol{B}}_{1}),\sigma^{2}\boldsymbol{I}\otimes{\boldsymbol{I}})$ (denoted by $P_{{\boldsymbol{B}}_{1}}$) is $\mathcal{K}(\mathcal{P}_{{\boldsymbol{B}}_{1}},\mathcal{P}_{{\boldsymbol{B}}_{2}})=\|{\boldsymbol{X}}{\boldsymbol{B}}_{1}-{\boldsymbol{X}}{\boldsymbol{B}}_{2}\|_{F}^{2}/({2\sigma^{2}})$. Let $P_{\boldsymbol{0}}$ be $\mathcal{N}(\boldsymbol{0},\sigma^{2}\boldsymbol{I}\otimes{\boldsymbol{I}})$. By the assumption again, for any ${\boldsymbol{B}}:r({\boldsymbol{B}})\leq r$, we have

[TABLE]

where we used $\rho({\boldsymbol{B}}_{1},{\boldsymbol{B}}_{2})\leq(r/2)(J+m-r/2)\leq r(J+m-r)$. Therefore,

[TABLE]

Combining (A.1) and (A.2) and choosing a sufficiently small value of $\gamma$, we can apply Theorem 2.7 of Tsybakov, (2009) to get the desired lower bound.

Case (ii): $(q\wedge J+m-r)r\geq J\log(ep/J)$ and $q<J$. Let $\underline{\kappa}=\sigma_{\min}^{2}({\boldsymbol{X}})$ and $\overline{\kappa}=\sigma_{\max}^{2}({\boldsymbol{X}})$. Note that $\sigma_{\min}({\boldsymbol{X}})$ is the smallest nonzero singular-value of ${\boldsymbol{X}}$ (and so the restricted condition number assumption can hold even when $p>n$). Assume that the SVD of ${\boldsymbol{X}}$ is given by ${\boldsymbol{U}}{\boldsymbol{D}}{\boldsymbol{V}}^{T}$ with ${\boldsymbol{D}}$ of size $q\times q$. Define ${\boldsymbol{Z}}:={\boldsymbol{U}}{\boldsymbol{D}}$ and $\tilde{\mathcal{S}}(r)=\{{\boldsymbol{A}}\in\mathbb{R}^{q\times m}:r({\boldsymbol{A}})\leq r\}$. Then for any estimator $\hat{\boldsymbol{B}}$, and $\hat{\boldsymbol{A}}={\boldsymbol{V}}^{T}\hat{\boldsymbol{B}}$, we have

[TABLE]

because $\tilde{\mathcal{S}}(r)\subset\{{\boldsymbol{V}}^{T}{\boldsymbol{B}}:{\boldsymbol{B}}\in S(J,r)\}$ under $q<J$. The new design matrix ${\boldsymbol{Z}}$ has $q$ columns, and satisfies $\underline{\kappa}\|{\boldsymbol{A}}\|_{F}^{2}\leq\|{\boldsymbol{Z}}{\boldsymbol{A}}\|_{F}^{2}\leq\overline{\kappa}\|{\boldsymbol{A}}\|_{F}^{2}$ for any ${\boldsymbol{A}}$. Repeating the argument in (i) gives the result.

Case (iii): $J\log(ep/J)\geq(J+m-r)r$ and $J\leq q$. Construct

[TABLE]

where $R=\frac{\sigma}{{\overline{\kappa}^{1/2}}}(\frac{\log(ep/J)}{m})^{1/2}$ and $\gamma>0$ is a small constant. Clearly, ${\mathcal{B}}^{3}(J)\subset{S}(J,r)$. By Stirling’s approximation, $\log|\mathcal{B}^{3}(J)|\geq\log{p\choose J}\geq J\log(p/J)\geq cJ\log(ep/J)$ for some universal constant $c$. Due to Lemma A.3 in Rigollet and Tsybakov, (2011), there exists a subset ${\mathcal{B}}^{30}(J)\subset{\mathcal{B}}^{3}(J)$ such that

[TABLE]

for some universal constants $c_{1},c_{2}>0$. The afterward treatment follows the same lines as in (i) and the details are omitted.

Appendix B Proof of Theorem 2

The theorem can be derived from Theorem 2 of She, (2017). For the sake of completeness, we provide its proof, which will be used for proving other theorems. Recall that a random variable $\xi$ is sub-Gaussian if $\,\mathbb{P}(|\xi|\geq t)\leq Ce^{-ct^{2}}$ for any $t>0$ and some constants $C,c>0$, and its scale is defined as $\sigma(\xi)=\inf\{\sigma>0:\,\mathbb{E}\{\exp(\xi^{2}/\sigma^{2})\}\leq 2\}$; ${\boldsymbol{\xi}}\in\mathbb{R}^{p}$ is a sub-Gaussian random vector with its scale bounded by $\sigma$, if $\langle{\boldsymbol{\xi}},{\boldsymbol{\alpha}}\rangle$ is sub-Gaussian and $\sigma(\langle{\boldsymbol{\xi}},{\boldsymbol{\alpha}}\rangle)\leq\sigma\|{\boldsymbol{\alpha}}\|_{2}$ for any ${\boldsymbol{\alpha}}\in\mathbb{R}^{p}$. In this proof, given any index set $\mathcal{J}\subset[p]:=\{1,\ldots,p\}$, ${\boldsymbol{\mathcal{P}}}_{{\boldsymbol{X}}_{\mathcal{J}}}$ is abbreviated to ${\boldsymbol{\mathcal{P}}}_{\mathcal{J}}$ when there is no ambiguity. Given any matrix ${\boldsymbol{A}}$, we use $cs({\boldsymbol{A}})$ and $rs({\boldsymbol{A}})$ to denote its column space and row space, respectively. Because $P_{o}({\boldsymbol{B}})=\sigma^{2}[\{q\wedge J({\boldsymbol{B}})+m-r({\boldsymbol{B}})\}r({\boldsymbol{B}})+J({\boldsymbol{B}})\log\{ep/J({\boldsymbol{B}})\}]$ only depends on $J({\boldsymbol{B}})$ and $r({\boldsymbol{B}})$, we also denote it by $P_{o}(J({\boldsymbol{B}}),r({\boldsymbol{B}}))$. The optimality of $\hat{\boldsymbol{B}}$ implies that for any ${\boldsymbol{B}}={\boldsymbol{B}}_{l}$

[TABLE]

Let ${\boldsymbol{\Delta}}=\hat{\boldsymbol{B}}-{\boldsymbol{B}}$, $\hat{\mathcal{J}}=\mathcal{J}(\hat{\boldsymbol{B}})$, $\mathcal{J}=\mathcal{J}({\boldsymbol{B}})$, $J=J({\boldsymbol{B}})$, $\hat{J}=J(\hat{\boldsymbol{B}})$, $r=r({\boldsymbol{B}})$, $\hat{r}=r(\hat{\boldsymbol{B}})$. Denote the orthogonal projection onto the row space of ${\boldsymbol{X}}_{{\mathcal{J}}}{\boldsymbol{B}}_{{\mathcal{J}}}$ by ${\boldsymbol{\mathcal{P}}}_{rs}$ and its orthogonal complement by ${\boldsymbol{\mathcal{P}}}_{rs}^{\perp}$. Decompose ${\boldsymbol{X}}{\boldsymbol{\Delta}}$ as follows

[TABLE]

Then $\|{\boldsymbol{X}}{\boldsymbol{\Delta}}\|_{F}^{2}=\|{\boldsymbol{\mathcal{P}}}_{\mathcal{J}}{\boldsymbol{X}}{\boldsymbol{\Delta}}{\boldsymbol{\mathcal{P}}}_{rs}\|_{F}^{2}+\|{\boldsymbol{\mathcal{P}}}_{\mathcal{J}}^{\perp}{\boldsymbol{X}}_{\hat{\mathcal{J}}}\hat{\boldsymbol{B}}_{\hat{\mathcal{J}}}{\boldsymbol{\mathcal{P}}}_{rs}\|_{F}^{2}+\|{\boldsymbol{X}}_{\hat{\mathcal{J}}}\hat{\boldsymbol{B}}_{\hat{\mathcal{J}}}{\boldsymbol{\mathcal{P}}}_{rs}^{\perp}\|_{F}^{2}$ and so

[TABLE]

Lemma 1**.**

(She,, 2017, Lemma 4)* Suppose $\operatorname{\mbox{vec}\,}({\boldsymbol{E}})$ is sub-Gaussian with mean zero and scale bounded by $\sigma$. Given ${\boldsymbol{X}}\in\mathbb{R}^{n\times p}$, $1\leq J\leq p$, $1\leq r\leq J\wedge m$, define $\Gamma_{J,r}=\{{\boldsymbol{\Delta}}\in\mathbb{R}^{n\times m}:\|{\boldsymbol{\Delta}}\|_{F}\leq 1,r({\boldsymbol{\Delta}})\leq r,cs({\boldsymbol{\Delta}})\subset cs({\boldsymbol{X}}_{\mathcal{J}})\mbox{ for some }\mathcal{J}:|\mathcal{J}|\leq J\}$. Let*

[TABLE]

Then for any $t\geq 0$,

[TABLE]

where $L,C,c>0$ are universal constants.

Applying the lemma to the first term on the right-hand side of (B.3), we have for any $a,b,a^{\prime}>0$,

[TABLE]

We claim that that $\,\mathbb{E}R^{2}\leq C\sigma^{2}$ when choosing $4b>a$. In fact,

[TABLE]

Similarly, noticing that $r({\boldsymbol{X}}_{\hat{\mathcal{J}}}\hat{\boldsymbol{B}}_{\hat{\mathcal{J}}}{\boldsymbol{\mathcal{P}}}_{rs}^{\perp})\leq r(\hat{\boldsymbol{B}}_{\hat{\mathcal{J}}})\leq\hat{r}$ and $CS({\boldsymbol{X}}_{\hat{\mathcal{J}}}\hat{\boldsymbol{B}}_{\hat{\mathcal{J}}}\allowbreak{\boldsymbol{\mathcal{P}}}_{rs}^{\perp})\allowbreak\subset CS({\boldsymbol{X}}_{\hat{\mathcal{J}}})$, we have $\langle{\boldsymbol{E}},{\boldsymbol{X}}_{\hat{\mathcal{J}}}\hat{\boldsymbol{B}}_{\hat{\mathcal{J}}}{\boldsymbol{\mathcal{P}}}_{rs}^{\perp}\rangle-\|{\boldsymbol{\mathcal{P}}}_{\hat{\mathcal{J}}}{\boldsymbol{X}}{\boldsymbol{\Delta}}{\boldsymbol{\mathcal{P}}}_{rs}^{\perp}\|_{F}^{2}/a-bLP_{o}(\hat{J},\hat{r})\leq\|{\boldsymbol{\mathcal{P}}}_{\mathcal{J}}{\boldsymbol{X}}{\boldsymbol{\Delta}}{\boldsymbol{\mathcal{P}}}_{rs}\|_{F}^{2}/{a^{\prime}}\allowbreak+{a^{\prime}}R^{2}/4$. To handle the second term in (B.3), we introduce the following lemma.

Lemma 2**.**

(She,, 2017, Lemma 5)* Suppose $\operatorname{\mbox{vec}\,}({\boldsymbol{E}})$ is sub-Gaussian with mean zero and $\psi_{2}$-norm bounded by $\sigma$. Given ${\boldsymbol{X}}\in\mathbb{R}^{n\times p}$, $1\leq J\leq p$, $1\leq J^{\prime}\leq p$, $1\leq r\leq J\wedge m$, define $\Gamma_{J^{\prime},J,r}=\{{\boldsymbol{\Delta}}\in\mathbb{R}^{n\times m}:\|{\boldsymbol{\Delta}}\|_{F}\leq 1,r({\boldsymbol{\Delta}})\leq r,cs({\boldsymbol{\Delta}})\subset cs({\boldsymbol{\mathcal{P}}}_{{\mathcal{J}}^{\prime}}^{\perp}{\boldsymbol{\mathcal{P}}}_{{\mathcal{J}}})\mbox{ for some }{\mathcal{J}}^{\prime},{\mathcal{J}}\subset[p]\mbox{ satisfying }|{\mathcal{J}}^{\prime}|\allowbreak\leq J^{\prime},|{\mathcal{J}}|\leq J\}$. Let*

[TABLE]

Then for any $t\geq 0$,

[TABLE]

where $L,C,c>0$ are universal constants.

Because $P_{o}^{\prime\prime}(J,\hat{J},\hat{r})\leq P_{o}(J,r)+P_{o}(\hat{J},\hat{r})$,

[TABLE]

where $\,\mathbb{E}{R^{\prime}}^{2}\leq C\sigma^{2}$ and $L$ is a sufficiently large constant. In summary, we have

[TABLE]

Choose constants $a$, $a^{\prime}$, $b$, $b^{\prime}$, and $A$ sufficiently large such that $({1}/{a}+{1}/{a^{\prime}})(1+{1}/{b^{\prime}})<{1}/{2}$, $4b>a$, and $A>2bL$. Setting ${\boldsymbol{B}}={\boldsymbol{B}}_{l}$ (which is nonzero), and taking expectation on both sides, we obtain the inequality as desired.

Appendix C Proof of Theorem 3

We take (3.8) as an instance to show the conclusion. First, $1/(1-\delta(\hat{\boldsymbol{B}}))^{2}\geq 1/(1-\delta(\hat{\boldsymbol{B}}))$ due to $0\leq\delta(\hat{\boldsymbol{B}})<1$. Let $A_{0}$ and $A$ be chosen such that $2A_{0}>A$ and so $\delta({\boldsymbol{B}}^{*})<1/2$. Then $1/(1-2\delta({\boldsymbol{B}}^{*}))\geq 1/(1-\delta({\boldsymbol{B}}^{*}))^{2}$.

Let $h({\boldsymbol{B}};A)=1/\{mn-AP_{o}({\boldsymbol{B}})\}$. Based on the previous facts and the optimality of $\hat{\boldsymbol{B}}$, we have

[TABLE]

and so $\|{\boldsymbol{Y}}-{\boldsymbol{X}}\hat{\boldsymbol{B}}\|_{F}^{2}\leq\|{\boldsymbol{Y}}-{\boldsymbol{X}}{\boldsymbol{B}}^{*}\|_{F}^{2}\,h({\boldsymbol{B}}^{*};2A)/h(\hat{\boldsymbol{B}};A)$ since $h(\hat{\boldsymbol{B}};A)>0$. It follows that

[TABLE]

The stochastic term $2\langle{\boldsymbol{E}},{\boldsymbol{X}}\hat{\boldsymbol{B}}-{\boldsymbol{X}}{\boldsymbol{B}}^{*}\rangle$ can be decomposed and bounded in the same way as in the proof of Theorem 2, except that we use the high-probability form results here. For example, for term $\langle{\boldsymbol{E}},{\boldsymbol{\mathcal{P}}}_{\mathcal{J}}{\boldsymbol{X}}{\boldsymbol{\Delta}}{\boldsymbol{\mathcal{P}}}_{rs}\rangle$ in (B.3), Lemma 1 shows that for any constants $a,b,a^{\prime}>0$ satisfying $4b>a$, the following event

[TABLE]

occurs with probability at least $1-\sum_{J=1}^{p}\sum_{r=1}^{m\wedge J}C\exp[-c\{(2({{b}/{a}})^{1/2}-1\}^{2}L\allowbreak P_{o}(J,r)/\sigma^{2}]$ or $1-Cp^{-c}$ for a sufficiently large value of $L$. Repeating the analysis in the proof of Theorem 2 shows

[TABLE]

with probability at least $1-Cp^{-c}$ for some $c,C>0$.

Let $\gamma$ and $\gamma^{\prime}$ be constants satisfying $0<\gamma<\ 1,\gamma^{\prime}>0$. On $\mathcal{E}=\{(1-\gamma){mn\sigma^{2}}\leq\|{\boldsymbol{E}}\|_{F}^{2}\leq(1+\gamma^{\prime}){mn\sigma^{2}}\}$ , we have

[TABLE]

From Laurent and Massart, (2000), the complement of $\mathcal{E}$ occurs with probability at most $C^{\prime}\exp(-c^{\prime}mn)$ (with $c^{\prime},C^{\prime}$ dependent on $\gamma,\gamma^{\prime}$). With $A_{0}$ large enough, we can choose $a,a^{\prime},b,A$ such that $({1}/{a}+{1}/{a^{\prime}})<{1}/{2}$, $4b>a$, and $4bL\leq(1-\gamma)A$. The conclusion results, and it is easy to see that it holds for the fractional form as well.

Noticing that $1/(1-\delta)\geq\exp(\delta)\geq 1+\delta\geq 1/(1-\delta/2)$ for any $0\leq\delta<1$, we can show the same conclusion for the log-form and the plug-in form, the proof of which follows the same lines.

Appendix D Proofs of Theorem 4 and Corollary 1

For notational clarity, we prove the conclusion for $K=2$ only. The proof in the general case follows the same lines. Let $\bar{r}=r$ and ${\boldsymbol{O}}={\boldsymbol{S}}{\boldsymbol{U}}$ (which has $r$ orthogonal columns). Denote the overall data by $(\bar{\boldsymbol{X}},\bar{\boldsymbol{Y}})$. Let $({\boldsymbol{X}},{\boldsymbol{Y}})$ be the training data and $(\tilde{\boldsymbol{Y}},\tilde{\boldsymbol{X}})$ be the test data in the first round. Then $\hat{\boldsymbol{B}}={\boldsymbol{O}}\hat{\boldsymbol{C}}$ and $\hat{\boldsymbol{C}}=({\boldsymbol{X}}_{{\boldsymbol{O}}}^{T}{\boldsymbol{X}}_{{\boldsymbol{O}}})^{-1}{\boldsymbol{X}}_{{\boldsymbol{O}}}^{T}{\boldsymbol{Y}}$ (where ${\boldsymbol{X}}_{{\boldsymbol{O}}}$ is short for ${\boldsymbol{X}}{\boldsymbol{O}}$ throughout the proof). Also, from the model assumption, we can write ${\boldsymbol{B}}^{*}={\boldsymbol{O}}^{*}{\boldsymbol{C}}^{*}$ and ${\boldsymbol{X}}{\boldsymbol{B}}^{*}={\boldsymbol{X}}_{{\boldsymbol{O}}^{*}}{\boldsymbol{C}}^{*}$ for some ${\boldsymbol{C}}^{*}$, where ${\boldsymbol{O}}^{*}$ is the $s$-pattern of ${\boldsymbol{B}}^{*}$.

First, we calculate the extra-sample test error on $(\tilde{\boldsymbol{Y}},\tilde{\boldsymbol{X}})$:

[TABLE]

where we used $\,\mathbb{E}[\bar{\boldsymbol{E}}\bar{\boldsymbol{E}}^{T}]=m\sigma^{2}{\boldsymbol{I}}$ and the independence between $\bar{\boldsymbol{E}}$ and $\bar{\boldsymbol{X}}$. The second term can be simplified further:

[TABLE]

Therefore, the total cross-validation error is given by

[TABLE]

Similarly, the overall training error is

[TABLE]

It remains to analyze $\,\mathbb{E}[\|\tilde{\boldsymbol{X}}{\boldsymbol{B}}^{*}-\tilde{\boldsymbol{X}}{\boldsymbol{O}}({\boldsymbol{X}}_{\boldsymbol{O}}^{T}{\boldsymbol{X}}_{\boldsymbol{O}})^{-1}{\boldsymbol{X}}_{\boldsymbol{O}}^{T}{\boldsymbol{X}}{\boldsymbol{B}}^{*}\|_{F}^{2}]]+\,\mathbb{E}[\|{\boldsymbol{X}}{\boldsymbol{B}}^{*}-\allowbreak{\boldsymbol{X}}{\boldsymbol{O}}(\tilde{\boldsymbol{X}}_{\boldsymbol{O}}^{T}\tilde{\boldsymbol{X}}_{\boldsymbol{O}})^{-1}\allowbreak\tilde{\boldsymbol{X}}_{\boldsymbol{O}}^{T}\tilde{\boldsymbol{X}}{\boldsymbol{B}}^{*}\|_{F}^{2}]]-\,\mathbb{E}[\|\bar{\boldsymbol{X}}{\boldsymbol{B}}^{*}-\bar{\boldsymbol{X}}_{\boldsymbol{O}}(\bar{\boldsymbol{X}}_{\boldsymbol{O}}^{T}\bar{\boldsymbol{X}}_{\boldsymbol{O}})^{-1}\bar{\boldsymbol{X}}_{\boldsymbol{O}}^{T}\bar{\boldsymbol{X}}{\boldsymbol{B}}^{*}\|_{F}^{2}]].$

Define ${\boldsymbol{\mathcal{P}}}_{1}={\boldsymbol{\mathcal{P}}}_{{\boldsymbol{O}}}\cap{\boldsymbol{\mathcal{P}}}_{{\boldsymbol{O}}^{*}}$, ${\boldsymbol{\mathcal{P}}}_{2}={\boldsymbol{\mathcal{P}}}_{1}^{\perp}\cap{\boldsymbol{\mathcal{P}}}_{{\boldsymbol{O}}^{*}}(={\boldsymbol{\mathcal{P}}}_{{\boldsymbol{O}}^{*}}\setminus({\boldsymbol{\mathcal{P}}}_{{\boldsymbol{O}}}\cap{\boldsymbol{\mathcal{P}}}_{{\boldsymbol{O}}^{*}}))=\mathfrak{O}\mathfrak{O}^{T}$. Then ${\boldsymbol{B}}^{*}={\boldsymbol{O}}^{*}{\boldsymbol{C}}^{*}=({\boldsymbol{\mathcal{P}}}_{1}+{\boldsymbol{\mathcal{P}}}_{2}){\boldsymbol{B}}^{*}$, and so

[TABLE]

With this fact, and using the shorthand notation ${\boldsymbol{X}}_{1}={\boldsymbol{X}}{\boldsymbol{O}}$, ${\boldsymbol{X}}_{2}={\boldsymbol{X}}{\mathfrak{O}}$, ${\boldsymbol{B}}_{2}^{*}={\mathfrak{O}}^{T}{\boldsymbol{B}}^{*}$, we get

[TABLE]

and

[TABLE]

Introduce ${\boldsymbol{\Omega}}_{1}={\boldsymbol{X}}_{1}^{T}{\boldsymbol{X}}_{1}$, ${\boldsymbol{\Omega}}_{2}={\boldsymbol{X}}_{2}^{T}{\boldsymbol{X}}_{2}$, ${\boldsymbol{\Omega}}_{12}={\boldsymbol{X}}_{1}^{T}{\boldsymbol{X}}_{2}$. In the following, we investigate

[TABLE]

Define ${\boldsymbol{\Omega}}_{1}^{-1}{\boldsymbol{\Omega}}_{12}=({\boldsymbol{X}}_{1}^{T}{\boldsymbol{X}}_{1})^{-1}{\boldsymbol{X}}_{1}^{T}{\boldsymbol{X}}_{2}={\boldsymbol{R}}+{\boldsymbol{\Delta}}$, where  ${\boldsymbol{R}}=\,\mathbb{E}[({\boldsymbol{X}}_{1}^{T}{\boldsymbol{X}}_{1})^{-1}{\boldsymbol{X}}_{1}^{T}{\boldsymbol{X}}_{2}]$ (non-random) and $\,\mathbb{E}{\boldsymbol{\Delta}}=\boldsymbol{0}$. By symmetry, $\tilde{\boldsymbol{\Omega}}_{1}^{-1}\tilde{\boldsymbol{\Omega}}_{12}={\boldsymbol{R}}+\tilde{\boldsymbol{\Delta}}$ with $\,\mathbb{E}\tilde{\boldsymbol{\Delta}}=\boldsymbol{0}$. Represent $\bar{\boldsymbol{\Omega}}_{1}^{-1}\bar{\boldsymbol{\Omega}}_{12}$ by ${\boldsymbol{R}}+\bar{\boldsymbol{\Delta}}$ (but $\,\mathbb{E}\bar{\boldsymbol{\Delta}}$ may not be $\boldsymbol{0}$ in general).

Writing

[TABLE]

and using the additivity of

[TABLE]

(D.1) can be reduced to the sum of

[TABLE]

and

[TABLE]

In expectation, $2({\boldsymbol{R}}^{T}\tilde{\boldsymbol{\Omega}}_{1}-\tilde{\boldsymbol{\Omega}}_{21}){\boldsymbol{\Delta}}$ and $2({\boldsymbol{R}}^{T}{\boldsymbol{\Omega}}_{1}-{\boldsymbol{\Omega}}_{21})\tilde{\boldsymbol{\Delta}}$ vanish. From the data splitting manner and the row independence, we have

[TABLE]

where ${\boldsymbol{\Sigma}}_{1}=\,\mathbb{E}\bar{\boldsymbol{\Omega}}_{1}={\boldsymbol{O}}^{T}\,\mathbb{E}[\bar{\boldsymbol{X}}^{T}\bar{\boldsymbol{X}}]{\boldsymbol{O}}={\boldsymbol{O}}^{T}(n{\boldsymbol{\Sigma}}){\boldsymbol{O}}$. Finally,

[TABLE]

The conclusion follows.

To show the result in the corollary, use the fact that under the Gaussian design assumption, ${\boldsymbol{Z}}_{n-d}^{T}{\boldsymbol{Z}}_{n-d}$ follows a Wishart distribution $\mathcal{W}_{r}({\boldsymbol{I}},n-d)$, and so $({\boldsymbol{Z}}_{n-d}^{T}{\boldsymbol{Z}}_{n-d})^{-1}$ follows an inverse Wishart distribution and has mean ${\boldsymbol{I}}_{r}/(n-d-1)$.

Appendix E Rowwise error bound for support recovery

Given ${\boldsymbol{B}}=[{\boldsymbol{\beta}}_{1},\ldots,{\boldsymbol{\beta}}_{p}]^{T}\in\mathbb{R}^{p\times m}$, define $\|{\boldsymbol{B}}\|_{2,\infty}=\max_{1\leq j\leq p}\|{\boldsymbol{\beta}}_{j}\|_{2}$. Recall $J^{*}=J({\boldsymbol{B}}^{*})=\|{\boldsymbol{B}}^{*}\|_{2,0}$ and $r^{*}=r({\boldsymbol{B}}^{*})$.

Theorem 5**.**

Assume ${\boldsymbol{Y}}={\boldsymbol{X}}{\boldsymbol{B}}^{*}+{\boldsymbol{E}}$ where $\operatorname{\mbox{vec}\,}({\boldsymbol{E}})$ is sub-Gaussian with mean $\boldsymbol{0}$ and scale bounded by $\sigma$ and ${\boldsymbol{B}}^{*}\neq\boldsymbol{0}$. Let $\hat{\boldsymbol{B}}={\mathop{\rm arg\,min}}_{{\boldsymbol{B}}}\frac{1}{2}\|{\boldsymbol{Y}}-{\boldsymbol{X}}{\boldsymbol{B}}\|^{2}_{F}+P({\boldsymbol{B}};A)$ where $P({\boldsymbol{B}};A)=A\sigma^{2}P_{o}({\boldsymbol{B}})$ and $A$ is a constant. Assume that there exist $\kappa>0$ and large enough $A$ such that $\kappa nJ^{*}\|{\boldsymbol{B}}-{\boldsymbol{B}}^{*}\|_{2,\infty}^{2}/2\leq\|{\boldsymbol{X}}{\boldsymbol{B}}-{\boldsymbol{X}}{\boldsymbol{B}}^{*}\|_{F}^{2}/2+P({\boldsymbol{B}};A)+P({\boldsymbol{B}}^{*};A)$ for any ${\boldsymbol{B}}\in\mathbb{R}^{p\times m}$. Let ${\boldsymbol{B}}^{*}=[{\boldsymbol{\beta}}_{1}^{*},\ldots,{\boldsymbol{\beta}}_{p}^{*}]^{T}$, $\hat{\boldsymbol{B}}=[\hat{\boldsymbol{\beta}}_{1},\ldots,\hat{\boldsymbol{\beta}}_{p}]^{T}$ and $\hat{\mathcal{J}}=\mathcal{J}(\hat{\boldsymbol{B}})$. Then

[TABLE]

If, in addition, the minimum signal strength satisfies

[TABLE]

then with probability at least $1-Cp^{-c}$ for some constants $C,c>0$, ${\mathcal{J}}^{*}\subset\hat{\mathcal{J}}$ and $\|\hat{\boldsymbol{\beta}}_{j}\|_{2}$ for all $j\in\mathcal{J}^{*}$ are exactly the $J^{*}$ largest norms in $\|\hat{\boldsymbol{\beta}}_{j}\|_{2}$ ($1\leq j\leq p$).

Proof.

First, the stochastic term $2\langle{\boldsymbol{E}},{\boldsymbol{X}}\hat{\boldsymbol{B}}-{\boldsymbol{X}}{\boldsymbol{B}}^{*}\rangle$ can be decomposed and bounded in the same way as in the proof of Theorem 2, except that we use the high-probability form results here. For example, for term $\langle{\boldsymbol{E}},{\boldsymbol{\mathcal{P}}}_{\mathcal{J}}{\boldsymbol{X}}{\boldsymbol{\Delta}}{\boldsymbol{\mathcal{P}}}_{rs}\rangle$ in (B.3), Lemma 1 shows that for any constants $a,b,a^{\prime}>0$ satisfying $4b>a$, the following event

[TABLE]

occurs with probability at least $1-\sum_{J=1}^{p}\sum_{r=1}^{m\wedge J}C\exp[-c\{(2({{b}/{a}})^{1/2}-1\}^{2}L\allowbreak P_{o}(J,r)/\sigma^{2}]$ or $1-Cp^{-c}$ for a sufficiently large value of $L$. Repeating the analysis in the proof of Theorem 2 shows

[TABLE]

with probability at least $1-Cp^{-c}$ for some $c,C>0$, from which it follows that

[TABLE]

From the assumption,

[TABLE]

Hence we obtain

[TABLE]

Let $1/a+1/a^{\prime}=1/(2(1+\alpha))$, $4b>a$ and $A={2bL}(1+\alpha)$ with $\alpha>0$. Then

[TABLE]

Taking the optimal $\alpha=1$, we get (E.1).

Furthermore, under (E.2), for any $j\in{\mathcal{J}}^{*}$

[TABLE]

Combining it with $\|\hat{\boldsymbol{\beta}}_{j^{\prime}}\|_{2}^{2}\leq 8P({\boldsymbol{B}}^{*};A)/(\kappa nJ^{*})$ for any $j^{\prime}\in{{\mathcal{J}}^{*}}^{c}$ gives the conclusion. ∎