# A Two-Stage Regularization Method for Variable Selection and Forecasting in High-Order Interaction Model

Authors: Yao Dong; He Jiang
Journal: Complexity
 (2018)

---

## Abstract

Forecasting models with high-order interaction has become popular in many applications since researchers gradually notice that an additive linear model is not adequate for accurate forecasting. However, the excessive number of variables with low sample size in the model poses critically challenges to predication accuracy. To enhance the forecasting accuracy and training speed simultaneously, an interpretable model is essential in knowledge recovery. To deal with ultra-high dimensionality, this paper investigates and studies a two-stage procedure to demand sparsity within high-order interaction model. In each stage, square root hard ridge (SRHR) method is applied to discover the relevant variables. The application of square root loss function facilitates the parameter tuning work. On the other hand, hard ridge penalty function is able to handle both the high multicollinearity and selection inconsistency. The real data experiments reveal the superior performances to other comparing approaches.


---

## Body

## 1. Introduction

Sparse representation has attracted a large deal of attention from researchers in different scientific fields including financial engineering and risk management computational biology and machine learning community [1]. For instance, in spectrum data analysis, scientists are exploring how to select few exact frequencies on which the beam energy is concentrated to avoid lower leakage problem [2]. The main goal of demanding sparsity is to discover quite a few numbers of relevant variables from a large pool of candidate variables which are essential when establishing an interpretable forecasting model. A model with high precision, low model complexity, and sensitivity to the model parameters can not only avoid the expensive cost of excessive supply but also reduce the time cost, service interruption, and resource waste [3]. Furthermore, researchers have noticed that a linear additive model using main effects only cannot provide accurate forecasting results. This motivates them to turn their attentions to high-order interaction model which considers interaction terms between variables. Interaction model has advantage over a linear model which does not consider possible relationship between variables. In microarray data analysis, for example, biologists are particular interested in gene-gene interactions and gene-environment interactions than a single gene because interactions of single nucleotide polymorphisms (SNPs) are more important in cancer diagnosis [4]. However, new studies have also revealed the interactive effect of multiple environment factors. For example, a child with a poor quality environment would be more sensitive to a poor environment as an adult which ultimately led to higher psychological distress scores. This depicts a three-way interaction term between gene and environment. Furthermore, interaction with multiple genes is also considered to be essentially important in the molecular analysis [5]. This motivates us to consider a high-order interaction model which includes both two-way and three-way interaction terms in this paper:(1)y=b0+∑jbjxj+∑j,kϕjkxj⊙xk+∑j,k,lψjklxj⊙xk⊙xl+εwhere b0 is the intercept term which vanishes after centering, xj⊙xk⊙xl represents the Hadamard product of xj, xk, and xl which denotes the three-way interaction term. The noise term follows Gaussian distribution. bj, ϕjk, and ψjkl are corresponding coefficients for main effects, two-way, and three-way interaction terms, respectively. Obviously, this high-order interaction model is more complicate than two-way interaction model.Although high-order interactions play an important role in depicting the nonlinear relationship among variables, sparse representation of high-order interaction model is quite challenging since the number of interaction terms is huge (of orderO(p2) or even O(p3)) if there are p main effects in the model. For instance, even with 102 variables in the model, a high-order interaction model has to consider 104 or 106 variables which make variable selection procedure more difficult. Furthermore, computation cost is another issue because of the huge number of variables in the high-order interaction model. Even worse, severe multicollinearity exists among variables caused by the high dimensionality. There is a large body of literature on sparse representation methods. Subset selection methods based on L0 type penalty have been discovered for decades. These methods including Akaike information criterion (AIC) [6], Bayesian information criterion (BIC) [7], and mallow’s Cp [8] which are NP hard since all the possible combinations of variables are taken into account. Namely, if there arem variables, 2m subset models have to be checked. To improve the computational efficiency, sparse representation approaches based on convex optimization are proposed. For instance, Tibshrani (1996) proposed the least absolute shrinkage and selection operator (LASSO) which applies L1 type penalty to pursue sparsity. In wavelet study, it is famous as basis pursuit method [9]. It is more computational efficient than subset selection methods. Efron (2004) studied the least-angle regression (LAR) which solved the same problem as LASSO but in a forward selection manner. They have shown that the solution path of LAR is pretty similar as that of LASSO. Yuan and Lin (2006) investigated group LASSO method which selected group variables based on “all in and all out” principle. To overcome the shortcoming of LASSO and further improve the selection accuracy, Zou (2006) studied adaptive LASSO using weights for each variable to adjust the sparsity degree. Friedman (2012) summarized and provided a review of sparse regression and classification methods [10].However, convex optimization has drawbacks since it is difficult to discover all the relevant variables if the model coherence is sufficiently high. To overcome the difficulties, nonconvex penalties including smoothly clipped absolute deviation (SCAD) [11], minimax concave penalty (MCP) [12], and hard ridge (HR) [13, 14] are proposed to boost the accuracy. Ye et al. (2017) studied sparse least squares support vector machine with Lp norm function applied to select the relevant support vectors [15]. To facilitate parameter tuning work and increase the forecasting accuracy, Jiang (2018) studied square root nonconvex optimization (SRNO) and showed its superior performance [2] via simulation and real data studies.Sparse representation in the two-way interaction model has been studied for years. Bach (2008) explored hierarchical multiple learning and studied its selection property [16]. Zhao et al. (2009) proposed composite absolute penalty (CAP) which considers a mixture of Lp and Lq norm functions which are given by ϕjkp+bj,bk,ϕjkq [17]. Choi et al. (2010) studied SHIM method considering ϕjk=ρbjbk. They also designed a nonconvex optimization to perform sparsity representation. However, the number of features is only 10 which is pretty small. Thus, there is no guarantee that their method can work in the high dimensional setup [18]. Bickel et al. (2010) used a two-step algorithm to purse sparsity. In the first step, LASSO was applied and backward procedure was used in the second step [19]. Wu et al. (2010) proposed a stagewise procedure which screened the main effects in the first step and the selected variable will be cleaned based on t-test statistic [20]. Radchenko and Jame (2010) investigated VANISH constructing the main effects and interactions from two small sets of orthogonal functions [21]. Bien et al. (2013) studied hierarchical LASSO (HL) which considered a L1 optimization problem using a nonconvex constraint ϕj1≤bj [22]. She et al. (2014) studied GRESH method which enforced sparsity to main effects and two-way interactions simultaneously [23]. Hao et al. (2014) screened the main effects using forward regression in the first step and fitted the quadratic model with selected main effects and associated interactions [24]. Simon and Hastie (2015) proposed “Glinternet” to discover the two-way interaction terms between categorical variables and continuous variables [25]. Yan and Bien (2017) advocated VARXL method to implement sparse representation in a high dimensional time series data [26]. To the best of our knowledge, most of variable selection methods mentioned above only consider two-way interaction terms. Few studies have been done on the high-order interaction model with both two-way and three-way interaction terms. To tackle this issue, this paper explores sparse representation in a high-order interaction model. The main contributions of this paper are given as follows:(i)
A two-stage square root hard ridge method (TSSRHR) for sparse representation in high-order interaction model is proposed.(ii)
In computation, a fast and simple-to-implement algorithm is designed.(iii)
In theory, we have shown the prediction error bound and estimation error bound of the TSSRHR estimate.(iv)
Real data examples are shown to demonstrate the efficiency and superior performance of the proposed method over other existing competitors.The rest of this paper is organized as follows. Section2 represents the general framework of TSSRHR. Section 3 designs a fast and simple-to-implement algorithm with theoretical guarantee of its convergence and optimality. The theory part of TSSRHR is given in Section 4. In Section 5, real data analyses are shown to exhibit superior performance of the proposed method over other approaches. The conclusion of this paper is given in Section 6.


