# An Orthogonal Matching Pursuit Variable Screening Algorithm for High-Dimensional Linear Regression Models

Authors: Yanxi Xie; Yuewen Li; Victor Shi; Quan Lu
Journal: Scientific Programming
 (2022)

---

## Abstract

Variable selection plays an important role in data mining. It is crucial to filter useful variables and extract useful information in a high-dimensional setup when the number of predictor variablesd tends to be much larger than the sample size n. Statistical inferences can be more precise after irrelevant variables are moved out by the screening method. This article proposes an orthogonal matching pursuit algorithm for variable screening under the high-dimensional setup. The proposed orthogonal matching pursuit method demonstrates good performance in variable screening. In particular, if the dimension of the true model is finite, OMP might discover all relevant predictors within a finite number of steps. Throughout theoretical analysis and simulations, it is confirmed that the orthogonal matching pursuit algorithm can identify relevant predictors to ensure screening consistency in variable selection. Given the sure screening property, the BIC criterion can be used to practically select the best candidate from the models generated by the OMP algorithm. Compared with the traditional orthogonal matching pursuit method, the resulting model can improve prediction accuracy and reduce computational cost by screening out the relevant variables.


---

## Body

## 1. Introduction

Variable screening is an important technique in data mining. It captures informative variables by reducing the dimension in a high-dimensional setup when the number of predictor variablesd tends to be much larger than the sample size n. However, statistical inference is difficult to compute in ultra-high dimensional linear models before variable screening due to the computational complexity. It is necessary to remove the irrelevant variables from the model before statistical inference. The core idea is to screening out the informative variables with the aim of building a relevant model for future prediction. By removing most irrelevant and redundant variables from the data, variable selection helps improve the performance of learning models in terms of obtaining higher estimation accuracy [1]. Then the AIC [2] or BIC [3] can be applied to further guarantee the accuracy of the relevant model.The focus of this article is on ultra-high dimensional linear models, in which the number of predictor variablesd tends to be much larger than the sample size n. In particular, the number of covariates may increase at an exponential rate. Such linear models have gained a lot of attention in practical areas, such as sentiment analysis and finance. Existing techniques in the past literature include forward selection [1], least absolute shrinkage and selection operator (Lasso) [4], smoothly clipped absolute deviation penalty (SCAD) [5], etc. These efforts have been devoted to the challenging ultra-high dimensionality problem, which is motivated by contemporary applications such as bioinformatics, genomics, finance, etc. In other words, it is becoming a major issue to investigate the existence of complex relationships and dependencies in data with the aim of building a relevant model for inference. A practically attractive approach is to first use a quick screening procedure to reduce the dimensionality of the covariates to a reasonable scale, for example below the sample size, and then apply variable selection techniques such as LASSO and SCAD in the second stage.Motivated by the current studies on variable screening approaches in ultra-high dimensional linear models, it is of our interests in showing the screening consistency property of the OMP under certain conditions, by restricting the technical conditions stated in Wang [6] and hence selecting a subset of predictors, which includes all relevant predictors, to ensure variable screening results.The rest of this article is organized as follows: Section2 provides the literature review on current variable screening methods. Section 3 demonstrates a variable screening algorithm based on the OMP. Furthermore, the asymptotic results of the estimators are studied. Section 4 examines via simulation that our proposed technique exhibits desired sample properties and can be useful in practical applications. Finally, Section 5 concludes the article and provides some future research direction. In particular, the proof of the asymptotic theories and lemmas can be found in Appendix.


## 2. Literature Review