## 2. Two-Stage Square Root Hard Ridge Method

We investigate and study the square root hard ridge (SRHR) method which combines the benefits of square root loss function and nonconvex hard ridge penalty function for high-order interaction selection. We describe the two-stage SRHR algorithm at first. At stage 1, only main effects are selected by SRHR with all of high-order interactions including two-way and three-way interactions kept out. The selected variables are put in the set A. At stage 2, the selected set in stage 1.A is expanded by adding both associated two-way and three-way interactions within A. SRHR is applied again to implement variable selection task.(i)
Stage 1. DefineG=[p]=1,…,p which represents all the main effects. Select the important effects using SRHR method by considering the following optimization problem:(2)minb∈Rp∑i=1nyi-b0-∑j∈pbjxj21/2+λ1221+η1∑j∈p1bj≠0+η12∑j∈pbj2whereλ1 and η1 are two regularization parameters which are supposed to be determined. The index for selected variables are given in A and the corresponding solution path is given by Sλ1,η1(1).(ii)
Stage 2. UpdateG=A∪(j,k);j∈A,k∈A∪(j,k,l);j∈A,l∈A,k∈A. Select the important main effects, two-way and three-way interaction using SRHR method again by considering the following problem (for notation simplicity, xj⊙xk and xj⊙xk⊙xl are denoted by xjxk and xjxkxl, respectively):(3)minb,ϕ,ψ∑i=1nyi-b0-∑j∈Abjxj-∑j,k∈A⊗Aϕjkxjxk-∑j,k,l∈A⊗A⊗Aψjklxjxkxl21/2+λ2221+η2∑j∈A1bj≠0+∑j,k∈A⊗A1ϕjk≠0+∑j,k,l∈A⊗A⊗A1ψjkl≠0+η22∑j∈Abj2+∑j,k∈A⊗Aϕjk2+∑j,k,l∈A⊗A⊗Aψjkl2where1(·) denotes the indicator function, ⊗ represents the Cartesian product, and λ2 and η2 are two regularization parameters which need to be selected carefully. The solution path is provided by Sλ2,η2(2) and the selected variable are given in B.A penalized sparse representation method considers the sum of two components: the loss function and the penalty function. The loss function often applies the negative log-likelihood of the error and in Gaussian setup, it is a squared error loss function which involves the noise levelσ which is difficult to estimate. Specifically, in a Gaussian linear regression setting y=Xβ+ε, the negative likelihood function of y is -logfyXβ=y-Xβ22/σ2+hσ+c, where h(σ) is a function of σ and c is a constant. SRHR uses square root error loss function which can facilitate the parameter tuning work by avoiding estimating the noise level σ or using other forms of estimates [27]. Hard ridge penalty has two parts; it is obvious that the first part which is L0 norm penalized the number of nonzero elements using the regularization parameter λ. The second part is L2 type ridge penalty. It can compensate the shrinkage given by the L0 norm and is able to handle the high collinearity of the high-order interaction model. Therefore, SRHR combines the advantages of square root loss function and hard ridge penalty.The choices of regularization parameters (λ1,λ2,η1, and η2) play an important role in balancing a model’s in sample fit and out-of-sample forecasting ability. We grid λ1 starting from the largest value λ1,max=maxXTy such that all the coefficients are shrunk to zeros. Then let the smallest value λ1,min be a small proportion of λ1,max, that is, λ1,min=0.05λ1,max and the grid values are generated between λ1,max and λ1,min exponentially. For λ2, the same strategy is also applied. The optimal values for η1 and η2 are obtained from a grid of values 1e-5,1e-4,1e-3,1e-2,1e-1,1e+1.To choose the optimal regularization parameters, extended Bayesian information criterion (EBIC) [28] and high dimensional Bayesian information criteria (HDBIC) [29] are applied in this paper. Particularly, EBIC is defined as(4)EBICA=1n-A∑i=1nyi-b^0-∑j∈Ab^jxj2+n-1Alog⁡n+log⁡p,and HDBIC is defined as(5)EBICB=1n-B∑i=1nyi-b^0-∑j∈Ab^jxj-∑j,k∈A⊗Aϕ^jkxjxk-∑j,k,l∈A⊗A⊗Aψ^jklxjxkxl2+n-1Blog⁡n+logA+A2+A3where A and B are represent the cardinalities of A and B. In the first stage, λ1 and η1 will be selected by the EBIC since large number of main effects are considered. It worth mentioning that different from [24], we did not apply BIC in the second stage. Although the number of variables is reduced dramatically, we still use HDBIC which is a high dimensional version of BIC because both two-way and three-way interaction terms are considered which still cause the dimensions of the model to be high.


## 3. Algorithm

In this section, for notation simplicity, two design matrices are designed as follows:(6)X-=x12,⋯,x1xp,xpx1,⋯,xp2∈Rn×p2,and(7)X˘=x13,x12x2⋯,x12xp,⋯,xp2x1,⋯,xp3∈Rn×p3Furthermore, let X~=X,X-,X˘. Clearly, X- represents the design matrix for two-way interaction terms, X˘ denotes the design matrix for three-way interaction terms, and X~ is the overall design matrix for main effects, two-way interaction terms, and three-way interaction terms. For the coefficients to be estimated, let Φ and Ψ be the coefficient matrix for two-way and three-way interaction terms, respectively. Then (1) can be written as(8)y=b0+Xb+X-vecΦ+X˘vecΨ+εwhere vec(A)=a11,a12,⋯,ann for any A=aijn×n. Define Ω=b,vec(Φ)T,vec(Ψ)TT. To design an algorithm for TSSRHR, we consider the following optimization problem in the first stage:(9)minby-Xb2+λ1221+η1b0+η12b22and in the second stage, we consider(10)minΩy-X~Ω2+λ2221+η2Ω0+η22Ω22To solve (9) and (10) accurately and efficiently, threshold rule is defined at first.Definition (threshold rule). A threshold rule is a real-valued function Θ(t;λ) defined for -∞<t<∞ and 0≤λ<∞ such that(1)
Θ(-t;λ)=-Θ(t;λ);(2)
Θ(t;λ)≤-Θ(t′;λ) for t≤t′;(3)
limt→+∞Θ(t;λ)=∞; and(4)
0≤Θ(t;λ)≤t for 0≤t≤∞.From the definition, it is easy to see thatΘ(t;λ) is an odd monotone unbounded shrinkage rule for t, at any λ. A vector version of Θ is defined component-wise if either t or λ is replaced by a vector. Now define the HR thresholding rule, which is closely associated with the HR penalty, as (11)ΘHRt;λ,η=t1+η1t≥λwhere η>0, λ>0 are two regularization parameters. Using HR thresholding function, (9) and (10) can be solved based on the following iterations:(12)The1ststage:bj+1=Θ→HRbj+XTXbj-y;λ1y-Xbj2,η1y-Xbj2and(13)The2ndstage:Ωk+1=Θ→HRΩk+X~TX~Ωk-y;λ2y-X~Ωk2,η2y-X~Ωk2where Θ→HR is a vector version of HR thresholding function with t replaced by a vector.Algorithm1 exhibits the details of TSRHR algorithm with ω(-1)=0, ω(0)=1, b(0)=b(-1), and Ω(0)=Ω(-1)defined. Consider the nonconvex optimization is not easy to solve using ADMM algorithm [30]. HR thresholding function is applied in the algorithm. The basic structure of this algorithm follows coordinate decent algorithm. The error tolerance of tol1, tol2 and maximum number of iterations m1 and m2 are determined based on trial and error. Usually, the default values of tol1, tol2 and m1 and m2 are 1e-4, 1e-4, 500, and 500, respectively. TSSRHR algorithm has some advantages. Firstly, it does not involve any matrix inversion. Secondly, ΘHR(t;λ,η) is applied to enforce elementwise sparsity. The step sizes τ1 and τ2 are able to guarantee the convergence of the algorithm. The theoretical results show that τ1≥X2 and τ2≥X~2 provide satisfactory results. The step sizes do not depend on line search which takes a lot of computation time. The TSSRHR algorithm is simple to implement and computational efficient. To further speed up the convergence, a fast iterative shrinkage thresholding algorithm (FISTA) [31] is used to obtain the TSSRHR estimate. The following theorem provides the theoretical guarantee of convergence of Algorithm 1.Algorithm 1:The TSRRHR algorithm.
Inputs:Z=(X,y)=x~1,y1,⋯,x~N,yN:x~k∈Rp,X∈Rn×p,y∈Rn the dataset m 1: maximum number of iterations used in the SRHR algorithm of the first stage; m 2: maximum number of iterations used in the SRHR algorithm of the second stage; t o l 1: the error tolerance used in the SRHR algorithm of the first stage; t o l 2: the error tolerance used in the SRHR algorithm of the second stage; τ 1: scale parameter used in the SRHR algorithm of the first stage; τ 2: scale parameter used in the SRHR algorithm of the second stage; Output:The forecasting test error and selected pattern; ( 1 ) Randomly divide the original data Z into the training dataset (Xtrain,ytrain) and test dataset(Xtst,ytst). The first stage using SRHR algorithm: ( 2 ) Generate grid values of λ1G=λ1,uu=1s1 and η1G=η1,vv=1t1. ( 3 )  for u = 1 to s 1 ( 4 )  for v = 1 to t 1 ( 5 ) Initialization: j←0,b(j)←0 ( 6 ) Scaling: Xtrain←Xtrain/τ1,ytrain←ytrain/τ1 ( 7 )  while b ( j + 1 ) - b ( j ) 2 < t o l 1 or j>m1 do ( 8 )  Step  1. ξ1(j)←b(j)+ω(j-1)/ω(j)b(j)-b(j-1) ( 9 )  Step  2. ζ1(j)←ξ2(j)+XtrainTytrain-Xtrainξ1(j) ( 10 )  Step  3. b(j+1)←Θ→HRζ1(j);λ1,uy-X~b(j)2/τ12,η1,vy-X~b(j)2/τ12 ( 11 )  Step  4: ω(j+1)=1+1+4ω(j)2/2 ( 12 )  end while ( 13 )  end for ( 14 )  end for ( 15 ) Obtain the solution path A and the corresponding sparsity patterns based on EBIC criterion. The second stage using SRHR algorithm: ( 16 ) Generate the high-order interaction model (X~train,ytrain) and (X~tst,ytst) based on the sparsity pattern A. ( 17 ) Generate grid values of λ2G=λ2,uu=1s2 and η2G=η2,vv=1t2. ( 18 )  for u = 1 to s 2 ( 19 )  for v = 1 to t 2 ( 20 ) Initialization: k←0,Ω(k)←0 ( 21 ) Scaling: X~train←X~train/τ2,ytrain←ytrain/τ2 ( 22 )  while Ω ( k + 1 ) - Ω ( k ) 2 < t o l 2 or k>m2 do ( 23 ) Step 1. ξ2(k)←Ω(k)+ω(k-1)/ω(k)Ω(k)-Ω(k-1) ( 24 ) Step 2. ζ2(k)←ξ2(k)+X~trainTytrain-X~trainξ2(k) ( 25 )  Step 3. Ω(k+1)←Θ→HRζ2(k);λ2,uy-X~Ω(k)2/τ22,η2,vy-X~Ω(k)2/τ22 ( 26 ) Step 4: ω(k+1)=1+1+4ω(k)2/2 ( 27 )  end while ( 28 )  end for ( 29 )  end for ( 30 )  Obtain the solution path B and update the sparsity patterns G using HDBIC criterion. ( 31 )  Calculate the test error using test dataset (X~tst,ytst)