In the content of variable selection, screening approaches have gained a lot of attention besides the penalty approaches such as Lasso [4] and SCAD [5]. When the predictor dimension is much larger than the sample size, the story changes drastically in the sense that the conditions for most of the Lasso-type algorithms cannot be satisfied. Therefore, to conduct model selection in the high-dimensional setup, variable screening is a reasonable solution.Sure independence screening (SIS), which is proposed by Fan and Lv [7], has gained popularity under the condition when the number of predictor variables d tends to be much larger than the sample size n. Sure screening means a property that all the important variables are selected after applying a variable screening procedure with probability tending to 1. It is desired to have a dimensionality reduction method with the sure screening property. There are three facts why sure screening is of great importance and usage when dimension d is larger than sample size n, which is clearly stated in Fan and Lv [7]. First of all, the design matrix X is rectangular, having more columns than rows. In this case, the matrix XTX is giant in dimension and singular. The maximum spurious correlation between a covariate and a response can be large due to the dimensionality and the fact that an unimportant predictor can be highly correlated with the response variable owing to the presence of important predictors associated with the predictor. In addition, the population covariance matrix Σ may become ill conditioned as n grows, and it makes variable selection difficult. Third, the minimum nonzero absolute coefficient βj may decay with n and fall close to the noise level, say, the order logd/n−1/2. Hence, in general, it becomes challenging to estimate the sparse parameter vector β accurately when d≫n.To solve the abovementioned difficulties in variable selection, Fan and Lv [7] proposed a simple sure screening method using componentwise regression or equivalently correlation learning, to reduce dimensionality from high to moderate scale that is below sample size. Below is the description of the SIS method.Letω=ω1,…,ωdT be a d-vector that is obtained by componentwise regression, that is(1)ω=XTy,where the n×d data matrix X is first standardized columnwise. For any given γ∈0,1, we sort the d componentwise magnitudes of the vector ω in a descending order and define a submodel(2)Mγ=1≤i≤d:ωiis among the firstγnlargest of all,where γn denotes the integer part of γn. It shrinks the full model 1,2,…,d down to a submodel Mγ with size γn smaller than the sample size n. This correlation learning ranks the importance of features according to their marginal correlations with the response variable. Moreover, it is called the independence screening because each feature is used independently as a predictor to decide the usefulness for predicting the response variable. The computational cost of SIS is of order Ond.With dimension reduced accurately from high to below sample size, variable selection can be improved on both speed and accuracy, and can then be accomplished by a well-developed method such as SCAD, Lasso, or adaptive Lasso [8, 9], denoted by SIS-SCAD, SIS-Lasso, or SIS-AdapLasso, respectively. Moreover, sure screening property has been proven in Fan and Lv [7]. Intuitively, the core idea of SIS is to select the variables by two stages. In the first stage, an easy-to-implement method is used to remove the least important variables. In the second stage, a more sophisticated and accurate method is applied to reduce the variables further.Though SIS enjoys sure screening property and is easy to be applied, it has several potential problems. First of all, if there is an important predictor jointly correlated but marginally uncorrelated with the response variable, it is not selected by SIS and thus cannot be included in the estimated model. Second, similar to Lasso, SIS cannot handle the collinearity problem between predictors in terms of variable selection. Third, when there are some unimportant predictors which are highly correlated with the important predictors, these unimportant predictors can have higher chance of being selected by SIS than other important predictors that are relatively weakly related to the response variable. In all, these three potential issues can be carefully treated when some extensions of SIS are proposed. In particular, iterative SIS (ISIS) is designed to overcome the weakness of SIS.ISIS works in two steps. In the first step, a subset ofk1 variables A1=Xi1,…,Xik1 is selected by using an SIS-based model selection method such as SIS-SCAD or SIS-Lasso methods. There is an n-vector of residuals from regressing the response y over Xi1,…,Xik1. In the second step, the residuals are treated as the new response variable and the previous step is repeated to the remaining d−k1 variables. It returns a subset of k2 variables A2=Xj1,…,Xjk2. Fitting the residuals from the previous step on X1,…,XdA1 can significantly weaken the prior selection of those unimportant variables that are highly correlated with the response through their relations with Xi1,…,Xik1. In addition, the second step also makes those important variables which are missed out in the first step possible to be selected. Iteratively, the second step is iterated until l disjoint subsets A1,…,Al are obtained with the union A=∪j=1lAj has a size γn.If SIS is used to select only one variable at each iteration, that isAi=1, ISIS is equivalent to orthogonal matching pursuit (OMP) [10], which is a greedy algorithm for variable selection. This is discussed in the study by Barron and Cohen [11].Kim [12] proposed a filter ranking method using the elastic net penalty with sure independence screening (SIS) on resampling technique to overcome the overfitting and high-performance computational issues. It is demonstrated via extensive simulation studies that SIS-LASSO, SIS-MCP, and SIS-SCAD with the proposed filtering method achieve superior performance of not only accuracy but also true positive detection compared to those with the marginal maximum likelihood ranking (MMLR) method.Another very popular yet classical variable screening method is the forward regression (FR). As one type of important greedy algorithms, FR’s theoretical properties have been investigated in the studies by Barron and Cohen [11], Donoho and Stodden [13], and Wang [6]. In particular, Wang [6] investigated FR’s screening consistency property, under an ultra-high dimensional setup, by introducing the four technical conditions.There are a few comments on those four technical conditions introduced in the study by Wang [6]. First of all, the normality assumption has been popularly used in the past literature for theory development. Second, the smallest and largest eigenvalues of the covariance matrix Σ need to be properly bounded. This bounded condition together with the normality assumption ensures the sparse Riesz condition (SRC) defined in the study by Zhang and Huang [14]. Third, the standard L2 norm of the regression coefficients β is bounded above by some proper constant. It guarantees that the signal-to-noise ratio is convergent. Moreover, the minimum value of the nonzero βjs needs to be bounded below. This constraint on the minimal size of the nonzero regression coefficient ensures that relevant predictors can be correctly selected. Otherwise, if some of the nonzero coefficients converge too fast, they cannot be selected consistently. Last but not least, logd is bounded above in the order of nξ for some small constant ξ. This condition allows the predictor dimension d to diverge to infinity at an exponential fast speed, which implies that the predictor dimension can be substantially larger than the sample size n.Under the assumption that the true modelT exists, Wang [6] introduces the FR algorithm in the aim of discovering all relevant predictors consistently. The main step of FR algorithm is the iterative forward regression part. Consider the case where k−1 relevant predictors have been selected accordingly. Then the next step is to construct a candidate model that include one more predictor that belongs to the full set but excluding the selected k−1 predictors and calculate the residual sum of squares based on the constructed candidate model. This step is repeated for each predictor that belongs to the full set but excluding the selected k−1 predictors, and all the residual sum of squares are recorded accordingly. The minimum value of all the recorded residual sum of squares are found, and the kth relevant predictor is updated based on the index of the corresponding minimum residual sum of squares. A detailed algorithm can be found in the study by Wang [6].Wang [6] showed the theoretical proof that FR can identify all relevant predictors consistently, even if the predictor dimension is considerably larger than the sample size. In particular, if the dimension of the true model is finite, FR might discover all relevant predictors within a finite number of steps. In other words, sure screening property can be guaranteed under the four technical conditions. Given the sure screening property, the recently proposed BIC of Chen and Chen [3] can be used to practically select the best candidate from the models generated by the FR algorithm. The resulting model is good in the sense that many existing variable selection methods, such as Adaptive Lasso and SCAD, can be applied directly to increase the estimation accuracy.The extended Bayes information criterion (EBIC) proposed by Chen and Chen [3] is suitable for large model spaces. It has the following form:(3)BICM=logσ^M2+n−1Mlogn+2logd,where M is an arbitrary candidate model with M≤n, σ^M2=n−1RSSM=n−1YTIn−HMY, and HM=XMXMTXM−1XMT. We then select the best model S^=Sm^, where m^=argmin1≤m≤nBICSm.EBIC, which includes the original BIC as a special case, examines both the number of unknown parameters and the complexity of the model space. The model in Chen and Chen [3] is defined to be identifiable if no model of comparable size other than the true submodel can predict the response almost equally well. It has been shown that EBIC is selection consistent under some mild conditions. It also handles the heavy collinearity problem for the covariates. Furthermore, EBIC is easy to implement due to the fact the extended BIC family does not require a data adaptive tuning parameter procedure.Other screening approaches include tournament screening (TS) [15], sequential Lasso [16], quantile-adaptive model-free variable screening [17], and conditional screening [18]. When P≫n, the tournament screening possesses the sure screening property to reduce spurious correlation. Furthermore, the asymptotic properties of sequential Lasso for feature selection in linear regression models with ultra-high dimensional feature spaces are investigated. The advantage of sequential Lasso is that it is not restricted by the dimensionality of the feature space. Quantile-adaptive model-free variable screening has two distinctive features, allowing the set of active variables to vary across quantiles and overcoming the difficulty in specifying the form of a statistical model in a high-dimensional space. Baranowski [19] proposed a workflow representation for scheduling, provenance, or visualization to resolve variable and method dependencies and evaluated the performance of screening properties. Samudrala [20] proposed a parallel algorithm by identifying key components for dimensionality reduction of large-scale data. It shows better performance for dimension reduction compared to the existing methods. Chen [21] proposed a model-free feature screening method when the censored response and error-prone covariates both exist. An iterative algorithm is developed in the presence of the censored response and error-prone covariates. In addition, we also develop the iteration method to improve the accuracy of selecting all important covariates. Choudalakis [22] proposed appropriate numerical methods for parameter estimation under the high-dimensional setup. A thorough comparison is considered among existing methods for both coefficient estimations and variable selection for supersaturated designs. Xu et al. [23–25] proposed several multi-objective robust optimization models for MDVRPLS in refined oil distribution. Ren [26] proposed an asymmetric learning to hash with variable bit encoding algorithm (AVBH) to solve the high-dimensional data problem, and a real data application is applied for the finite performance of the proposed AVBH algorithm. We proposed a parallel framework for dimensionality reduction of large-scale data. We also identified key components underlying the spectral dimensionality reduction techniques and proposed their efficient parallel implementation.