## 4. Prediction and Estimation

Notice thatP(Ω)=λ2Ω0/2(1+η)+ηΩ22/2 does not follow the triangle inequality P(a+b)≤P(a)+P(b) but follows the following inequality P(a+b)≤2P(a)+P(b) for any vectors a,b. This poses critically challenges to the derivation the forecasting and estimation error bounds. Namely, the ways of performing error rate analysis on square root lasso [27] are not applicable to SRHR because of its nonconvexity. Therefore, a novel approach is urgently needed to derive the forecasting and estimation error bounds.Assumption 1 (overfitting assumption).
Suppose thaty-X~Ω^2=0, then Ω^ does the tasks as the Ordinary Least Square (OLS) estimator Ω^OLS=(X~TX~)-TX~Ty.Assumption 2 (weak decomposability assumption).
A normM in Rp is called weakly decomposable for an index set J⊂[p], if there exists a norm MJcon Rp-J such that(14)MΩ≥MΩJ+MJcΩJc∀Ω∈RpFurthermore, a set J is called allowed set if M is a weakly decomposable norm for this set.We can show that hard ridge penalty is weakly decomposable as follows: defineM(Ω)≔λ2Ω0/2(1+η)+ηΩ22 and MJc(Ω)≔λ2Ω0/2(1+η)+ηΩ22. Then we can obtain weak decomposability for any set J(15)MΩJ+MJcΩJc≤MΩThis is due to the weak decomposability property of the penalty and the penalty.Definition 3 (M-effective sparsity).
For an allowed set of a normM, and L>0 a constant, the M eigenvalue is defined as (16)δML,J≔minX~ΩJ-X~ΩJc2:MJcΩJc≤LMΩJThen theM-effective sparsity is defined as(17)ΓM2L,J≔1δM2L,J
M-effective sparsity is a variant of regularity condition used to establish the nonasymptotic oracle inequality, such as Restricted Eigenvalue condition [32], Mutual Coherence condition [33], Compatible Condition [34], Cone invertibility factor [35], and Restricted invertibility factor [36]. With Assumptions 1 and 2 satisfied, based on the Theorem 1 in [37], we can obtain the forecasting and estimation error bounds of the TSSRHR estimator in the following remark.Remark. Choose(18)λ=C12ntδ+2+logp-Jand(19)η=C22ntδ+2+log⁡p-Jwhere t=log⁡2/α with α>2/en/4 and δ=1-2t2/n. Under the overfitting assumption and weak decomposability for M(·)≔λ2·0/2(1+η)+η·/2, the following inequalities hold with probability 1-α(20)X~Ω^-X~Ω2≤C3λ2ΓM2L,J∗and(21)Ω^-Ω22≤C4λ2ΓM2L,J∗where C1,C2,C3, and C4 are some positive constants.This remark provides prediction error bound and estimation error bound for the TSSRHR estimate. The proof of remark follows Theorem 1 and Corollaries 1 and 2 in [37] directly. Thus, the details are omitted here.