## 3. Main Results

Orthogonal matching pursuit (OMP) is an iterative greedy algorithm that selects at each step the column which is most correlated with the current residuals. The selected column is then added into the set of selected columns. Inspired by the idea of the FR algorithm in Wang (2009), it is shown that under some proper conditions, OMP can enjoy the sure screening property in the linear model setup.


### 3.1. Model Setup and Technical Conditions

LetXi,yi be the observation collected from the ith subject 1≤i≤n, where yi∈R1 is the response and Xi=Xi1,…,XidT∈Rd is the high-dimensional predictor with d>n and covXi=Σ. Moreover, β=β1,…,βdT is the regression coefficient. In matrix representation, the design matrix is X∈Rn×d and the response vector is y∈Rn. Consider the linear regression model as(4)y=Xβ+ε.Without loss of generality, it is assumed that the data are centered, that is the columns ofX are orthonormal and yis are conditionally independent given the design matrix X. Equivalently, EXij=0 and VarXij=1. Moreover, the error term ϵ are independently and identically distributed with mean zero and finite variance σ2. A model fitting procedure produces the vector of coefficients β^=β^1,…,β^dT.Before the main result for the screening property of OMP is presented, four technical conditions are needed as follows:


#### 3.1.1. Assumption 1. Technical Conditions

(C1) Normality assumption. Assume thatε follows the normal distribution.(C2) Covariance matrix.λminA and λmaxA represent, respectively, the smallest and largest eigenvalues of an arbitrary positive definite matrix A. We assume that there exist a positive constant 1<c<∞, such that c−1<λminΣ≤λmaxΣ<c.(C3) Regression coefficients. We assume thatβ≤Cβ for some constant Cβ>0, where . denotes the standard L2 norm. Also assume βmin≥νβn−ξmin for some ξmin>0, with βmin=minj∈Tβj.(C4) Divergence speed ofd and d0. We assume logd=Onξ and d0=Onξ0. In other words, there exists constants ξ, ξmin, and ν, such that logd≤νnξ, d0≤νnξ0, and ξ+6ξ0+12ξmin<1.


### 3.2. OMP Algorithm

Under the assumption that the true setT exists, our main objective is to discover all relevant predictors consistently. To this end, we consider the following OMP algorithm (Algorithm1):Algorithm 1: The OMP algorithm.
Step 1 (Initialization). SetS0=∅. Set the residual r0=y.Step 2 (Forward selection).(i) (2.1) Evaluation. In thekth step k≥1, we are given Sk−1. Then, for every j∈F∖Sk−1, we computeXjTrk−1,whererk−1=In−HSk−1y, HSk−1=XSk−1XSk−1TXSk−1−1XSk−1T is the projection onto the linear space spanned by the elements of XSk−1 and In∈Rn×n is the identity matrix.(ii) (2.2) Screening. We then findak=argmaxj∈F∖Sk−1XjTrk−1,and updateSk=Sk−1∪ak accordingly.Step 3 (Solution path). Iterating Step 2 forn times leads a total of n nested candidate models. We then collect those models by a solution path S=Sk:1≤k≤n with Sk=a1,…,ak.


### 3.3. Theoretical Properties

To prove Theorem1, the following lemmas are needed. For convenience, we define Σ^=n−1XTX. Moreover, Σ^M and ΣM are the submatrices of Σ^ and Σ (corresponding to M), respectively.Lemma 1.
If we have1/c<λminΣ<λmaxΣ<c, then eigenvalues of the submatrix of Σ, that is Σ11, are also bounded. Suppose Σ11 contains m variables, that is Σ11 is a m×m matrix. Moreover, with probability tending to 1, we have(5)1/c<λminΣ^11<λmaxΣ^11<c,as long as m=On1−ξ/3.The Proof of Lemma 1 can be found in Appendix.Lemma 2.
Assume conditions (C1), (C2), (C4), andm=On1−ξ/3. Then, with probability tending to 1, we have(6)τmin≤minM≤mλminΣ^M≤maxM≤mλmaxΣ^M≤τmax.Note that proof of Lemma2 can be found in the study by Wang (2009), with only slight changes.Before the theorems are established, we follow Wang (2009)’s idea on screening consistency of a solution path and define the solution pathS to be screening consistent, if(7)PT⊂Sk∈Sfor some1≤k≤n⟶1.Then OMP’s screening consistency can be formally established by the following theorem.Theorem 1.
Under model (4) and conditions (C1)–(C4), we have as n⟶∞(8)PT⊂SKnξ0+4ξmin⟶1,where the constant K=2τmaxτmin−2Cβ2νβ−4ν is independent of n, the constants τmax, τmin, Cβ, νβ, and ν are defined in conditions (C2)–(C4), and t is the smallest integer no less than t.Theorem1 proves that within Knξ0+4ξmin steps, all relevant predictors will be identified by the OMP algorithm. This number of steps is much smaller than the sample size n under condition (C4). In particular, if the dimension of the true model is finite with ξ0=ξmin=0, only a finite number of steps are needed to discover the entire relevant variable set.Furthermore, Theorem1 provides a theoretical basis for OMP, which enables us to empirically select the best model from S. On the other hand, the solution path S contains a total of n nested models. To further select relevant variables from the solution path S, the following BIC (Chen and Chen, 2008) is considered,(9)BICM=logσ^M2+n−1Mlogn+2logd,where M is an arbitrary candidate model with M≤n, σ^M2=n−1RSSM=n−1YTIn−HMY, and HM=XMXMTXM−1XMT. We then select the best model S^=Sm^, where m^=argmin1≤m≤nBICSm. We typically do not expect S^ to be selection consistent (i.e., PS^=T⟶1). However, we are able to show that S^ is indeed screening consistent.Theorem 2.
Under model (4) and conditions (C1)–(C4), then as n⟶∞(10)PT⊂S^⟶1.Definekmin=min1≤k≤nk:T⊂Sk. By Theorem 2.1, we know that kmin satisfies kmin≤Knξ0+4ξmin with probability tending to 1. Therefore, our aim is to prove that Pm^≤kmin⟶0 as n⟶∞. Then the theorem conclusion follows. Equivalently, it suffices to show that(11)Pmin1≤k<kminBICSk−BICSk+1>0⟶1.A detailed proof can be found in Appendix.


## 3.1. Model Setup and Technical Conditions

LetXi,yi be the observation collected from the ith subject 1≤i≤n, where yi∈R1 is the response and Xi=Xi1,…,XidT∈Rd is the high-dimensional predictor with d>n and covXi=Σ. Moreover, β=β1,…,βdT is the regression coefficient. In matrix representation, the design matrix is X∈Rn×d and the response vector is y∈Rn. Consider the linear regression model as(4)y=Xβ+ε.Without loss of generality, it is assumed that the data are centered, that is the columns ofX are orthonormal and yis are conditionally independent given the design matrix X. Equivalently, EXij=0 and VarXij=1. Moreover, the error term ϵ are independently and identically distributed with mean zero and finite variance σ2. A model fitting procedure produces the vector of coefficients β^=β^1,…,β^dT.Before the main result for the screening property of OMP is presented, four technical conditions are needed as follows:


### 3.1.1. Assumption 1. Technical Conditions

(C1) Normality assumption. Assume thatε follows the normal distribution.(C2) Covariance matrix.λminA and λmaxA represent, respectively, the smallest and largest eigenvalues of an arbitrary positive definite matrix A. We assume that there exist a positive constant 1<c<∞, such that c−1<λminΣ≤λmaxΣ<c.(C3) Regression coefficients. We assume thatβ≤Cβ for some constant Cβ>0, where . denotes the standard L2 norm. Also assume βmin≥νβn−ξmin for some ξmin>0, with βmin=minj∈Tβj.(C4) Divergence speed ofd and d0. We assume logd=Onξ and d0=Onξ0. In other words, there exists constants ξ, ξmin, and ν, such that logd≤νnξ, d0≤νnξ0, and ξ+6ξ0+12ξmin<1.


## 3.1.1. Assumption 1. Technical Conditions

(C1) Normality assumption. Assume thatε follows the normal distribution.(C2) Covariance matrix.λminA and λmaxA represent, respectively, the smallest and largest eigenvalues of an arbitrary positive definite matrix A. We assume that there exist a positive constant 1<c<∞, such that c−1<λminΣ≤λmaxΣ<c.(C3) Regression coefficients. We assume thatβ≤Cβ for some constant Cβ>0, where . denotes the standard L2 norm. Also assume βmin≥νβn−ξmin for some ξmin>0, with βmin=minj∈Tβj.(C4) Divergence speed ofd and d0. We assume logd=Onξ and d0=Onξ0. In other words, there exists constants ξ, ξmin, and ν, such that logd≤νnξ, d0≤νnξ0, and ξ+6ξ0+12ξmin<1.


## 3.2. OMP Algorithm

Under the assumption that the true setT exists, our main objective is to discover all relevant predictors consistently. To this end, we consider the following OMP algorithm (Algorithm1):Algorithm 1: The OMP algorithm.
Step 1 (Initialization). SetS0=∅. Set the residual r0=y.Step 2 (Forward selection).(i) (2.1) Evaluation. In thekth step k≥1, we are given Sk−1. Then, for every j∈F∖Sk−1, we computeXjTrk−1,whererk−1=In−HSk−1y, HSk−1=XSk−1XSk−1TXSk−1−1XSk−1T is the projection onto the linear space spanned by the elements of XSk−1 and In∈Rn×n is the identity matrix.(ii) (2.2) Screening. We then findak=argmaxj∈F∖Sk−1XjTrk−1,and updateSk=Sk−1∪ak accordingly.Step 3 (Solution path). Iterating Step 2 forn times leads a total of n nested candidate models. We then collect those models by a solution path S=Sk:1≤k≤n with Sk=a1,…,ak.


## 3.3. Theoretical Properties

To prove Theorem1, the following lemmas are needed. For convenience, we define Σ^=n−1XTX. Moreover, Σ^M and ΣM are the submatrices of Σ^ and Σ (corresponding to M), respectively.Lemma 1.
If we have1/c<λminΣ<λmaxΣ<c, then eigenvalues of the submatrix of Σ, that is Σ11, are also bounded. Suppose Σ11 contains m variables, that is Σ11 is a m×m matrix. Moreover, with probability tending to 1, we have(5)1/c<λminΣ^11<λmaxΣ^11<c,as long as m=On1−ξ/3.The Proof of Lemma 1 can be found in Appendix.Lemma 2.
Assume conditions (C1), (C2), (C4), andm=On1−ξ/3. Then, with probability tending to 1, we have(6)τmin≤minM≤mλminΣ^M≤maxM≤mλmaxΣ^M≤τmax.Note that proof of Lemma2 can be found in the study by Wang (2009), with only slight changes.Before the theorems are established, we follow Wang (2009)’s idea on screening consistency of a solution path and define the solution pathS to be screening consistent, if(7)PT⊂Sk∈Sfor some1≤k≤n⟶1.Then OMP’s screening consistency can be formally established by the following theorem.Theorem 1.
Under model (4) and conditions (C1)–(C4), we have as n⟶∞(8)PT⊂SKnξ0+4ξmin⟶1,where the constant K=2τmaxτmin−2Cβ2νβ−4ν is independent of n, the constants τmax, τmin, Cβ, νβ, and ν are defined in conditions (C2)–(C4), and t is the smallest integer no less than t.Theorem1 proves that within Knξ0+4ξmin steps, all relevant predictors will be identified by the OMP algorithm. This number of steps is much smaller than the sample size n under condition (C4). In particular, if the dimension of the true model is finite with ξ0=ξmin=0, only a finite number of steps are needed to discover the entire relevant variable set.Furthermore, Theorem1 provides a theoretical basis for OMP, which enables us to empirically select the best model from S. On the other hand, the solution path S contains a total of n nested models. To further select relevant variables from the solution path S, the following BIC (Chen and Chen, 2008) is considered,(9)BICM=logσ^M2+n−1Mlogn+2logd,where M is an arbitrary candidate model with M≤n, σ^M2=n−1RSSM=n−1YTIn−HMY, and HM=XMXMTXM−1XMT. We then select the best model S^=Sm^, where m^=argmin1≤m≤nBICSm. We typically do not expect S^ to be selection consistent (i.e., PS^=T⟶1). However, we are able to show that S^ is indeed screening consistent.Theorem 2.
Under model (4) and conditions (C1)–(C4), then as n⟶∞(10)PT⊂S^⟶1.Definekmin=min1≤k≤nk:T⊂Sk. By Theorem 2.1, we know that kmin satisfies kmin≤Knξ0+4ξmin with probability tending to 1. Therefore, our aim is to prove that Pm^≤kmin⟶0 as n⟶∞. Then the theorem conclusion follows. Equivalently, it suffices to show that(11)Pmin1≤k<kminBICSk−BICSk+1>0⟶1.A detailed proof can be found in Appendix.