## 5. Real World Data

In this section, twinsine signal data including 8 scenarios and microarray gene expression data about inbred mouse are treated as real data examples to test the performances of sparse representation methods. It worth mentioning both of these two real data examples are ultra-high dimensional data when three-way interaction terms are considered. Comparisons will be made between the proposed two-stage hard ridge (TSSRHR) and other two-stage sparse representation methods including two-stage LASSO (TSLASSO), two-stage SCAD (TSSCAD), and two-stage square root lasso (TSSRL) in terms of forecasting accuracy, sparsity, and computational efficiency. The unknown parameters in TSLASSO and TSSCAD are determined using 10-fold cross validation. TSSRL selects the unknown parameter based on a theoretical choice:λ=1.1G-1(1-0.05/2K) where K is the number of variables and G is the cumulative distribution function of Gaussian distribution. Here we remind that TSSRHR selects the regularization parameters by EBIC and HDBIC (cf. Section 2). The forecasting model will be trained with 75% of the samples and the remaining 25% will be used to calculate the test errors. The entire procedure will be repeated for 150 times and the median of evaluation criteria including number of nonzero elements (NZ), mean square error (MSE), root mean square error (RMSE), and computational cost (CC) will be reported.


### 5.1. Twinsine Signal Data

Spectral estimation studying the distribution over frequencies is widely applied in speech coding and radar solar signal processing. Nevertheless, it is critically difficult and challenging when there are a large number of frequencies which causes high resolution problem. This type of problem is called super resolution spectral estimation. To avoid the power leakage and noisy observation caused by the super resolution, a novel method is urgently needed to choose the important frequencies which plays a dominant role in super-resolution estimation. The twinsine signal is generated from the target function:(22)yt=A1cos⁡2πf1t+ϕ1+A2cos⁡2πf2t+ϕ2+εtwhere A1=2, A2=3, ϕ1=π/3, f1=0.25HZ, and ε(t) follows white Gaussian noise with variance given by σ2. The frequence resolution is chosen as 0.002HZ to detect and distinguish the two sinusoidal components. fk=fmax·k/D for k=0,1,⋯,D and the frequency atoms are given by Xcos=cos(2πtifk)1≤i≤n,1≤k≤D and Xsin=sin(2πtifk)1≤i≤n,1≤k≤D. In our experiments, we will consider the following 8 different scenarios:(i)
Scenario 1:(n,D,σ2)=(100,300,1)(ii)
Scenario 2:(n,D,σ2)=(100,500,1)(iii)
Scenario 3:(n,D,σ2)=(100,1000,1)(iv)
Scenario 4:(n,D,σ2)=(100,5000,1)(v)
Scenario 5:(n,D,σ2)=(100,300,5)(vi)
Scenario 6:(n,D,σ2)=(100,500,5)(vii)
Scenario 7:(n,D,σ2)=(100,1000,5)(viii)
Scenario 8:(n,D,σ2)=(100,5000,5)(ix)
Scenario 9:(n,D,σ2)=(500,300,1)(x)
Scenario 10:(n,D,σ2)=(1000,300,1)The results of forecasting, variable selection, and computation are summarized in Table1 and Figures 1, 2, 3, and 4. It is not difficult to observe that TSSRHR obtains the most accurate outcome in all the scenarios including the most challenging scenario when D=5000 and σ2=5. TSSRL outperforms TSLASSO and TSSCAD method by delivering lower MSE and its forecasting results are comparable with TSSRHR in scenario 6. From the prospective of model sparsity, TSSRL achieves the most interpretable model followed by TSSRHR. TSLASSO and TSSCAD tend to select a large number of variables which can be found in all the scenarios except scenario 5 and scenario 6. For instance, TSLASSO selects 31458 variables and TSSCAD chooses 29113 variables in scenario 8. A large number of nuisance variables are included in the models given by TSLASSO and TSSCAD which will affect the forecasting results negatively. In terms of computational efficiency, TSSCAD spends more computational time in scenarios 1-2 and scenarios 5-6 when the dimensions are moderate. As the dimensions increase, all of the forecasting approaches take more computation cost. For instance, TSSRL took 1215.93 seconds and 1529.38 seconds in scenarios 4 and 8 when D=5000. TSLASSO runs fastest by using the convex penalty. However, notice that whenD is large, although TSSRL also uses the convex penalty, it is not computational efficient since the design of the TSSRL is based on COORD algorithm which runs much slower than coordinate descent algorithm. In scenario 9 and 10 where the sample sizes are large, TSSRHR still performs better than other methods. Same as other scenarios, TSSRL gives the most sparse model with fewer variables. Also notice that the train speed of TSSRHR is acceptable. Therefore, TSSRHR outperforms other forecasting methods in terms of forecasting accuracy, model sparsity, and computation speed.Table 1
The statistical performance of TSSRHR, TSLASSO, TSSCAD, and TSSRL. The performances are measured by number of nonzero elements of main effects, two-way interactions, three-way interactions, MSE (×100), RMSE (×100), and computational costs (CC, in seconds, unless otherwise specified).
Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC Scenario 1 Scenario 2  TSSRHR 3.00 1.00 3.00 93.47 96.68 10.62 TSSRHR 3.00 1.00 5.00 82.71 91.13 5.82 TSLASSO 5.00 59.50 3505.50 145.72 119.13 6.94 TSLASSO 7.00 46.00 3596.00 156.48 124.12 4.12 TSSCAD 2.00 1.00 8.00 127.40 112.92 18.77 TSSCAD 14.00 21.50 3678.00 185.49 132.81 20.73 TSSRL 3.00 0.00 1.00 102.71 101.36 5.26 TSSRL 2.00 0.00 4.00 86.78 91.98 8.13  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC Scenario 3 Scenario 4  TSSRHR 2.00 1.00 6.00 80.21 89.38 11.72 TSSRHR 3.00 0.00 0.00 87.92 93.13 20.46 TSLASSO 12.00 20.00 4796.00 137.54 117.13 4.32 TSLASSO 35.00 0.00 31458.00 264.13 161.81 11.12 TSSCAD 4.00 0.00 15.00 101.92 100.83 19.72 TSSCAD 19.00 0.00 6661.00 237.35 152.51 52.15 TSSRL 2.00 0.00 1.00 88.86 93.53 25.43 TSSRL 2.00 0.00 1.00 96.92 97.43 1215.93  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC Scenario 5 Scenario 6  TSSRHR 3.00 0.00 4.00 413.85 200.94 4.21 TSSRHR 3.00 0.00 4.00 418.72 205.19 4.92 TSLASSO 6.00 7.00 94.00 533.48 228.12 3.98 TSLASSO 5.00 2.00 50.50 438.59 207.14 4.01 TSSCAD 3.00 2.50 6.00 470.83 215.94 20.15 TSSCAD 6.00 1.00 94.00 473.21 216.83 25.13 TSSRL 3.00 0.00 1.00 468.10 215.18 2.13 TSSRL 2.00 0.00 3.00 428.13 206.89 3.72  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC Scenario 7 Scenario 8  TSSRHR 2.00 0.00 3.00 405.13 201.98 10.98 TSSRHR 3.00 0.00 2.00 407.95 201.34 28.13 TSLASSO 5.00 0.00 77.50 436.82 207.92 4.35 TSLASSO 16.00 0.00 1195.50 870.94 294.48 7.19 TSSCAD 9.00 28.00 451.50 434.12 206.15 30.95 TSSCAD 28.00 0.00 29113 454.79 211.83 72.49 TSSRL 2.00 0.00 1.00 430.93 205.73 3.43 TSSRL 2.00 0.00 1.00 445.23 209.89 1529.38  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC Scenario 9 Scenario 10  TSSRHR 3.00 1.00 4.00 147.33 121.38 14.03 TSSRHR 3.00 1.00 3.00 142.62 119.42 19.66 TSLASSO 3.00 1.00 12.00 194.07 139.31 4.09 TSLASSO 4.00 0.00 19.00 183.93 135.62 6.28 TSSCAD 3.00 4.00 18.00 195.77 139.92 8.73 TSSCAD 3.00 1.00 15.00 165.67 128.71 14.62 TSSRL 3.00 0.00 3.00 183.71 135.54 2.84 TSSRL 3.00 0.00 3.00 174.79 132.21 3.91Figure 1
Boxplots of compared forecasting methods in Scenarios 1-10.Figure 2
Sparsity pattern of compared forecasting methods in Scenarios 1-4.Figure 3
Sparsity pattern of compared forecasting methods in Scenarios 5-8.Figure 4
Sparsity pattern of compared forecasting methods in Scenarios 9-10.