## 4. Numerical Analysis

### 4.1. Simulation Setup

For a reliable numerical comparison, the following three simulation examples on OMP algorithm are presented, to examine the finite performances of the screening consistent property of OMP. For each parameter setup, a total ofN=100 simulation replications are conducted.LetS^k=j:β^jk≠0 be the model selected in the kth simulation replications and the corresponding average model size =100−1∑kS^k. Recall T represents the true model, and the coverage probability=100−1∑kIS^k⊃T, which measures how likely all relevant variables is discovered by the method, is evaluated. This defined coverage probability characterizes the screening property of a particular method.To characterize the capability of a method in producing sparse solutions, we define(12)Percentage of Correct Zeros%=1d−d01100∑k=1100∑j=1dIβ^jk=0×Iβj=0.To characterize the method’s underfitting effect, we further define(13)Percentage of Incorrect Zeros%=1d01100∑k=1100∑j=1dIβ^jk=0×Iβj≠0.If all sparse solutions are correctly identified for all irrelevant predictors and no sparse solution is mistakenly produced for all relevant variables, the true model is perfectly identified, that isS^k=T. To measure such a performance, we define the percentage of correctly fitted (%)=100−1∑kIS^k=T, which characterizes the selection consistency property of a particular method.As we need to know which variables are truly relevant or irrelevant, we create sparse regression vectors by settingβi=0 for all i=1,…,d, except for a chosen set T of coefficients, where βi is defined in advance for every 1≤i≤d0. Moreover, the noise vector ϵ1,…,ϵn is chosen i.i.d. N0,1. Note that all the simulation runs are conducted in MATLAB.Example 1. (independent predictors).
This is an example borrowed from Fan and Lv [7]. Xi is generated independently according to a standard multivariate normal distribution. Thus, different predictors are mutually independent. n,d,d0=100,5000,8 with βj=−1Uj4lognn+Zj, where Uj is a binary random variable with PUj=0.4 and Zj is a standard normal random variable.Example 2. (autoregessive correlation).
Xi is generated from a multivariate normal distribution with mean 0 and CovXij1,Xij2=0.5j1−j2. This is called an autoregressive type correlation structure. Such type of correlation structure might be useful if a natural order exists among the predictors. As a consequence, the predictors with large distances in order are expected to be mutually independent approximately. This is an example from Tibshrani [4] with n,d,d0=100,5000,3. In addition, the first, fourth, and seventh components of β are set to be 3, 1.5, and 2, respectively.Example 3. (grouped variables).
Xi is generated by the following rule. Xij=3/20Z1+17/20ϵx,j for j∈1,…,d0, Xij=19/20Z2+1/20ϵx,j for j∈d0+1,…,d0+5, and Xij=ϵx,j otherwise, where Z1∼N0,1, Z2∼N0,1, and ϵx,j∼N0,1 are independent. This creates within-group correlations of ρij=0.15 for i,j∈1,…,d0 and ρij=0.95 for i,j∈d0+1,…,d0+5. This example presents an interesting scenario where a group of significant variables are mildly correlated and simultaneously a group of insignificant variables are strongly correlated. The settings are similar to those in Example 2. n,d,d0=100,5000,3. In addition, the three nonzero components of β are set to be 3, 1.5, and 2, respectively.


### 4.2. Simulation Results for OMP Screening Consistent Property