### 5.2. Microarray Gene Expression Data

The inbred mouse microarray gene expression dataset with 31 female mice and 29 male mice is considered to demonstrate the effectiveness of the proposed TSSRHR method. The gene expression values of 22,690 genes are measured by each array and the response gene is a continuous phenotypic variable measured by stearoyl-coenzyme desaturase (scd1) with probe set ID 1415965 at. The remaining 22,689 probe set IDs are regarded as candidate variables. We are seeking a good model which is able to describe the relationship between the response variable and other candidate genes. The performances of forecast models on inbred mouse microarray gene expression data are shown in Table2. It worth mentioning that a gene filtering preprocessing procedure is applied to the inbred mouse data to reduce the expensive computational cost. TSSRHR delivers the best forecasting performances by giving the lowest MSE and RMSE followed by TSSRL, TSLASSO, and TSSCAD. TSLASSO selects a model with 6311 three-way interaction terms while TSSCAD only selects 3 variables which is fewer than those of TSSRL and TSSRHR. This can be verified in Figure 5 which shows the selection pattern of the compared methods. The computational cost of TSSRL is almost 4 times of TSSRHR. Comparing with TSSRHR and TSSCAD, TSLASSO and TSCAD spend less computational cost. Overall, TSSRHR provides the lowest forecasting error with sufficient model sparsity and acceptable computational time.Table 2
The statistical performances of TSSRHR, TSLASSO, TSSCAD, and TSSRL. The performances are measured by number of nonzero elements of main effects, two-way interaction, three-way interactions, three-way interactions, MSE (×100), RMSE (×100), and computational costs (CC, in seconds, unless otherwise specified).
Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC TSSRHR 1.00 1.00 7.00 8.61 29.33 0.48 TSLASSO 1.00 63.00 6311.00 136.66 116.89 0.22 TSSCAD 1.00 1.00 1.00 142.75 119.48 0.24 TSSRL 0.00 0.00 14.00 10.98 33.13 1.50Figure 5
Sparsity pattern of compared forecasting methods in the inbred mouse gene expression data.


## 5.1. Twinsine Signal Data

Spectral estimation studying the distribution over frequencies is widely applied in speech coding and radar solar signal processing. Nevertheless, it is critically difficult and challenging when there are a large number of frequencies which causes high resolution problem. This type of problem is called super resolution spectral estimation. To avoid the power leakage and noisy observation caused by the super resolution, a novel method is urgently needed to choose the important frequencies which plays a dominant role in super-resolution estimation. The twinsine signal is generated from the target function:(22)yt=A1cos⁡2πf1t+ϕ1+A2cos⁡2πf2t+ϕ2+εtwhere A1=2, A2=3, ϕ1=π/3, f1=0.25HZ, and ε(t) follows white Gaussian noise with variance given by σ2. The frequence resolution is chosen as 0.002HZ to detect and distinguish the two sinusoidal components. fk=fmax·k/D for k=0,1,⋯,D and the frequency atoms are given by Xcos=cos(2πtifk)1≤i≤n,1≤k≤D and Xsin=sin(2πtifk)1≤i≤n,1≤k≤D. In our experiments, we will consider the following 8 different scenarios:(i)
Scenario 1:(n,D,σ2)=(100,300,1)(ii)
Scenario 2:(n,D,σ2)=(100,500,1)(iii)
Scenario 3:(n,D,σ2)=(100,1000,1)(iv)
Scenario 4:(n,D,σ2)=(100,5000,1)(v)
Scenario 5:(n,D,σ2)=(100,300,5)(vi)
Scenario 6:(n,D,σ2)=(100,500,5)(vii)
Scenario 7:(n,D,σ2)=(100,1000,5)(viii)
Scenario 8:(n,D,σ2)=(100,5000,5)(ix)
Scenario 9:(n,D,σ2)=(500,300,1)(x)
Scenario 10:(n,D,σ2)=(1000,300,1)The results of forecasting, variable selection, and computation are summarized in Table1 and Figures 1, 2, 3, and 4. It is not difficult to observe that TSSRHR obtains the most accurate outcome in all the scenarios including the most challenging scenario when D=5000 and σ2=5. TSSRL outperforms TSLASSO and TSSCAD method by delivering lower MSE and its forecasting results are comparable with TSSRHR in scenario 6. From the prospective of model sparsity, TSSRL achieves the most interpretable model followed by TSSRHR. TSLASSO and TSSCAD tend to select a large number of variables which can be found in all the scenarios except scenario 5 and scenario 6. For instance, TSLASSO selects 31458 variables and TSSCAD chooses 29113 variables in scenario 8. A large number of nuisance variables are included in the models given by TSLASSO and TSSCAD which will affect the forecasting results negatively. In terms of computational efficiency, TSSCAD spends more computational time in scenarios 1-2 and scenarios 5-6 when the dimensions are moderate. As the dimensions increase, all of the forecasting approaches take more computation cost. For instance, TSSRL took 1215.93 seconds and 1529.38 seconds in scenarios 4 and 8 when D=5000. TSLASSO runs fastest by using the convex penalty. However, notice that whenD is large, although TSSRL also uses the convex penalty, it is not computational efficient since the design of the TSSRL is based on COORD algorithm which runs much slower than coordinate descent algorithm. In scenario 9 and 10 where the sample sizes are large, TSSRHR still performs better than other methods. Same as other scenarios, TSSRL gives the most sparse model with fewer variables. Also notice that the train speed of TSSRHR is acceptable. Therefore, TSSRHR outperforms other forecasting methods in terms of forecasting accuracy, model sparsity, and computation speed.Table 1
The statistical performance of TSSRHR, TSLASSO, TSSCAD, and TSSRL. The performances are measured by number of nonzero elements of main effects, two-way interactions, three-way interactions, MSE (×100), RMSE (×100), and computational costs (CC, in seconds, unless otherwise specified).
Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC Scenario 1 Scenario 2  TSSRHR 3.00 1.00 3.00 93.47 96.68 10.62 TSSRHR 3.00 1.00 5.00 82.71 91.13 5.82 TSLASSO 5.00 59.50 3505.50 145.72 119.13 6.94 TSLASSO 7.00 46.00 3596.00 156.48 124.12 4.12 TSSCAD 2.00 1.00 8.00 127.40 112.92 18.77 TSSCAD 14.00 21.50 3678.00 185.49 132.81 20.73 TSSRL 3.00 0.00 1.00 102.71 101.36 5.26 TSSRL 2.00 0.00 4.00 86.78 91.98 8.13  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC Scenario 3 Scenario 4  TSSRHR 2.00 1.00 6.00 80.21 89.38 11.72 TSSRHR 3.00 0.00 0.00 87.92 93.13 20.46 TSLASSO 12.00 20.00 4796.00 137.54 117.13 4.32 TSLASSO 35.00 0.00 31458.00 264.13 161.81 11.12 TSSCAD 4.00 0.00 15.00 101.92 100.83 19.72 TSSCAD 19.00 0.00 6661.00 237.35 152.51 52.15 TSSRL 2.00 0.00 1.00 88.86 93.53 25.43 TSSRL 2.00 0.00 1.00 96.92 97.43 1215.93  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC Scenario 5 Scenario 6  TSSRHR 3.00 0.00 4.00 413.85 200.94 4.21 TSSRHR 3.00 0.00 4.00 418.72 205.19 4.92 TSLASSO 6.00 7.00 94.00 533.48 228.12 3.98 TSLASSO 5.00 2.00 50.50 438.59 207.14 4.01 TSSCAD 3.00 2.50 6.00 470.83 215.94 20.15 TSSCAD 6.00 1.00 94.00 473.21 216.83 25.13 TSSRL 3.00 0.00 1.00 468.10 215.18 2.13 TSSRL 2.00 0.00 3.00 428.13 206.89 3.72  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC Scenario 7 Scenario 8  TSSRHR 2.00 0.00 3.00 405.13 201.98 10.98 TSSRHR 3.00 0.00 2.00 407.95 201.34 28.13 TSLASSO 5.00 0.00 77.50 436.82 207.92 4.35 TSLASSO 16.00 0.00 1195.50 870.94 294.48 7.19 TSSCAD 9.00 28.00 451.50 434.12 206.15 30.95 TSSCAD 28.00 0.00 29113 454.79 211.83 72.49 TSSRL 2.00 0.00 1.00 430.93 205.73 3.43 TSSRL 2.00 0.00 1.00 445.23 209.89 1529.38  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC  Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC Scenario 9 Scenario 10  TSSRHR 3.00 1.00 4.00 147.33 121.38 14.03 TSSRHR 3.00 1.00 3.00 142.62 119.42 19.66 TSLASSO 3.00 1.00 12.00 194.07 139.31 4.09 TSLASSO 4.00 0.00 19.00 183.93 135.62 6.28 TSSCAD 3.00 4.00 18.00 195.77 139.92 8.73 TSSCAD 3.00 1.00 15.00 165.67 128.71 14.62 TSSRL 3.00 0.00 3.00 183.71 135.54 2.84 TSSRL 3.00 0.00 3.00 174.79 132.21 3.91Figure 1
Boxplots of compared forecasting methods in Scenarios 1-10.Figure 2
Sparsity pattern of compared forecasting methods in Scenarios 1-4.Figure 3
Sparsity pattern of compared forecasting methods in Scenarios 5-8.Figure 4
Sparsity pattern of compared forecasting methods in Scenarios 9-10.