Finite sample performance of OMP screening consistent property is investigated based on the abovementioned three examples in Section4.1. Simulation results are presented in Table 1.Table 1
Simulation summary for OMP with (n, d) = (100, 5000).
d0Coverage probability (%)Percentage of correct zeros (%)Percentage of incorrect zeros (%)Percentage of correctly fitted (%)Average model sizeIP810099.91.6947.94AC310010001003GV39695.92343.84First of all, simulation results for the independent predictor example are in good performance in terms of screening selection consistency for OMP. In other words, we have 100% coverage probability, which means all relevant variables can be discovered by OMP method. In addition, 94% of correctly fitted denotes that BIC selects the true set of variables correctly 94 times out of 100 simulation replications. This result is not surprising since Zhang [14] pointed out that OMP can select features or variables consistently under a certain irrepresentable condition. Furthermore, the percentage of correct zeros and the percentage of incorrect zeros are 99.9% and 1.6%, respectively. Last but not least, the average model size is 7.94, which is slightly below d0=8.Furthermore, simulation results for autoregressive correlation example are in very good performance in terms of screening selection consistency for OMP. Both of the coverage probability and the percentage of correctly fitted are 100%. Especially 100% of correctly fitted denotes BIC selects the true set of variables correctly 100 times out of 100 simulation replications. This is good news since the number of nonzeroβjd0 is 3, which is a very sparse representation given d=5000. On top of that, the percentage of correct zeros and the percentage of incorrect zeros are 100% and 0%, respectively. Last but not least, the average model size is 3. Therefore, it seems that our OMP algorithm works pretty well under this autoregressive correlation setup with a sparse representation of β.Third, simulation results for grouped variables example are in worst performance among all the three examples in terms of screening selection consistency for OMP. However, the performance itself is not bad. Coverage probability is 96%, meaning that not all the relevant predictors can be discovered by OMP algorithm in some of the simulation replications. In addition, 34% of correctly fitted denotes that BIC selects the true set of variables correctly only 34 times out of 100 simulation replications. On top of that, the percentage of correct zeros and the percentage of incorrect zeros are 95.9% and 2%, respectively. Last but not least, the average model size is 3.84.Besides a summary of simulation results of OMP algorithm, three plots are presented in Figure1. For each of the three examples, one particular plot of number of variables included in the final model versus BIC is extracted for reference. These graphs are not representable as a whole; however, they do provide trends of BIC casewise. Take BIC of Example 1 as an example. Please refer to Figure 1(a). BIC decreases as the number of variables included in the model increases. BIC reaches a local minimum when the number of variables included hits d0=8. Thereafter, BIC increases as the number of variables included in the model increases until the number of variables gets near to the sample size. This is not surprising since BIC decreases as the model complexity increases. Similar trends can be observed for Example 2 and Example 3. Please refer to Figures 1(b) and 1(c).Figure 1
OMP simulation results of BIC Trends:X-axis for number of predictors included in the model and Y axis for BIC. (a): Example 1, (b): Example 2, and (c): Example 3.
(a)(b)(c)One possible suggestion for OMP algorithm is that after the screening process withn candidate models, only n/2 BIC candidate models for minimum BIC values are compared. By doing so, computational time can be saved without loss of correctness of screening consistent property of OMP.In conclusion, finite simulation performances in terms of screening selection consistency for OMP are good under all the three examples. Those performances support our theories proposed in Section3.3.


## 4.1. Simulation Setup

For a reliable numerical comparison, the following three simulation examples on OMP algorithm are presented, to examine the finite performances of the screening consistent property of OMP. For each parameter setup, a total ofN=100 simulation replications are conducted.LetS^k=j:β^jk≠0 be the model selected in the kth simulation replications and the corresponding average model size =100−1∑kS^k. Recall T represents the true model, and the coverage probability=100−1∑kIS^k⊃T, which measures how likely all relevant variables is discovered by the method, is evaluated. This defined coverage probability characterizes the screening property of a particular method.To characterize the capability of a method in producing sparse solutions, we define(12)Percentage of Correct Zeros%=1d−d01100∑k=1100∑j=1dIβ^jk=0×Iβj=0.To characterize the method’s underfitting effect, we further define(13)Percentage of Incorrect Zeros%=1d01100∑k=1100∑j=1dIβ^jk=0×Iβj≠0.If all sparse solutions are correctly identified for all irrelevant predictors and no sparse solution is mistakenly produced for all relevant variables, the true model is perfectly identified, that isS^k=T. To measure such a performance, we define the percentage of correctly fitted (%)=100−1∑kIS^k=T, which characterizes the selection consistency property of a particular method.As we need to know which variables are truly relevant or irrelevant, we create sparse regression vectors by settingβi=0 for all i=1,…,d, except for a chosen set T of coefficients, where βi is defined in advance for every 1≤i≤d0. Moreover, the noise vector ϵ1,…,ϵn is chosen i.i.d. N0,1. Note that all the simulation runs are conducted in MATLAB.Example 1. (independent predictors).
This is an example borrowed from Fan and Lv [7]. Xi is generated independently according to a standard multivariate normal distribution. Thus, different predictors are mutually independent. n,d,d0=100,5000,8 with βj=−1Uj4lognn+Zj, where Uj is a binary random variable with PUj=0.4 and Zj is a standard normal random variable.Example 2. (autoregessive correlation).
Xi is generated from a multivariate normal distribution with mean 0 and CovXij1,Xij2=0.5j1−j2. This is called an autoregressive type correlation structure. Such type of correlation structure might be useful if a natural order exists among the predictors. As a consequence, the predictors with large distances in order are expected to be mutually independent approximately. This is an example from Tibshrani [4] with n,d,d0=100,5000,3. In addition, the first, fourth, and seventh components of β are set to be 3, 1.5, and 2, respectively.Example 3. (grouped variables).
Xi is generated by the following rule. Xij=3/20Z1+17/20ϵx,j for j∈1,…,d0, Xij=19/20Z2+1/20ϵx,j for j∈d0+1,…,d0+5, and Xij=ϵx,j otherwise, where Z1∼N0,1, Z2∼N0,1, and ϵx,j∼N0,1 are independent. This creates within-group correlations of ρij=0.15 for i,j∈1,…,d0 and ρij=0.95 for i,j∈d0+1,…,d0+5. This example presents an interesting scenario where a group of significant variables are mildly correlated and simultaneously a group of insignificant variables are strongly correlated. The settings are similar to those in Example 2. n,d,d0=100,5000,3. In addition, the three nonzero components of β are set to be 3, 1.5, and 2, respectively.