## 5.2. Microarray Gene Expression Data

The inbred mouse microarray gene expression dataset with 31 female mice and 29 male mice is considered to demonstrate the effectiveness of the proposed TSSRHR method. The gene expression values of 22,690 genes are measured by each array and the response gene is a continuous phenotypic variable measured by stearoyl-coenzyme desaturase (scd1) with probe set ID 1415965 at. The remaining 22,689 probe set IDs are regarded as candidate variables. We are seeking a good model which is able to describe the relationship between the response variable and other candidate genes. The performances of forecast models on inbred mouse microarray gene expression data are shown in Table2. It worth mentioning that a gene filtering preprocessing procedure is applied to the inbred mouse data to reduce the expensive computational cost. TSSRHR delivers the best forecasting performances by giving the lowest MSE and RMSE followed by TSSRL, TSLASSO, and TSSCAD. TSLASSO selects a model with 6311 three-way interaction terms while TSSCAD only selects 3 variables which is fewer than those of TSSRL and TSSRHR. This can be verified in Figure 5 which shows the selection pattern of the compared methods. The computational cost of TSSRL is almost 4 times of TSSRHR. Comparing with TSSRHR and TSSCAD, TSLASSO and TSCAD spend less computational cost. Overall, TSSRHR provides the lowest forecasting error with sufficient model sparsity and acceptable computational time.Table 2
The statistical performances of TSSRHR, TSLASSO, TSSCAD, and TSSRL. The performances are measured by number of nonzero elements of main effects, two-way interaction, three-way interactions, three-way interactions, MSE (×100), RMSE (×100), and computational costs (CC, in seconds, unless otherwise specified).
Nz(b) Nz(Φ) Nz(Ψ) MSE RMSE CC TSSRHR 1.00 1.00 7.00 8.61 29.33 0.48 TSLASSO 1.00 63.00 6311.00 136.66 116.89 0.22 TSSCAD 1.00 1.00 1.00 142.75 119.48 0.24 TSSRL 0.00 0.00 14.00 10.98 33.13 1.50Figure 5
Sparsity pattern of compared forecasting methods in the inbred mouse gene expression data.


## 6. Conclusion

Two-way interaction model has been widely used in numbers of scientific fields. However, researchers find that three-way interaction terms play more important role than two-way interaction terms when establishing forecasting model. To this end, this paper studies high-order interaction model which has both two-way and three-way interaction terms. To establish an interpretable model and reduce the computation time, sparse representation method using a two-stage square root hard ridge (TSSRHR) is proposed and studied. Based on the property of two-stage method, the number of interactions will decrease dramatically so that the computational issue is solved. To implement SRHR in each stage, a simple and efficient algorithm is designed. Furthermore, the prediction and estimation error bounds of TSSRHR are also exhibited using two regularity assumptions. For real data analysis, a twinsine signal data and microarray gene expression data are provided to show the forecasting and selection performances of TSSRHR comparing with other stagewise sparse representation methods. A hierarchical variable selection problem on the high-order interaction model will be conducted for the future research.



---

*Source: 2032987-2018-11-11.xml*