## 4.2. Simulation Results for OMP Screening Consistent Property

Finite sample performance of OMP screening consistent property is investigated based on the abovementioned three examples in Section4.1. Simulation results are presented in Table 1.Table 1
Simulation summary for OMP with (n, d) = (100, 5000).
d0Coverage probability (%)Percentage of correct zeros (%)Percentage of incorrect zeros (%)Percentage of correctly fitted (%)Average model sizeIP810099.91.6947.94AC310010001003GV39695.92343.84First of all, simulation results for the independent predictor example are in good performance in terms of screening selection consistency for OMP. In other words, we have 100% coverage probability, which means all relevant variables can be discovered by OMP method. In addition, 94% of correctly fitted denotes that BIC selects the true set of variables correctly 94 times out of 100 simulation replications. This result is not surprising since Zhang [14] pointed out that OMP can select features or variables consistently under a certain irrepresentable condition. Furthermore, the percentage of correct zeros and the percentage of incorrect zeros are 99.9% and 1.6%, respectively. Last but not least, the average model size is 7.94, which is slightly below d0=8.Furthermore, simulation results for autoregressive correlation example are in very good performance in terms of screening selection consistency for OMP. Both of the coverage probability and the percentage of correctly fitted are 100%. Especially 100% of correctly fitted denotes BIC selects the true set of variables correctly 100 times out of 100 simulation replications. This is good news since the number of nonzeroβjd0 is 3, which is a very sparse representation given d=5000. On top of that, the percentage of correct zeros and the percentage of incorrect zeros are 100% and 0%, respectively. Last but not least, the average model size is 3. Therefore, it seems that our OMP algorithm works pretty well under this autoregressive correlation setup with a sparse representation of β.Third, simulation results for grouped variables example are in worst performance among all the three examples in terms of screening selection consistency for OMP. However, the performance itself is not bad. Coverage probability is 96%, meaning that not all the relevant predictors can be discovered by OMP algorithm in some of the simulation replications. In addition, 34% of correctly fitted denotes that BIC selects the true set of variables correctly only 34 times out of 100 simulation replications. On top of that, the percentage of correct zeros and the percentage of incorrect zeros are 95.9% and 2%, respectively. Last but not least, the average model size is 3.84.Besides a summary of simulation results of OMP algorithm, three plots are presented in Figure1. For each of the three examples, one particular plot of number of variables included in the final model versus BIC is extracted for reference. These graphs are not representable as a whole; however, they do provide trends of BIC casewise. Take BIC of Example 1 as an example. Please refer to Figure 1(a). BIC decreases as the number of variables included in the model increases. BIC reaches a local minimum when the number of variables included hits d0=8. Thereafter, BIC increases as the number of variables included in the model increases until the number of variables gets near to the sample size. This is not surprising since BIC decreases as the model complexity increases. Similar trends can be observed for Example 2 and Example 3. Please refer to Figures 1(b) and 1(c).Figure 1
OMP simulation results of BIC Trends:X-axis for number of predictors included in the model and Y axis for BIC. (a): Example 1, (b): Example 2, and (c): Example 3.
(a)(b)(c)One possible suggestion for OMP algorithm is that after the screening process withn candidate models, only n/2 BIC candidate models for minimum BIC values are compared. By doing so, computational time can be saved without loss of correctness of screening consistent property of OMP.In conclusion, finite simulation performances in terms of screening selection consistency for OMP are good under all the three examples. Those performances support our theories proposed in Section3.3.


## 5. Conclusion and Future Research

To conclude, this article shows the theoretical proof that OMP can identify all relevant predictors consistently, even if the predictor dimension is considerably larger than the sample size. In particular, if the dimension of the true model is finite, OMP might discover all relevant predictors within a finite number of steps. In other words, sure screening property can be guaranteed under the four technical conditions. Given the sure screening property, the recently proposed BIC of Chen and Chen (2008) can be used to practically select the best candidate from the models generated by the OMP algorithm. The resulting model is good in the sense that many existing variable selection methods, such as adaptive Lasso and SCAD, can be applied directly to increase the estimation accuracy. Compared with the traditional orthogonal matching pursuit method, the resulting model can improve prediction accuracy and reduce computational cost by screening out the relevant variables.The abovementioned variable selection procedure only considers the fixed effect estimates in the linear models. However, in real life, a lot of existing data have both the fixed effects and random effects involved. For example, in the clinic trials, several observations are taken for a period of time for one particular patient. After collecting the data needed for all the patients, it is natural to consider random effects for each individual patient in the model setting since a common error term for all the observations is not sufficient to capture the individual randomness. Future research may include random effects in the model by imposing penalized hierarchical likelihood algorithm for accurate variable selection.



---

*Source: 6446903-2022-08-01.xml*