Loann D Desboulets
email: loann.desboulets@univ-amu.fr

A Review on Variable Selection in Regression Analysis

Keywords: Variable Selection, Automatic Modelling, Sparse Models JEL Classification: C50,C59

In this paper, we investigate on 39 Variable Selection procedures to give an overview of the existing literature for practitioners. "Let the data speak for themselves" has become the motto of many applied researchers since the amount of data has significantly grew. Automatic model selection have been raised by the search for data-driven theories for quite a long time now. However while great extensions have been made on the theoretical side still basic procedures are used in most empirical work, eg. Stepwise Regression. Some reviews are already available in the literature for variable selection, but always focus on a specific topic like linear regression, groups of variables or smoothly varying coefficients. Here we provide a review of main methods and state-of-the art extensions as well as a topology of them over a wide range of model structures (linear, grouped, additive, partially linear and non-parametric). We provide explanations for which methods to use for different model purposes and what are key differences among them. We also review two methods for improving variable selection in the general sense.

Introduction

When building a statistical model the question of which variables to include arise very often. In practice is it true almost all the time. This can come from ignorance, competing theories, or whatever. Practitioners have now at their disposal a wide range of technologies to solve this issue. Literature on this topic started with Stepwise Regression [START_REF] Breaux | On stepwise multiple linear regression[END_REF] and Autometrics [START_REF] Hendry | Recent developments in the theory of encompassing[END_REF], moving to more advanced procedures from which the most famous are the Non Negative Garrotte [START_REF] Breiman | Better subset regression using the nonnegative garrote[END_REF], the Least Angle and Shrinkage Selection Operator (hereafter LASSO, [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]) and the Sure Independence Screening [START_REF] Fan | Statistical methods with varying coefficient models[END_REF].

Many papers are available for empiricists to get an overview of the existing methods. [START_REF] Fan | A selective overview of variable selection in high dimensional feature space[END_REF] reviews most of the literature on linear and generalized models. A large part is devoted to penalized methods and algorithmic solutions, also the optimal choice of the parameter penalty is discussed. [START_REF] Breheny | Penalized methods for bi-level variable selection[END_REF] and [START_REF] Huang | A selective review of group selection in high-dimensional models[END_REF] gave a complete review of selection procedures in grouped variables models with great technical comparisons, especially in terms of rate of convergence. [START_REF] Castle | Evaluating automatic model selection[END_REF] compared Autometrics to a wide range of other methods (Stepwise, Akaike Information Criterion, LASSO, etc. 1 ) in terms of prediction accuracy under orthogonality of the regressors, with a particular attention given to dynamic models. In the same spirit as our paper [START_REF] Park | Varying coefficient regression models: a review and new developments[END_REF] gave a very recent review of variable selection procedures but dealing only with varying-coefficient models. [START_REF] Fan | Sure independence screening[END_REF] provided a comprehensive review in the context of Sure Independence Screening major improvements. We can also cite more focusing papers like [START_REF] Fu | Penalized regressions: the bridge versus the lasso[END_REF] who compared the Bridge and the LASSO theoretically but also empirically both through simulation and real data. [START_REF] Epprecht | Comparing variable selection techniques for linear regression: Lasso and autometrics[END_REF] that compared Autometrics and the LASSO according to prediction and selection accuracy.

The contribution of this paper is threefold. First, 39 procedures are considered, these are listed and clearly classified. Secondly, we establish a topology of procedures under different model structures. We consider major ones: Linear Models, Grouped Variables Models, Additive Models, Partial Linear Models and Non-parametric Models. Thirdly, we describe and compare state-of-the-art papers in the literature. We give contexts where each procedure should be used, to which specific problem they answer and compare them on this ground. In this sense any practitioner with enough knowledge in Statistics can refer to our paper as a methodological guide for doing variable selection. It gives a wider view of existing technologies than the other reviews we mentioned.

The paper is organized as follows, in section 2 we introduce the three main categories of variable selection procedures and we provide a typological table of these ones on the ground of model structures. Descriptions as well as comparisons are discussed along the sections 3 to 7. Each of these sections focuses on a particular model structure. Section 8 is devoted to two general methods for improving model selection, both can be applied on all the procedures presented accross the paper. In the final section we make few critics on actual procedures and give insight on future area of research.

Typology of Procedures

In this section we propose a typology of state-of-the-art selection methods in many different frameworks, there are many types of models that can be considered. For this aim, Table 1 provides the classification of statistical methods available in the literature and that will be discussed in the paper. From the latter we have determined 3 main categories of algorithms: Then there are Penalty-based procedures. Imposing a constraint on parameters directly inside estimation encourages sparsity among them. For instance LASSO and Ridge belongs to this category.

The lasts are Screening procedures, they are not all designed to do selection intrinsically but rather ranking variables by importance. The main advantage is that it applies more easily to very large dimensional problems, when number of regressors is diverging with the number of observations (eg. cases with p >> n). This is mainly true because it considers additive models (linear or not) and so variables can be treated separately.

We discuss this distinction more deeply in subsections below and give brief description of their main features.

Tests-Based

This is the most classical way of handling variable selection in statistical models. It was also the first attempt of variable selection. Everything started with Stepwise Regression [START_REF] Breaux | On stepwise multiple linear regression[END_REF], one of the latest of this kind is Autometrics [START_REF] Hendry | Recent developments in the theory of encompassing[END_REF] 2 . We focus on Stepwise Regression and Autometrics for two reasons. The first is that Stepwise Regression is the most well-known and the most widespread method for choosing variables in a model. Despite it dates back to 1967 many empiricists still practice it. The second is that Autometrics has integrated many features of Econometrics to achieve the highest degree of completeness for an automatic procedure. Authors have considered endogeneity, non-linearities, unit-roots and many others, trying to overcome most issues a statistician can face.

Stepwise Regression is the most simple and most straightforward way of doing model selection by just retrieving 2 Even though it started in 1987 there are still improvements nowadays. insignificant variables (backward approach) or adding significant ones (forward approach) based on some statistical criterion. Therefore it is pretty easy to use it empirically because implementation is straightforward.

However in several situations this does not ensure consistent selection. Its selection properties have been investigated in [START_REF] Wang | Forward regression for ultra-high dimensional variable screening[END_REF]. On the contrary Autometrics is a complete philosophy of modelling, but comes at the cost of a quite complex algorithm and many tuning parameters are required, making its use more difficult for non-expert.

Penalty-Based

Thanks to the work of [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] it became a quite common strategy in empirics. This kind of methods involves applying a penalty on parameters to encourage sparsity (i.e. some are set exactly to zero).

Sparsity is a necessary condition for situations of unidentifiability ie. where p > n. Such a problem can be solved using penalties on parameters to make inference possible. These parameters can come from parametric models or from non-parametric models, so penalty based method can be applied on both structures. This kind of procedure started with the Non Negative Garrote (NNG of [START_REF] Breiman | Better subset regression using the nonnegative garrote[END_REF]) in an Ordinary Least Squares framework up now to much complex model structure like varying coefficients and two-ways interactions ANOVA non-parametric models. The idea of producing sparse models is a convenient way of integrating a test inside the estimation. Inference of such models requires the prior assumption that some variables are not relevant, this is the test part, and penalty-based methods helps estimating the coefficients, this is the inference part. So both procedures are merged in an unified framework giving rise to a novel conception of statistical modelling. Maybe the most famous in this category is the LASSO of [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF].

Screening-Based

Screening is actually the most effective way of dealing with very high dimensional features (large p). Few other selection methods can be as computationally efficient as these ones. However Screening often does not perform model selection itself, it rather ranks variables. To do so they have to be mixed up with other procedures, in the literature they are mainly penalty-based. Even if it does not select variables reducing the candidate set is an important aspect of the variable selection and screening methods are powerful in this task. In this respect it is worth mentioning the Sure Independence Screening (hereafter SIS, [START_REF] Fan | Sure independence screening for ultrahigh dimensional feature space[END_REF]) that is the first of this kind.

Screening makes the use of a ranking measure, either linear or not so it can be applied in both frameworks. Some may rely on specific models (like a linear model) while others are model-free. The major differences among procedures in this category relies on the choice of the ranking measure. Correlation coefficients are the first coming to mind, these are mainly used. One limitation in screening is that they usually treats variables by pairs to compute their measure of association, so every effect is considered as additive and does not correct for the presence of interactions effects. This is not necessarily true, especially in the non-parametric settings. Sophisticated correlations such as distance correlation or canonical kernel correlation are employed in a multivariate framework and account for such interactions even if they do not model them explicitly. However in this case they may loose their computational efficiency compared to independence screening ones. As said before, a brief review of some SIS methods can be found in [START_REF] Fan | Sure independence screening[END_REF].

Linear Models

We began with the first model structure: the linear model. It is described as:

y = Xβ + ε (1)
The variable to be explained y (sometimes also called the output, the response or the dependent variable) is a one dimensional vector of length n, corresponding to the number of observations. The matrix X contains the explanatory variables (sometimes also called the inputs, the regressors or the independent variables) of length n and dimension p which is the number of candidate variables. Therefore the one dimensional vector β of length k contains the parameters of interest. The residuals (sometimes called also the error term) are denoted ε, even if it could be of interest we do not solely focus on their properties and consequences on variable selection in this paper. This notation will be held constant throughout the paper. Notice that all of the three methodologies are able to handle linear models, while this is not necessarily true for other structures (e.g. additive models).

Testing

Stepwise Regression [START_REF] Breaux | On stepwise multiple linear regression[END_REF]) is the first model selection procedure. This approach have been motivated when statisticians started to consider model uncertainty. This means that among p variables we can possibly construct 2 p models, so we should maybe take them all into account. To test all possibilities we have to compute "all-subsets". This cannot be achieved for large p. In order to overcome this problem and reduce the search, stepwise regression investigates only a subset of all possible regressions with the hope to end with the true model. There exist two approaches: Backward and Forward. Either the process starts from a null model (only an intercept) and introduces variables one by one, this is the forward step. Or it starts from the full model (all variables) and deletes them one by one, this is the backward step. One improvement is also to consider both. Usually the selection within each step is made according to some criterion. One consider all one-variable increments from the actual model and choose the best move according to this criterion, it might be the lowest p-value, highest adjusted R 2 , lowest Mallow's C p , lowest AIC, lowest prediction error, leave-one-out cross validation, etc.

You can imagine any criterion to perform this job, but the main issue arising from Stepwise Regression does not come from the choice of the criterion. Interesting critics [START_REF] Doornik | Econometric model selection with more variables than observations[END_REF] arise from the developers of Autometrics.

The main one is the lack of search. Stepwise regression proceeds step by step along a single path. Then, there is no backtesting. That is the procedure never considers testing again variables that have been removed after each step. Such an idea is present using the forward backward combination but it is restricted to the previous step only.

Obviously they are not the only one to express admonitions about Stepwise Regression. We can mention many papers from [START_REF] Hurvich | The impact of model selection on inference in linear regression[END_REF], [START_REF] Steyerberg | Stepwise selection in small data sets: a simulation study of bias in logistic regression analysis[END_REF], [START_REF] Whittingham | Why do we still use stepwise modelling in ecology and behaviour[END_REF] or [START_REF] Flom | Stopping stepwise: Why stepwise and similar selection methods are bad, and what you should use[END_REF] where they all prove biased estimation and inconsistent selection of Stepwise Regression.

However even if used as a selection method it behaves poorly, used a screening method it showed better results.

This has been developed by [START_REF] Wang | Forward regression for ultra-high dimensional variable screening[END_REF] and is detailed in the next subsection.

On the other side is Autometrics, an algorithm for model selection developed by [START_REF] Hendry | Recent developments in the theory of encompassing[END_REF] under the famous Theory of Encompassing and the LSE (London School of Economics) type of Econometrics.

This method has been created as early as 1987 and is still under development. The basis of its methodology is the General-to-Specific approach. Theory of Encompassing states that the researcher should start from a very large model (called the GUM: General Unrestricted Model) encompassing all other possible models and then reduce it to a simpler but congruent specification. This idea is somehow related to the backward specification in

Stepwise Regression. His work is an automation of standard way of testing for relevance in Econometrics such as t-tests and F-tests and major concerns deal with power of tests, repeated testing but also outliers detection, non-linearities, high dimensional features (with p > n) and parameter invariance.

Tests come with some hyperparameters specifying the size of the battery of tests (t-tests, F-tests, normality checks, etc.).

Repeated testing occurs when a variable that has been deleted under a certain specification that has now changed is reintroduced and tested again. The absence of such a thing in Stepwise Regression is a severe drawback and the main reason why it fails pretty often.

Non-linearities are handled using Principal Components Analysis (see [START_REF] Castle | A low-dimension portmanteau test for non-linearity[END_REF]) that makes the design matrix orthogonal. Such a decomposition allows to introduce squares and cubics of the transformed variables which are linear combination of the original ones. Orthogonality limits the number of non-linear terms since it already accounts for interaction using components. In simple words a polynomial of degree d with p variable results in ( d+p d ) -1 terms, while their methods reduces to d × p which is very much less. It is advocated that it can reproduce non-linear functions often met in Economics and Social Sciences. However the class of functions that it can reproduce may be restricted compared to standard non-parametric methods 3 .

High dimensional features and non-identifiability (p > n) of the GUM is solved in a very simple way called "Block Search". Regressors are divided in different blocks until the size of each block is lower than p. Then tests are applied in each block, some variables are discarded, the remaining blocks are merged and the process continues. This idea is based on the fact that the methodology is still consistent under separability. This idea is quite similar to the Split-and-Conquer methodology of Chen and Xie (2014) to solve ultra-high dimensional problems.

Outliers can be detected using the Impulse Indicator Saturated Selection (IIS) developed by [START_REF] Hendry | Selecting a regression saturated by indicators[END_REF].

This is in the same spirit as the Block-Search (or Split-and-Conquer) approach defined previously. A set of indicator is added to the GUM for every observation, and tests are applied in a Block-Search manner to remove observations that are not consistent with the model, identified as outliers.

3

This should be investigated more deeply, to the best of our knowledge no papers have tried to compare their non-linear regression to the very well-known non-parametric procedures like Kernels or Splines. An obvious link can be made with Projection Pursuit Regression (PPR), in this respect we claim that Autometrics may be a special case of PPR.

Stepwise Regression and Autometrics are serial procedures where selection and estimation are performed sequentially. In some sense Penalty-Based methods aim at performing both at the same time. One can view penalty-based procedures as the direct implementation of tests inside inference.

Penalty

Penalty-based methods can be divided in two categories: penalties on the norm and concave ones. The shape of the penalty may have a great influence on the selected set of variables and their estimates. Sparse model is achieved because we reduce nearly zero coefficients to zero in estimation. The penalty parameter plays the role of a threshold but in a non-orthogonal framework. To understand better the origins of these penalty one should refer to threshold methods in [START_REF] Kowalski | Thresholding rules and iterative shrinkage/thresholding algorithm: A convergence study[END_REF]. For that reason the penalty also introduces shrinkage of the coefficients, making them biased. The literature is focused on the choice of the penalty in terms of selection consistency and bias properties.

Norm Penalties

There are almost as many methods as there are norms, but generally the objective is to solve:

min β y -Xβ 2 2 + λ β 2 γ (2)
Each methods applies to different L γ norms.

• SparseStep: γ = 0

• LASSO : γ = 1 • Ridge : γ = 2
This methodology is gathered in the more general Bridge estimator [START_REF] Frank | A statistical view of some chemometrics regression tools[END_REF]) that considers any value for γ, but the authors did not say how to solve the problem. The advantage of Ridge [START_REF] Hoerl | Ridge regression: Biased estimation for nonorthogonal problems[END_REF] is that it has an analytical solution. However the solution is not sparse so it does not select variables (only shrinkage). The Least Absolute Shrinkage and Selection Operator (LASSO, [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]) does because the L 1 norm is singular at the origin. However both give bias estimates because they apply shrinkage to the coefficients. The zero norm used in SparseStep (van den Burg et al. 2017) is the counting norm, they penalize directly the number of non-zero elements in β , not their values (no shrinkage). Usually constraints on the number of non-zero elements require high computational costs (exhaustive search over the model space). Here they use an easy even though very precise continuous approximation from de Rooi and Eilers (2011) and that turns the problem into something computationally tractable. [START_REF] Meinshausen | High-dimensional graphs and variable selection with the lasso[END_REF] shown that LASSO tends to select noise variables using a penalty parameter optimally chosen for prediction. For this reason [START_REF] Zou | The adaptive lasso and its oracle properties[END_REF] developed AdaLASSO (Adaptive LASSO). His paper proved that the optimal estimation rate is not compatible with consistent selection. Moreover even sacrificing the estimation rate does not ensure that the LASSO will select the right variables with positive probability. This phenomenon is highlighted through a necessary condition on the covariance matrix of the regressors that cannot always be satisfied using the LASSO with a single penalty parameter. Therefore he introduced adaptive weights to the LASSO to make it consistent with variable selection.

min β y -Xβ 2 2 + λ wβ 1 (3)
The latest improvement on linear models is to allow for interactions terms. Even if it is possible, only adding them into a LASSO is not an efficient procedure because it greatly extends the dimensionality of the design matrix. The idea of the Strong Heredity Interaction Model (SHIM, [START_REF] Choi | Variable selection with the strong heredity constraint and its oracle property[END_REF]) is to add interactions only if main effects are selected also (strong heredity property), this greatly reduces the search space and provides an efficient way of doing ANOVA-types of models. They consider a reparametrization of the two-ways interactions models:

y = Xβ + p ∑ j=1 ∑ k = j γ jk β j β k x j x k (4)
Introducing main effect parameters β on top of cross-effects γ ensures that the interaction will be non-zero if and only if both main effects are non-zeros. The problem is a composite LASSO of the following form:

min β ,γ y -Xβ 2 2 + λ β β 1 + λ γ γ 1 (5)
Solutions to these problems are numerous. Usually either it reduces to the LASSO and then algorithms like Least Angle Regression (LARS) of [START_REF] Efron | Least angle regression[END_REF] are employed. Otherwise iterative algorithms like the Local Quadratic Approximation Fan and Li ( 2001) can be used.

Concave Penalties

Norm penalties are very standard and easy to work with but there exists also other types of penalties. Thus we can consider penalties in a very general framework:

min β y -Xβ 2 2 + p λ (β ) (6)
The difference will then lie in the choice of p λ (β ).

• NonNegativeGarotte:

p λ (β ) = nλ p ∑ j=1 1 - λ β 2 j + (7)
• SCAD :

p λ (β ) =        λ , if |β | ≤ λ aλ -|β | a-1 , if λ < |β | < aλ 0, if |β | ≥ aλ (8)
• MCP :

p λ (β ) = λ |β | -x 2 2γ , if |β | ≤ γλ 0.5γλ 2 , if |β | > γλ (9)
The Non Negative Garotte [START_REF] Breiman | Better subset regression using the nonnegative garrote[END_REF] was the first penalty of this kind, but because it has bad properties (especially variables selection inconsistency) it was rapidly abandoned. SCAD (Smoothly Clipped Absolute Deviation, [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF]) was the first penalty method that was consistent, continuous and unbiased for large values of β . MCP (Minimax Convex penalty,Zhang et al. ( 2010)) has little difference with SCAD in terms of selected variables. A comparative study between them can be found in [START_REF] Zhang | Penalized linear unbiased selection[END_REF].

One thing with penalty method is that there are always some penalty parameters (eg λ in LASSO) that have to be chosen. Usually they are set to optimal values according to some General Cross Validation (GCV) criterion or out-of-sample predictions. This is crucial because results can be very sensitive to the choice of these parameters.

SCAD is more robust to this problem thanks to a bias-free property. 4

Screening

Another methodology in variable selection is Screening. In fact these are ranking methods that rely on some association measure between the dependent variable and the regressors. Very often this measure is taken to be bivariate allowing then an extremely fast analysis.

4

This is true only for large values of parameters, the reader can get intuitions of this phenomenon with threshold methods [START_REF] Kowalski | Thresholding rules and iterative shrinkage/thresholding algorithm: A convergence study[END_REF].

Regressor Based

The Sure Independence Screening (SIS, [START_REF] Fan | Sure independence screening for ultrahigh dimensional feature space[END_REF]) is the first of this kind and almost all methods are derived from it. It uses simple correlation on standardized variables : ω(x j , y) = x j ỹ and gives a ranking of the x j . The set M of relevant features is determined by a simple threshold:

M = {1 ≤ j ≤ p : | ω(x j , y)| is among the top d largest ones} ( 10
)
This set is reduced step by step until some moment. The method in itself does not select anything in fact, it just remove the less correlated features from the set of candidates, but we are left with a candidate set where selection has to apply. SIS needs a selection procedure in the end to obtain consistent results. The main advantage of the method is that when the number of variables p is very large compared to the number of observations n usual selection procedures tend to misbehave [START_REF] Fan | Sure independence screening for ultrahigh dimensional feature space[END_REF]. In their paper SIS has proven to lead to a set of candidates that is manageable for LASSO and others in order to have good properties. SIS allows for ultrahigh dimensional features, ultrahigh being defined as:

log(p) = O(n α ) with 0 ≤ α ≤ 1.
In this respect the screening properties of screening of Forward Regression [START_REF] Wang | Forward regression for ultra-high dimensional variable screening[END_REF]) have been investigated

and with little improvements proved to be consistent in variables selection. However it still requires a selection procedure in the end, Forward Regression is just used for the screening part that is ranking and reducing the set of candidates.

Because SIS may encounters issue for selecting weakly correlated variables (weak signal-to-noise ratio) [START_REF] Fan | Sure independence screening for ultrahigh dimensional feature space[END_REF] introduced Iterative conditional SIS that is applying correlation ranking but conditional on selected features. This is equivalent as looking through correlation between features and residuals from a model using primarily selected variables instead of correlation with the dependent variable. This idea can be related to former algorithms that were developed to infer the LASSO (eg. Forward Stagewise).

Covariance Based

The last approach is less common. The Covariate Assisted Screening Estimates (CASE, [START_REF] Ke | Covariate assisted screening and estimation[END_REF]) is a method that looks for sparse models but in the case where signals are rare and weak. All methods presented so far work well if β is sparse (so rare) and has high values (strong signals). In this case methods like SCAD are even bias-free. But if the signals are weak on top of rare then they won't manage to perform variable selection very well. The idea in CASE is to sparsify the covariance matrix of the regressors using a linear filter and then look for models inside this sparse covariance matrix using tests and penalties. Drawbacks are the choice of the filter that is problem dependent and the power of the tests.

To improve on CASE when regressors are highly correlated, giving a very dense covariance structure, Factor Adjusted-Covariate Assisted Ranking (FA-CAR, Ke and Yang (2017)) proposes using PCA to sparsify it. This is in line with selecting appropriately the filter in CASE when the problem to solve includes strong collinearity. In fact the covariance is assumed to have a sparse structure, hidden by latent variables. These are estimated by PCA and then removed from the variables. The process does not change anything for the equation and the parameters to be estimated does not require more technology than the simple OLS on the transformed decorrelated variables.

The main issue is to select the number of latent variables to be removed, this can be done via cross-validation for instance, still it remains difficult.

Grouped Models

Depending on the application the model can come in a group structure form of the type:

y = G ∑ g=1 ∑ j∈g β j x j + ε (11)
which can be rewritten in matrix-grouped notation:

y = G ∑ g=1 X g β g + ε (12)
Within this framework there are 2 main possibilities. One can look for which group to be selected or which variable is more relevant in which group. The former is referred to as single-level selection (sparse between group estimates) and the latter as bi-level selection (sparse between and within group estimates). Technical reviews of selection procedures with grouped variables can be found in [START_REF] Breheny | Penalized methods for bi-level variable selection[END_REF] and [START_REF] Huang | A selective review of group selection in high-dimensional models[END_REF].

Penalty

Single-level

The concept of group-penalty was introduced in Yuan and Lin (2006) (groupLASSO) in a LASSO framework. The objective is to solve a modified LASSO:

min β y - G ∑ g=1 ∑ j∈g β j x j 2 2 + λ G ∑ g=1 c g (β g R g β g ) 1/2 (13)
The parameters c g are used to adjust for the group sizes in order to have selection consistency. The parameter λ controls for the penalty. The choice of R g that weights each coefficients within the group is still challenging.

A solution is to take R g = (X g X g )/n the Gram matrix of the grouped variables X g . The effect is to scale the variables within groups and so make coefficients comparable in some sense. It can be easily shown that this lead to the LASSO solution with standardization of regressors when the group is formed with only one variable, such a thing is made pretty often empirically and is even advised by the LASSO's authors. An obvious extension is to take into account any penalty, providing the following objective:

min β y - G ∑ g=1 ∑ j∈g β j x j 2 2 + p( G ∑ g=1 β g R g ; c g λ , γ) (14)
Where p(•) can be taken to be the Bridge, the SCAD or the MCP criterion introducing then the groupBridge [START_REF] Huang | A group bridge approach for variable selection[END_REF]), the groupSCAD [START_REF] Wang | Group scad regression analysis for microarray time course gene expression data[END_REF] and the groupMCP [START_REF] Breheny | Penalized methods for bi-level variable selection[END_REF] respectively.

Bi-level

Improvements have been made on norm penalties by considering mixed norms like the ElasticNet [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]: min

β y -Xβ 2 2 + λ 1 β 1 + λ 2 β 2 2 ( 15
)
This method overcomes the issue of collinearity because it favours selection of correlated regressors simultaneously while LASSO tends to select only on out of them. In fact the EasticNet can be solved as a LASSO using slight modification of the LARS algorithm. Since it is a mix of Ridge and LASSO, parameters can be estimated by Ridge in a first step then apply the LASSO. A small correction due to the second penalty λ 2 is required. Originally the Elastic-net was not designed explicitly for grouped structure.

Also composite penalties have been considered in [START_REF] Breheny | Penalized methods for bi-level variable selection[END_REF] using the MCP criterion at both stages (between and within).

Since there is a great literature of reviews on these method (Breheny andHuang 2009, Huang et al. 2012) we do not spend time giving more details and advise readers interested in group models to have a look at them.

Additive Models

A step further in model structure complexity is to consider different non-parametric functions associated with each variables. The non-parametric additive model takes the following form:

y = p ∑ j=1 f j (x j ) + ε (16)

Penalty

The Sparse Additive Model (SpAM) of [START_REF] Ravikumar | Spam: Sparse additive models[END_REF] applies to this kind of models. The idea is simply to apply the LASSO to functions non-parametrically fitted with parametric coefficients coming in top of them. This is obviously the most natural extension of LASSO to the additive structure. The main program to solve is:

min β , f j y - p ∑ j=1 β j f j (x j ) 2 2 + λ p ∑ j=1 |β j | (17)
Even though the term ∑ p j=1 β j f j (x j ) remind us the very well-known Splines where the f j would be the basis functions, the authors claim that any non-parametric method can be used for fitting them. The solution is given in the form of a backfitting algorithm (Breiman and Friedman 1985). Another approach have been investigated by [START_REF] Meier | High-dimensional additive modeling[END_REF]: the penalized General Additive Model (penGAM). It applies to the same models as before but are especially designed for splines estimation. In the same spirit the individual functions are penalized, but since each function can be represented as the sum of linear combinations of basis functions. It turns out to be a groupLASSO problem. Their contribution is also to consider not only sparsity but also smoothness in the estimation. Because complex functions require many basis functions it is common in the splines settings to construct an over complete basis and then apply shrinkage on coefficients5 to have a smooth estimates, this is known as smoothing splines. This takes the form of a Ridge regression so it can be easily integrated inside the procedure. The main objective is to solve:

min f y -f (X ) 2 2 + J( f ) (18)
With the sparsity-smoothness penalty being:

J( f ) = λ 1 f j 2 + λ 2 ( f j (x)) 2 dx (19)
and because we can rewrite each f j (x) = ∑ K k=1 β j,k b j,k (x) as a sum of K basis b(•) then the problem can be written as:

min β y -Bβ 2 2 + λ 1 p ∑ j=1 β j B j B j β j + λ 2 β j Ω j β j (20)
Ω j composed of the inner products of the second derivatives of the basis functions.

Screening

In an equivalent manner on the screening side the Non-parametric Independence Screening procedure (NIS) has been introduced by [START_REF] Fan | Nonparametric independence screening in sparse ultra-high-dimensional additive models[END_REF] as a natural extension to SIS. Instead of marginal linear correlation they use the concept of "marginal utility", already defined in [START_REF] Fan | Ultrahigh dimensional feature selection: beyond the linear model[END_REF] for generalized linear models, and here set this marginal utility to be the sum of squared marginal residuals resulting from a non-parametric additive model.

ω j = n ∑ i=1 y i -f j (x i, j ) 2 (21)
The latter, with f j (x i j ) obtained by splines6 , gives a ranking of variables in the same way as SIS:

M = {1 ≤ j ≤ p : ω j > δ } (22)
Where δ is a predefined threshold. Usually this step does not ensures selection consistency so they rely on a external procedure, namely SpAM or penGAM. Because of the problem of weak signals Iterative Conditional SIS has been discussed exactly the same as Iterative Conditional SIS was for SIS, that is applying NIS on residuals, conditionally on primarily selected variables. It is worth mentionning the work of [START_REF] Zhang | Correlation rank screening for ultrahigh-dimensional survival data[END_REF] who developed Correlation Ranked SIS (CR-SIS). The main purpose is to allow for any monotonic transformation of y by using its cumulative distribution as the dependent variable.

ω j = Cov( f j (x j ), G(y)) 2 G(y) = 1 n n ∑ i=1 I(y i ≤ y) (23)
The resulting model is less restricted allowing a non-linear response.

Partial Linear Models

A Partial Linear model takes the form:

y = X 1 β + g(X 2 ) + ε (24)
An important feature of these models is to assume two sets of variables. The X matrix is divided into X 1 and X 2 of dimension p1 and p2 respectively. The motivation behind this is to say that linearity is satisfactory enough for some variables and treating these ones non-parametrically result in a loss of efficiency. So one should divide the regressors according to their link function either it is parametric (X 1 ) or not (X 2 ). This section is divided in two parts. The first one will concern Partial Linear models in their general form. Because a great literature has focused on smoothly varying-coefficients the second part will focus only on them.

Standard

Penalty

The Double-Penalized Least Squares Estimator (DPLSE) of [START_REF] Ni | Automatic model selection for partially linear models[END_REF] is a method for selection of variables and selection between parametric and non-parametric parts. A penalty is imposed on the parametric part to select variables and splines are used for non-parametric estimation. Since in the splines settings one can rewrite this function as a linear combination of basis expansion:

g = [J, X 2 ]δ + Ba (25)
with J the unit vector a are the parameters of the basis expansion B and δ is the overall parameter on X 2 . The SCAD penalty is then applied on the vector β * = [β , δ ] This can be viewed as a composite penalty where the key idea is to write everything as linear and perform usual model selection. Partial Splines with Adaptive penalty (PSA) of [START_REF] Cheng | Sparse and efficient estimation for partial spline models with increasing dimension[END_REF] try to achieve a sparse parametric part while having a non-parametric part aside using a combination of Adaptive LASSO on the parametric part and Penalized Splines for the non-parametric. Therefore the problem to solve is:

min β , f y -X 1 β -f (X 2 ) 2 + λ 1 1 0 ( f (X 2 )) 2 dX 2 + λ 2 p ∑ j=1 |β j | | β j | γ (26)
We remark the last term is exactly the penalty from the adaptive LASSO. This is in line with DPLSE, adding a smoothness penalty on top of the procedure. In this respect it is worth mentioning the Penalized Estimation with Polynomial Splines (PEPS) of [START_REF] Lian | Separation of covariates into nonparametric and parametric parts in high-dimensional partially linear additive models[END_REF]. The same objective is achieved in a quite similar fashion. The only difference is that the penalty is not adaptive:

min β , f y -Bβ 2 + nλ 1 p ∑ j=1 w 1, j β j A j + nλ 2 p ∑ j=1 w 2, j β j D j ( 27
)
Basis expansion is contained in B therefore exploiting once again the linear transformation provided in splines, just like DPLSE introduced it. The whole thing is turned as a linear model on which penalties are applied to achieve sparsity β j A j = ∑ k β j,k B k (x j ) and linear parts are recovered from the smoothness penalty

β j D j = ∑ k β j,k B k (x j ) .
In the end there is little difference between the 3 procedures. All exploits the linearity provided by splines.

PEPS improves on DPLSE adding a smoothness penalty and PSA improves on PEPS making the penalty adaptive to achieve better selection consistency.

Varying Coefficients

Another usual structure for modelling is the semi-varying coefficient model, written as:

y = X 1 β + X 2 α(Z) + ε (28)
The coefficients α associated to each x j∈2 are supposed to vary smoothly along another variable Z. This can be seen as a particular case of previous models where g(•) has the specific varying coefficient form.

Penalty

The methods in this section do not use the semi-structure form, they work only with the varying-coefficient part.

y = Xβ (Z) + ε (29)
The Kernel LASSO of [START_REF] Wang | Shrinkage estimation of the varying coefficient model[END_REF] deals with this problem in the spirit of groupLASSO. min

β n ∑ t=1 n ∑ i=1 {y i -X i β (Z t )} 2 K h (Zt -Z i ) + p ∑ j=1 λ j β j ( 30
)
The penalty enforces the procedure to reduce estimated varying coefficients close to zero to true zeros in a single-level group fashion. Another improvement in this setting is the Adaptive Semi-Varying Coefficients (AdaSVC) of [START_REF] Hu | Adaptive semi-varying coefficient model selection[END_REF]. Instead of all coefficients varying smoothly one may think that some don't (hence semi-varying). To avoid the loss of efficiency introduced by non-parametric estimation when the true underlying coefficient is constant the latter have to be identified. Their method can simultaneously identify and estimate such a model. Selection is done only over constant regressors. They do not consider sparsity as in Kernel LASSO. The idea is to impose a group penalty on the estimated varying-coefficients such that the penalty enforces nearly constant coefficients to be truly constant. Their penalty is in line with the FusedLASSO of [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF]. The main idea is that nearly constant coefficients will become constant in a grouped fashion. The objective is to solve:

min β n ∑ t=1 n ∑ i=1 {y i -X i β (Z t )} 2 K h (Zt -Z i ) + p ∑ j=1 λ j b j (31)
with the penalty applied on a different norm than the Kernel LASSO:

b j = n ∑ t=2 (β j (Z t ) -β j (Z t-1 )) 2 1/2
(32)

Testing

The Semi-Parametric Generalized Likelihood Ratio Test (SP-GLRT) of [START_REF] Li | Variable selection in semiparametric regression modeling[END_REF]. It applies to semi-varying coefficients model. The purpose is both to identify relevant variables and whether if they belong to the non-linear or the linear component. The likelihood can be written as:

L(α, β ) = l(α, β ) -n p ∑ j=1 p λ j (|β j |) (33)
The two parts are estimated alternatively conditionally on the other. Then they introduce a novel generalized likelihood ratio test:

T G LR = r K {R(H 1 ) -R(H 0 )} (34) with R(H 1 ) = Q(X 1 β + X 2 α(Z), y) (35)
The conditional likelihood under H 1 : at least one coefficient from the non-parametric part is non-zero.

R(H 0 ) = Q(X 1 β , y) (36)
The conditional likelihood under H 0 : the variable does not appear in the non-parametric part. where the conditional likelihood is given by:

Q(µ, y) = y µ s -y V (s) ds ( 37
)
The test is then evaluated using a Monte Carlo or Bootstrap methods to empirically estimates distribution of the statistics since the theoretical degrees of freedom tends to infinity preventing from a parametric test.

This has to be noticed because this is one of the first attempt of introducing non-parametric and therefore automatic tests inside a selection procedure. While methods like Autometrics and Stepwise Regression relies on parametric tests, SP-GLRT uses data-driven tests to construct the model. This idea of exploiting the data themselves to conduct tests is certainly not new, but it was in model selection. This idea is the core of methodologies for improving model selection in section 8.

Non-Parametric Models

A fully non-parametric model takes the form of:

y = f (X ) + ε (38)
Where f (•) is any multivariate function, linear or not, additive or not. This framework is very general therefore making it complicated for estimation. The most well known drawback is the Curse of Dimensionality. Briefly, it states that the number of observations required for estimation of this function grows exponentially with the dimension of X : p. It is already complicated to fit such a (maybe very non-linear) function non-parametrically in a reduced dimension, thus looking for a sparse representation is necessary when dealing with large p.

This time the different methods differ under several aspects. Testing ones like MARS shares similarities with Stepwise Regression for example, in an ANOVA Splines settings. Penalty ones uses ANOVA models also, the reason is that it limits interactions terms and gets closer to an additive model, this is indeed very common when dealing with fully non-parametric regression. The screening based ones can be divided in two categories: some make the use of generalized correlations to avoid using a model (DC-SIS,HSIC-SIS,KCCA-SIS,Gcorr) 7

while others rely on a specific model ex ante (MDI,MDA,RODEO). 8

Penalty

The Variable selection using Adaptive Nonlinear Interaction Structures in High dimensions (VANISH) of [START_REF] Radchenko | Variable selection using adaptive nonlinear interaction structures in high dimensions[END_REF]. It is very similar to the SHIM of [START_REF] Choi | Variable selection with the strong heredity constraint and its oracle property[END_REF] but in a non-linear framework. In order to approach the complexity of the function it uses an ANOVA-type of model defined as:

f (X ) = p ∑ j=1 f j (x j ) + ∑ j<k f j,k (x j , x k ) + ... + ε (39)
Where f j are the main effects, f j,k are the two-way interactions and so on. Their approach is closely related to the penGAM of [START_REF] Meier | High-dimensional additive modeling[END_REF] generalized to include interaction terms 9 but with a different penalty. The authors say that the penalty shouldn't be the same for main effect than for two-way interactions. They advocate the fact that ceteris paribus including an interaction term add more regressors than a main effect and thus that they are less interpretable. So interactions should be more penalized. Therefore this condition is a little bit different from the "strong heredity constraint" introduced in Choi et al. ( 2010). The objective is to solve:

min f y -f (X ) 2 2 + τ 2 J( f ) (40)
With:

J( f ) = λ 1 p ∑ j=1 f j 2 + ∑ k = j f j,k 2 1/2 + λ 2 p ∑ j=1 p ∑ k= j+1 f j,k (41)
The penalty is written so that the first part penalizes additional regressors while the second penalizes interactions occurring without main effects. In the SHIM there was no possibility for that. Here this constraint is released but a stronger penalty can be applied to restrict interactions without main effects, which are less interpretable.

Another approach for fitting this type of models is the Component Selection and Smoothing Operator (COSSO) of [START_REF] Lin | Component selection and smoothing in multivariate nonparametric regression[END_REF]. It differs from VANISH in the penalty function. The key idea is to use a penalty term written in terms of a sum of Reproducible Kernel Hilbert Space (RHKS) norms. In a model with only two-way interactions it would be:

J( f ) = p(p-1)/2 ∑ α=1 P α f 2 (42)
This time the penalty is not designed to take into account the structure of the resulting model. There is no desire to limit interactions. Since the heredity constraint is not present as before the model authors of VANISH claim it has trouble with high dimensional settings. Nevertheless the heredity constraint can obviously be inadequate in some applications where only interactions matter, in this type of settings COSSO is more advisable than VANISH.

7

Respectively Distance Correlation-SIS, Hilbert Schmidt Independence Criterion-SIS, Kernel Canonical Correlation Analysis and the Generalized Correlation.

8

Respectively Mean Decrease Impurity, Mean Decrease Accuracy and the Regularization Of Derivative Expectation Operator.

9

They also introduce it as SpIn (SpAM with INteractions) in their paper but claim that interactions would then not be treated efficiently.

Testing

Introduced by [START_REF] Friedman | Multivariate adaptive regression splines[END_REF] the Multivariate Adaptive Regression Splines is a method for building non-parametric fully non-linear ANOVA sparse models (39). The model is written in terms of splines as:

f (x) = K ∑ k=1 c k B k (x) (43)
The basis functions B k are taken to be hinge functions. The form of these functions makes the model piecewise linear.

B k (x, α, β ) = β max(0, α + x) (44)
Therefore α can be considered as "knots" like in standard splines. The β are parameters on which selection will occur through a pretty complicated algorithm. The building process is quite comparable to the one of usual Regression Trees and Stepwise Regression. Starting from a null model a forward step search over all possible variables and determines by least squares the parameter β (thus it creates a new hinge function) and over all possible values where to add a knot α that reduces best the residuals sum of squares10 . This process goes until some stopping criterion is met. All combinations have to be taken into account, therefore it is computationally intractable for high interactions effects. Friedman advises to limit the number of interactions m to a small value like 2 such that the model can be build in a reasonable time. Selection of variables is part of the building process.

If using a fit based criterion like the sum of squares residuals, variables are selected only if they bring enough explanatory power during the search. The same thing applies for Regression Trees on non-parametric models. In this sense MARS is closely related to Stepwise Regression. Also MARS is available with a backward approach, and a combination of both. This method is mainly used to fit high dimensional non-linear functions because since it is piecewise linear, it does not suffer much from the Curse of Dimensionality. However its selection consistency can be directly linked to the way variables are selected in trees, this is discussed in the next subsections. Used directly MARS is more like a non-linear version of Stepwise Regression using piecewise functions.

Screening

Model-free

In the screening literature of non-parametric methods we find a bunch of papers that deals with the same core idea. They all define some association measure that generalizes usual linear correlation. Here is the list of them as well as the criteria they use. In fact these methods are quite nested within each other. Considering which one is the best is a question of computational complexity rather than in which case they apply. Otherwise it seems that the last one (KCCA) should be selected.

• DC-SIS [START_REF] Li | Feature screening via distance correlation learning[END_REF] The Distance Correlation is a generalization of the Pearson Correlation Coefficient in terms of norm distances. It can be written as:

ω j = dcov(x, y) dcov(x, x)dcov(y, y) (45)
Where:

dcov(x, y) 2 = E[ X -X Y -Y ] +E[ X -X ]E[ Y -Y ] -2E[E[ X -X ]E[ Y -Y ]] (46)
• HSIC-SIS [START_REF] Balasubramanian | Ultrahigh dimensional feature screening via rkhs embeddings[END_REF] The Hilbert Schmidt Independence Criterion generalizes the previous one as it defines a maximum distance metric in a RKHS space:

ω 2 (k) = E[k X (X, X )k Y (Y ,Y )] +E[k X (X, X )]E[k Y (Y ,Y )] -2E[E[k X (X, X )]E[k Y (Y ,Y )]] (47)
We recognize again the form of the usual correlation but this time written in terms of kernels. In order to avoid the choice of the bandwidths in kernels, they decided to use the sup of the criterion over a family of Kernel K.

γ = sup ω (k) : k ∈ K (48)
Empirically the ranking measure is simpler to compute:

ω = 1 n sup k X ,k Y trace(K X HK Y H) (49)
with H = I -(1/n)JJ , I being the n × n unit matrix and J the n × 1 unit vector.

• KCCA-SIS [START_REF] Liu | Ultrahigh dimensional feature selection via kernel canonical correlation analysis[END_REF] The Kernel Canonical Correlation Analysis is the last improvement in the field of Non-parametric Screening.

It encompasses SIS as it can handle non-linearities. Unlike DC-SIS it is scale-free and does not rely on the Gaussian assumption. However even if it shares many aspects of the HSIC-SIS it differs in one aspect: HSIC is based on maximum covariance between the transformations of two variables, while KCCA uses the maximum correlation between the transformations by removing the marginal variations. Their measure is defined as:

R Y X = Σ -1/2 YY Σ Y X Σ -1/2 XX ( 50
)
Because the covariance matrices may not be invertible they introduce a ridge penalty ε:

R Y X = (Σ YY + εI) -1/2 Σ Y X (Σ XX + εI) -1/2 (51)
The correlation measure is then defined as the norm of the correlation operator:

ω(ε) j = R Y X j (52)
Empirical estimates of covariance matrices Σ are obtained after singular decomposition of kernel matrices (the latter being the same as in HSIC). While bandwidths in kernels can be chosen optimally ex ante, ε has to be estimated via GCV over a grid of values.

For each one the variables are ranked along marginal association measures ω j between y and x j and one defines the set of relevant features after applying a threshold. The latter's value differs among them.

M = {1 ≤ j ≤ p : ω j ≥ δ } (53) • DC-SIS: δ = cn -k • HSIC-SIS: δ = cn -k • KCCA-SIS: δ = cn -k ε -3/2 with 0 ≤ k ≤ 1/2.
Another of the same kind is the Generalized Correlation Screening (Gcorr) of [START_REF] Hall | Using generalized correlation to effect variable selection in very high dimensional problems[END_REF] that was introduced as a more general method than NIS. The general correlation coefficient is used as the measure of non-linear relationship. It can be defined as:

ω j = sup h∈H ∑ n i=1 {h(x i, j ) -h j }(y i -ȳ) n ∑ n i=1 {h(x i, j ) 2 -h2 j } (54)
Then these estimates are tested using bootstrap confidence interval instead of threshold like the others usually do.

Finally significant ones are ranked. Even though their method seems very general, empirically h(•) are chosen to be polynomial functions. This can be restrictive in some situations and less non-parametric in some sense.

Model Based

The Regularization Of Derivative Expectation Operator (RODEO) of [START_REF] Lafferty | Rodeo: sparse, greedy nonparametric regression[END_REF], named in reference to the LASSO, applies in the framework of Multivariate Kernel Methods. In kernel regression a specific attention is given to the choice of the bandwidth. We recall that this hyperparameter defines the width of the support for the regression, the lower it is the less observations enter the local regression, leading to less bias but more variance and conversely for a high bandwidth. The authors here state that for variables that are important in the model the derivative of the estimated function with respect to the bandwidth h is higher than for useless variables. A change in bandwidth affects the estimation if the variable intervenes in the model, it affects the bias-variance trade-off. For an irrelevant variable a change in bandwidth has no effect since more or less observations does not change the fitted curve. For a Gaussian kernel we have:

∂ f h (x) ∂ h j = e (X W X ) -1 X ∂W ∂ h j (y -X β ) ∂W ∂ h j = W L j L j = 1 h 3 j diag((x 1, j -x j ) 2 , ..., (x n, j -x j ) 2 ) (55)
Note that it refers to a specific point in the sample x. The derivative is not computed over the whole sample. The authors propose an extension of local RODEO to a global procedure where the derivative is computed in every point and then averaged.

The idea is to exploit this derivative iteratively, starting from a high bandwidth value and adapted in each step according to a certain rate of decay. Important variables should have low bandwidth, so the derivative is greater and the bandwidth reduces more quickly. Variables then can be ranked according to the final value of their bandwidth. One can apply some threshold on these to end up with a sparse solution. In this respect RODEO can be classified as a screening procedure. RODEO is based on a full estimation via Kernel, therefore it suffers the Curse of Dimensionality mentioned earlier. RODEO may not be able to deal with high dimensional feature space.

A large part of the literature focuses on a quite restricted set of regression methods for doing selection such as Ordinary Least Squares for linear models, Splines and Kernels for non-linear ones. However there exists other ways for doing regression from which model selection procedures intuitively arise. In a Bayesian framework11 one will consider a collection of models called an Ensemble. There is a distribution of them and we are uncertain on which one is the truth12 . Still we can exploit this distribution accross these different models to assign probabilities to each variables since they may not all appear in every models. This idea has also been developed in the frequentist approach by [START_REF] Breiman | Random forests[END_REF] who introduced Random Forest. From an Ensemble of Regression Trees (called a Forest) he derived two types of variables importance measures : Mean Decrease Impurity (MDI) and Mean Decrease Accuracy (MDA). We recall briefly that a tree is constructed as a recursive partitioning over the sample space. Simple Regression Trees allows for constant estimation in subregions, this is closely related to the Nadaraya-Watson local constant kernel estimator. Splits are chosen according to an impurity criterion that describes the degree of similarity13 of the data in the partition.

MDI(x j ) = 1 N t ∑ T ∑ t∈T N t N i(t) - N t le f t N t i(t le f t ) - N t right N t i(t right ) (56)
The importance of variable j is computed as the average decrease in impurity among each node t in tree T . The idea is to show the decrease in impurity caused by a split in this variable. It is computed as the impurity in the node minus the sum of impurity in the child nodes weighted by their respective sizes. This gain is weighted in the end by the number of observations entering the node. MDI can be easily extended to an Ensemble of Trees (i.e. a Forest).

The second measure relies on the predictive power of the model instead of the impurity inside nodes. From a statistical point of view it is equivalent as focusing on out-of-sample fit rather than in-sample fit. Since it does not rely on an inside criterion it is only defined for a tree and therefore applies only for an ensemble of them.

V I T (x j ) = ∑ i∈B (T ) I(y i = y (T ) i ) |B (T ) | - ∑ i∈B (t) I(y i = y (T ) i,π j ) |B (T ) | (57) MDA(x j ) = ∑ T ∈F V I T (x j ) N T (58)
The importance of variable j is computed as the average decrease in accuracy among each tree T in the forest F. The idea is that if a variable is uninformative then the prediction accuracy should be unchanged under permutation. The difference between actual prediction and permuted prediction give sthe decrease in accuracy for each variable and the whole is a weighted average of each tree in the forest.

Improving on Variable Selection

This last section is devoted to general methodologies designed for improving model selection procedures.

Based on bootstrap or resampling, the core idea is to exploit randomness to account for uncertainty in the modelling. Usual model selection procedures may suffer from inconsistency under some conditions. For example we remember the LASSO where the regularization parameter λ can not be chosen optimally 14 so that it ensures correct identification. This has lead to the adaptive LASSO [START_REF] Zou | The adaptive lasso and its oracle properties[END_REF], but this problem can also be solved using these procedures.

Stability Selection

The Stability Selection (Stabsel) has been introduced by [START_REF] Meinshausen | Stability selection[END_REF] to improve on selection. Given a specific selection procedure a variable is said to be stable if its selection probability under subsampling 15 (number of times it has been selected among the random samples) exceeds a specified threshold δ . The selection probabilities for a variable j to belong to the set S λ of selected variables for a given regularization parameter λ is:

Π λ j = P( j ⊆ S λ ) (59)
The set of stable variables is then:

S Stable = { j : max λ ∈Λ Π λ j ≥ δ } (60)
This is given by the underlying selection procedure, it can be the LASSO or whatever, but the methodology aims at improving a procedure, not being one itself. Another way for randomness that is almost equivalent is to divide the sample in two non-overlapping parts of sizes n/2 and look for variables that are selected simultaneously in both. This is more computationally efficient. The threshold can be selected appropriately so that the expected number of false inclusion V is bounded.

E[V ] ≤ 1 2δ -1 q 2 λ p (61)
Thus one will ensure P(V > 0) ≤ α by setting for example: δ = 0.9 q λ = 0.8α p (62)

The results are then presented as stability paths: Π λ j as a function of λ . This is in contrast to regularization paths of LASSO: β j as a function of λ .

Extensions to Stabsel are proposed in [START_REF] Bach | Bolasso: model consistent lasso estimation through the bootstrap[END_REF] and [START_REF] Shah | Variable selection with error control: another look at stability selection[END_REF]. The first uses bootstrap with replacement instead of resampling without while the latter uses subsampling of complementary pairs.

Ranking-Based Variable Selection

The Ranking-Based Variable Selection (RBVS) of [START_REF] Baranowski | Ranking-based variable selection for high-dimensional data[END_REF] is a screening procedure based on bootstrap and permutation tests. Contrary to Stabsel it does not rely on any threshold nor any assumptions. Given a metric to assess the strength of the relationship denoted ω and then using the m-out-of-n bootstrap of Bickel et al. ( 2012) they construct a permutation ranking R. R = (R 1 , ..., R p ) satisfying ω R 1 ≥ ... ≥ ω R p (63)

The metric can be anything like the Pearson Correlation, the LASSO coefficients, etc. The probability of the set of the k top-ranked variables A k is defined as:

π(A k ) = P({R 1 , ..., R k } = A) (64)
This value is approximated with using the m-out-of-n bootstrap procedure involving random draws without replacements of the observations. In fact selection can be performed on the set of top-ranked variables A from which the number of terms k * can be determined automatically without threshold. The idea is not to look for a threshold δ that would cut in the ranking of ω. As an alternative they try to estimate k * as:

k * = argmin k=0,...,p-1 π(A k+1,m ) π(A k,m ) (65)
That is the number of terms for which the differences among the π(A) is the greater. This is equivalent to look for a threshold that best separates assuming there are two sets: the relevant and the irrelevant. It has the advantage of being totally non-parametric. Just like the SIS has its iterative counterpart they introduce the Iterative RBVS that accounts for marginally related variables with low Signal-to-Noise and for the multicollinearity problem.

Discussion

In this article, we provide a review for 39 state-of-the-art procedures to perform variable selection over a wide range of model structures going from the simple linear one to the complex non-parametric one. Procedures have been classified in three groups: Tests-Based, Penalty-Based and Screening-Based. They have been described and compared on the ground of model structures.

The main difference consists of modelling purposes and objectives rather than their strength as oracles. In an empirical work the choice between two strategies should rely on the form of the model, data specificities (collinearity,groups, etc..) and objectives (in other words understandability).

Selection consistency for widely used methods in empirical work have been discussed and several improvements were presented. Far beyond Stepwise Regression and the LASSO, empiricists have access to more advanced technologies that we claim are not much more complicated than the basic ones. The limits in main methods (LASSO, Stepwise Regression) are now well understood and various answers have come to light.

The area of model selection is still very investigated, much more now that amounts of data have become available. Nevertheless, methods for handling large number of variables are restricted in terms of model complexity. This is mainly due to the Curse of Dimensionality and it prevents from looking for very complex models in high dimensions. Sure Independence Screening is a powerful tool in linear models but have lower dataset capacities when it comes to non-linearities. Also, the literature is lacking from very complete algorithmic solutions. To the best of our knowledge, no statistical procedure have been developed to reach the level of completeness of Autometrics. Other methods are only parts of the statistical work and do not cover as many problems as Autometrics do.

•

  Tests-based • Penalty-based • Screening-based Originally, the first developed are based on statistical tests. The work was to automate standard tests in Econometrics (like testing residuals for normality, t-tests etc.) for choosing among candidate variables. It includes Stepwise Regression and Autometrics for example.

Table 1 .

 1 Topology of Variable Selection Methods

		Screening	Penalty	Testing
		SIS	SparseStep	Stepwise
		SFR	LASSO	Autometrics
		CASE	Ridge
	Linear	FA-CAR	BRidge SCAD
			MCP
			NNG
			SHIM
			gLASSO
	Group		gBridge gSCAD
			gMCP
			ElasticNet
		NIS	SpAM
	Additive	CR-SIS	penGAM
			kernelLASSO	SP-GLRT
			adaSVC
	Partial Linear		DPLSE
			PSA
			PEPS
		DC-SIS	VANISH	MARS
		HSIC-SIS	COSSO
	Non-Parametric	KCCA-SIS
		Gcorr
		MDI
		MDA
		RODEO

Usually the Ridge because it has an analytical solution.

Because of low computational costs, but it can be estimated with any non-parametric regression technology.

This is known as "Greedy Algorithms" where the optimal global solution is sought by taking optimal local solutions.

Which is out of the scope of this paper but still very important.

This relates obviously to the problem raised when discussing Stepwise Regression. Here the Ensemble is a subset of the model space.

In case of a regression: How well the subregion can be approximated by a constant.

Without replacement, random samples have to be non-overlapping.


Tom Eeckelaert
email: tom.eeckelaert@esat.kuleuven.ac.be

Trent Mcconaghy

Georges Gielen

Efficient Multiobjective Synthesis of Analog Circuits using Hierarchical Pareto-optimal Performance Hypersurfaces

An efficient methodology is presented to generate the Pareto-optimal hypersurface of the performance space of a complete mixed-signal electronic system. This Paretooptimal front offers the designer access to all optimal design solutions: starting from the performance specifications, a satisfactory point can a posteriori be selected on the hypersurface which immediately determines the final design parameters. Fast execution is guaranteed by using multi-objective evolutionary optimization techniques and hierachical decomposition. The presented method takes advantage of the Pareto hypersurfaces of the subblocks to generate the overall Pareto front. The hierarchical approach combines behavioral simulation with behavioral models at the higher levels, with SPICE simulations with transistor-level accuracy at the lowest level. Storing the performance data of all subblocks enables reuse for other systems later on.

Introduction

To keep up with the drastically increasing complexity of todays electronic systems, a lot of research has been done in the field of analog design automation [START_REF] Gielen | Computer-aided design of analog and mixed-signal integrated circuits[END_REF]. A lot of prototype CAD tools resulted from the research [START_REF] Rutenbar | Computer-Aided Design of Analog Integrated Circuits and System[END_REF]: IDAC , OASYS, OPTIMAN, AMGIE, ASTRX/OBLX, ANACONDA.

The most successful tools use optimization techniques to search for the optimal design solution according to the given specifications. Most of these tools however are limited to rather small building blocks. Because of the increasing complexity of electronic systems, there is a shift in the analog/mixed-signal research towards design methodologies for higher hierarchical design levels. By hierarchically decomposing the system in smaller less complex building blocks, dedicated simulators can be used at each hierarchical level. Examples are SD-OPT [START_REF] Medeiro | A vertically-integrated tool for automated design of Σ∆ modulators[END_REF] and DAISY [START_REF] Francken | A high-level simulation and synthesis environment for ∆Σ modulators[END_REF] for the high-level design of Sigma-Delta converters. However, these methodologies simulate at a higher hierarchical level without taking into account lower-level simulations. When enclosing simulations of lower-level building blocks, the design space explodes and these methodologies will fail to find an optimum solution. Moreover, because the majority of the existing tools is based on the optimization of one single objective function, the optimum that is found is depending on the a-priori choice of the weighting coefficients used to combine the different objectives. Therefore, the methodology presented in this paper is based on the hierarchical use of powerful Multiobjective Evolutionary Algorithms (MOEA) [START_REF] Van Veldhuizen | Multiobjective Evolutionary Algorithms: classifications, analysis, and new innovations[END_REF]:

• Through their inherent parallelism, these EA are well suited in handling the hierarchical decomposition of large systems.

• These multiobjective techniques are able to find the Pareto-optimal trade-off hypersurface between the different objective functions, in one single optimization run.

• To cope with the disadvantage of the exploding design space when including lower-level simulations of the building blocks, we will show that lower-level Pareto-fronts can be used to store the hierarchical performance information and make the whole approach very efficient.

• This also makes it possible to simulate each hierarchical building block by the best suitable simulator. This way, the use of SPICE-like simulators at the transistor level preserves the accuracy that these simulators provide.

• The generated Pareto-optimal data points at each level can be stored and modeled into a continuous Paretohypersurface. The stored data of the subblocks can be reused in future other systems without the need of resimulation.

The paper is organized as follows. Section 2 describes some definitions about multi-objective optimization and the algorithm that was selected for use in the presented methodology.

In section 3, the methodology is explained and choices in the implementation are discussed. Next, in section 4, the method is illustrated by generating the Pareto-optimal performancehypersurface of a ∆Σ converter is generated with and without the use of lower-level Pareto-optimal surfaces, and the two methods are compared. The final section, gives conclusions about the research.

Multi-objective Optimization

In this section first some basic definitions about multiobjective optimization (MOO) and Pareto optimality are given. Next, the Strength Pareto Evolutionary Algorithm (SPEA), which was used to implement the presented methodology, is briefly discussed. discussed.

Basic Definitions

In a general multi-objective optimization problem the goal is to find optimal values for the components of a vector function f (x) subject to some constraints e(x):

minimize/maximize y = f (x) = ( f 1 (x), f 2 (x),... , f n (x)) subject to e(x) = e 1 (x), e 2 (x),...,e k (x) ≤ 0, with x = (x 1 , x 2 ,... ,x m ) ∈ X, y = (y 1 , y 2 ,... ,y n ) ∈ Y (1)
where x is the decision vector, y is the objective vector, X is called the decision space and Y is called the objective space. The k constraints e(x) ≤ 0, set restrictions on the decision space. The multiple objective functions almost always conflict with one another. Therefore, for multi-objective problems, Pareto optimality is used. First the concept of Pareto dominance should be clear. For a maximization problem, a vector u = (u 1 ,... ,u m ) dominates a vector v = (v 1 ,... ,v m ) (u v) if and only if ∀i ∈ {1, 2,...,n} :

f i (u) ≥ f i (v)

and

(2) ∃ j ∈ {1, 2,...,n} :

f j (u) > f j (v).
Additionally, u covers v (u v) if and only if u v or u = v. A decision vector x ∈ X is now said to be Pareto-optimal if and only if there exists no z ∈ X for which f (z) dominates f (x).

Strength Pareto Evolutionary Algorithms

Evolutionary algorithms are stochastic optimization methods designed by analogy with natural evolution [START_REF] Zitzler | Evolutionary Algorithms for Multiobjective Optimization: Methods and Applications[END_REF]. Starting from a set of solutions (individuals) subsequent sets are composed by means of selection and variation.

We integrated the existing Strength Pareto Evolutionary Algorithm (SPEA) [START_REF] Zitzler | Multiobjective Evolutionary Algorithms: A Comparative Case Study and the Strength Pareto Approach[END_REF] in the implementation of our methodology. Characteristics of the SPEA are:

• The Pareto-optimal solutions found in each iterative optimization run, are stored in an external set. After each run, the external set is updated by comparison with the newly found Pareto-optimal solutions. This way, Pareto-optimal solutions don't get lost.

• The important fitness assignment step during selection uses the concept of Pareto dominance [START_REF] Gielen | Computer-aided design of analog and mixed-signal integrated circuits[END_REF]. Only the externally stored Pareto set (size N p ) determines the fitness of an individual. Domination by individuals in the population (size N) is irrelevant. In this step, each of the individuals in the population have to be compared with the individuals in the external Pareto set, for all objectives (M). This is the most time consuming step, resulting in a computing complexity of O(MNN p ).

• To reduce the size of the Pareto-optimal set, clustering is used in the original implementation. In particular, the average linkage method. At the first optimization step, each individual represents one cluster on its own. Then, in each step two clusters of individuals are chosen, using the nearest neighbor criterion, to merge into one cluster. After that, out of each cluster one individual is chosen to represent the cluster in the external set.

In our opinion it is better not to loose any Pareto-optimal solutions, so we didn't use the clustering technique. A disadvantage to this is the increasing individuals in the external set (O(MNN p )), but there is a gain from the eliminated clustering step.

System design space Exploration

In this section, first hierarchical decomposition and behavioral performance modeling, and the relevance to our methodology is discussed. In the second subsection it is discussed how the hierarchical performance information of the lowerlevel building blocks can be included in the generation of the Pareto-optimal hypersurface at a higher level.

Hierarchical decomposition and behavioral performance modeling

To cope with the increasing complexity of current analog systems, the systems have to be decomposed into smaller less complex subsystems. When the subsystems are still to complex to design, a second decomposition is performed. This systematic, hierarchical decomposition can go on until all subblocks are manageable for design. If at some hierarchical level a subblock has already been designed in the past, then the decomposition of the branch which comprises that block will end at the hierarchical level of that block. Figure 1 illustrates this. When at some branch in the decomposition, the system has never before been designed, then the lowest hierarchical level is the transistor level. At transistor level, the block can be simulated by a SPICE-like simulator to extract its performance behavior. If on the other hand a subsystem has already been used for another system, the past simulation results can be reused. The subblock is then said to be used as intellectual property (IP). After decomposition of the system, at each higher level, behavioral models have to be used to combine the performance information of lower levels. A lot of research has been published about behavioral performance modeling, and the majority of them is already implemented in some kind of highlevel simulator, e.g. VHDL-AMS. Our methodology provides the capability of using all possible simulators for a building block. The system is defined to the tool by means of simple description files. For each subsystem a description file has to be composed containing information about the simulator to be used, the design variables, the subsystems from which it uses performance information, and which performance variables have to be extracted from the output file of the simulator.

Modeling and reusing lower-level Pareto surfaces

As was mentioned in the previous subsection, performance trade-off information from previous simulation runs can be stored and reused when a subblock is used again in another system. A look-up table could be used, but it would even be best to find a mathematical model for the hypersurface from which information can be extracted without performing a simulation. In this way, only the very first time a subblock is used, transistor-level simulations are needed to create the data and generate the Pareto hypersurface, which can then by used for all later hierarchical synthesis runs in which that subblock is used. Advantages of this approach are:

• Algorithm speedup by saving time in the absence of lower-level simulation runs: either a lookup table used instead, or the mathematical model for the Pareto-optimal performance trade-off hypersurface has to be evaluated.

• When no lower-level performance information is used, the design variables of the high-level system, as well as the design variables from all the lower-level building blocks have to be determined during the optimization. This way the design space gets very large resulting in a huge population set (large N) and a huge externally stored set (large N p ) needed to sample that space. Through the use of lowerlevel Pareto-optimal performance hypersurfaces, the computing complexity (O(MNN p )) of the algorithm can thus be reduced.

The following derivation however shows that attention should be paid to the way the lower-level Pareto-optimal surfaces are generated: In a hierarchical design, performance characteristics from a lower level are used as design variables at the higher level. Suppose for the sake of simplicity a minimization problem with a two-dimensional decision space (x 1 , x 2 ) and a twodimensional objective space ( f 1 (x 1 , x 2 ), f 2 (x 1 , x 2 )). If it can be shown that there exists a dominated point (definitions in section 2.1) in the decision space which contributes to the Pareto-optimal curve in the objective space, then it is proven that essential information is lost when only the nondominated points (Pareto-optimal 2D-curve) in the decision space are used. Suppose a point (x * 1 , x * 2 ) in the decision space is dominated

⇒ ∃(x 1 , x 2 ) | x 1 ≤ x * 1 & x 2 ≤ x * 2 (3)
For this point not to contribute to a Pareto-optimal objective point,

( f * 1 , f * 2 ) = ( f 1 (x * 1 , x * 2 ), f 2 (x * 1 , x * 2 )) should be dominated. So, ∃( f 1 , f 2 ) | f 1 ≤ f * 1 & f 2 ≤ f * 2 (4)
Combining ( 3) and ( 4) results in the 2 following requirements: If ( 5) and ( 6) are not met, then the Pareto optimal performance trade-off curve (2D-hypersurface) contains points corresponding with dominated points in the performance space of a lower hierarchical building block. So, in order to use lower-level Pareto-optimal surfaces, a lower-level solution moving towards its Pareto-optimal surface, should result in a movement of its corresponding higherlevel solution towards the higher-level Pareto-optimal surface. As shown in the above, for both a minimization problem at the two hierarchical levels, the derivatives of the objective functions considered should be positive. For the performance parameters 'power' and 'area' for example, which are most typically used objectives in analog synthesis, this requirement is obviously met. As an other example, for the calculation of the Signal-to-Noise Ratio (SNR), calculated with the highlevel simulator depicted in [START_REF] Martens | A Model of Computation for Continuous-Time ∆Σ Modulators[END_REF], only the position of poles, zeros and the gain of the composing integrators is used from the lower level. Because the relationship between the gain of the different integrators and the SNR isn't as the requirements above state, gain should not be included as an objective at the integrator transistor level; it can be handled as a constraint.

1) δ f 1 δx 1 x * 1 > 0 & δ f 1 δx 2 x * 2 > 0 or if δ f 1 δx 1 x * 1 < 0 ⇒ δ f 1 δx 1 x * 1 < δ f 1 δx 2 x * 2 or if δ f 1 δx 2 x * 2 < 0 ⇒ δ f 1 δx 2 x * 2 < δ f 1 δx 1 x * 1 (5) 2) δ f 2 δx 1 x * 1 > 0 & δ f 2 δx 2 x * 2 > 0 or if δ f 2 δx 1 x * 1 < 0 ⇒ δ f 2 δx 1 x * 1 < δ f 2 δx 2 x * 2 or if δ f 2 δx 2 x * 2 < 0 ⇒ δ f 2 δx 2 x * 2 < δ f 2 δx 1 x * 1 ( 6
If on the other hand certain lower-level performance parameters do not influence higher-level objectives, then they can be included in the list of lower-level objectives as well. In the example of section 4 it will be shown that for performance parameters at higher and lower hierarchical levels that meet the requirements stated in the above derivation, significant benefit can be derived from using the lower-level Pareto-optimal surface instead of lower-level simulations.

The optimization framework

The hierarchical multi-objective synthesis method has been prototyped in a software framework. The framework contains three major fragments (Figure 2). In the initialization part, the complete hierarchically decomposed system is set up. As described in section 3.1, a description file has to be given for the complete system. This description file is then analyzed in the System Setup module and description files for the subsystems are systematically included during setup. The second major part, the simulation part, consists of a module (FileWriter) that uses the information extracted from the description files to generate a new input file for each externally used simulator. The generation of such an input file is based on a template file which contains variables for which values have to be substituted. These values are extracted from lowerlevel building blocks or given as design parameters during simulation startup. After writing the input files, the different external simulators are executed. When they are finished, a Simulation Extractor module extracts the needed performance and constraint variable values from the simulator outputs. All simulations are stored in files for later use (see sections 3.1 and 3.2). If on the other hand a flag is set to use previously stored simulations or a set of Pareto-optimal solutions, then no file is written and no external simulators are called. The Pareto-optimal solutions are stored in an ordered set to make it possible to mutate to nearby optimal solutions. The third fragment, implements the SPEA optimizer as discussed in section 2.2. In each optimization run, the design variables which can not be extracted from lower levels are generated by the MOEA for each subsystem in the hierarchy. Values for the design variables are chosen within their boundaries given in the description files. The complete systems performance variables are then determined by a call to the simulator. The externally stored Pareto set is updated and by assigning fitness values, individuals are selected to be used for mating (section 2.2) and generating the next set of design variables.

Experimental results

In this section we will illustrate the presented methodology by generating the Pareto-optimal performance-hypersurface of a third-order single-loop continuous-time delta-sigma modulator (Figure 3). We chose a simple hierarchical decomposition of one subblock for the ∆Σ converter: the integrator subblock is represented by a highspeed folded-cascode OTA (Figure 4), implemented in a 0.25 µm CMOS technology. For the high-level simulation of the ∆Σ, the Model of Computation from [START_REF] Martens | A Model of Computation for Continuous-Time ∆Σ Modulators[END_REF] was implemented in the MATLAB TM programming language. At transistor-level, HSpice was used. The information contained in the description file of the ∆Σ modulator is given below:

VSS VSS M2b M3b M1a M1b n2a n3 n2b nin_p nin_m n5 n4a n4b M2a M6 M7 M4a M5a M5b M4b n1 nout M3a i B1 i B2 V SS V DD i B3 i D1 i D2
• The performances chosen for the exploration are power consumption [W], silicon area [m 2 ], and signal-to-noise ratio (SNR) [dB].

• The design variables at the modulator level are the oversampling ratio (OSR), the input-signal frequency ( f 0 ), the sampling frequency ( f s ), and the input amplitude (A 0 ). Also the poles, zeros and gain of the integrators are needed but they get extracted from the lower level during the simulation.

• No constraint variables were set for the modulator.

For the description file of the OTA, the following information is given: • The design variables include all independent gate-source and drain-source voltages, transistor widths (operatingpoint driven formulation [START_REF] Leyn | An Efficient DC Root Solving Algorithm with Guaranteed Convergence for Analog Intergrated CMOS Circuits[END_REF]), and biasing currents, and the integrating capacitance. The supply voltages are fixed to Vdd = 1.25V and Vss = -1.25V.

•
• Constraints are drain-source and gate-source voltages to keep the transistors into saturation. Also the poles and zeros are given in the list of constraints, but without constraints, just to make sure they are determined for use at ∆Σ level.

To demonstrate the advantage of using hierarchical performance information in the form of Pareto-optimal hypersurfaces, two methods are compared here. One (A), where the design space of the ∆Σ modulator combines all design spaces of all the integrators and its own, and another method (B) where first a Pareto-optimal solution set is generated for the integrators and then this set is used during optimization of the In accordance to the definitions from [START_REF] Zitzler | Performance Assessment of Multiobjective Optimizers: An Analysis and Review[END_REF], the following relations between objective vectors (z i ) are used: First an optimization was started to find an approximation of the Pareto-optimal set for all 3 integrators. Keep in mind that this should normally only be done once, because the resulting set can be used for all three integrators and for future design as well. This is only done here to have different simulation data when we use method A. Table 1 shows that the fronts don't dominate each other. Method A consists of highlevel optimization while using all generated simulation results stored while generating the lower-level Pareto-optimal surfaces. Method B uses only final the Pareto-optimal sets of the integrators. The total amount of optimization generations at ∆Σ level was limited to 100. While a population limit of 100 individuals was set, the external Pareto set was not limited. For method A, a total of about 6300 good individuals were generated and an other 12200 individuals were discarded because they didn't meet the specs (e.g. transistors in saturation). For method B, a total of about 5400 good individuals were generated and about 12000 were discarded. Pareto-optimal performance information doesn't abundantly deteriorate the quality of the Pareto-optimal set (solution set of method B shouldn't be strictly dominated by the solution set of method A). As can be concluded from figure 5, method B even results in solution sets that dominate the solution sets of method A. The randomly generated initial population generated in method B is already slightly dominating the initial population in method A. This can be explained by the fact that the 'power' and 'area' performance parameters meet the requirements depicted in subsection 3.2, and because in method B only Pareto-optimal lower-level values are used the higher-level values of 'power' and 'area' are already Paretooptimal. Over the 100 generations, the solution set from method B almost always dominates the solution set of method A. For later generation however C (A, B) and C (A, B) curves are slightly rising again and C (B, A) and C (B, A) curves are decreasing. So, the solution sets of method A start to cath up with the solution sets of method B. Eventually both solutions should be at least as good because all the Pareto-optimal lower-level solutions used in method B are also contained in the lower-level solution set used by A. This only takes more generation runs.

z A ∈ A strictly dom- inates z B ∈ B (z A z B ) if

Conclusions

In this paper, an efficient methodology was presented for the generation of Pareto-optimal performance tradeoff surfaces of complex systems. The objective was to use hierarchical performance information from lower hierarchical subblocks in the form of Pareto-optimal performance hypersurfaces. It was shown that under certain requirements, these lower-level Pareto-optimal fronts can indeed be used to extract performance information of the subblocks avoiding time-consuming transistor-level simulations. Advantages of this approach are the limited dimension of the high-level design space, the storage memory reduction by modeling the Pareto-optimal performance trade-off surface, the significant speedup of the optimization algorithm, while the quality of the eventually generated high-level performance trade-off is maintained. The surfaces generated by this methodology can be used by an analog designer to a priori choose a performance configuration that best meets the given specifications and hereby immediately determines the values of all design parameters of all the subsystems.

Figure 1 .

 1 Figure 1. Decomposition of a complex system. Ending nodes in the branches represent subsystems for which no lower-level performance information is needed.

Figure 2 .

 2 Figure 2. Algorithm framework.

Figure 3 .

 3 Figure 3. A single-loop 3 th order Delta-Sigma modulator architecture.

Figure 4 .

 4 Figure 4. CMOS folded-cascode operational transconductance amplifier.

  The performances are power [W], area [m 2 ], UGF [Hz], phase margin (PM) [ • ], noise density [V 2 /Hz], slew rate [V/s]. They can all be extracted from the HSpice output file.

Int 1 2 3 C 1 C 2 C 3 C

 123123 Int Int = 1.00 C = 0.12 C = 0.16 Int = 0.00 C = 0.10 C = 0.14 C = 0.33 C = 1.00 C = 0.26 Int = 0.31 C = 0.00 C = 0.25 C = 0.20 C = 0.13 C = 1.00 Int = 0.18 C = 0.11 C = 0.00 ∆Σ level. To compare the two solution sets A and B, we selected the Coverage metric C (A, B) from [9]. This metric represents the amount of individuals in population B that are dominated by an individual in population A divided by the total amount of individuals in B. So, C (A, B) = 1 means that all individuals in B are dominated by or equal to individuals in A. C (A, B) = 0 means that none of the individuals in B are dominated.

Figure 5 Figure 5 .

 55 Figure 5. Comparison of ∆Σ Pareto-optimal fronts generated by method A and B, at different generations in the optimization process starting from the initial population, using the coverage metric C .

Table 1 . Comparison of the integrators Pareto-optimal fronts using C metrics.

 1

  all objective values in z A are better than in z B , z A dominates z B (z A z B ) if all objective values in z A are better or equal than in z B . Using these relations, a distinction is made between C (A, B) and C (A, B). For example, if 0 < C (A, B) < 1 and 0 < C (B, A) < 1 then both sets contain solutions which dominate solutions in the other set. But if for the same sets C (A, B) = 0 and C (B, A) > 0 it can be concluded that only B contains solutions which strictly dominate solutions in A, thus B is better then A in relation to approximating the Pareto-optimal front.

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05) 1530-1591/05 $ 20.00 IEEE

Proceedings of the Design, Automation and Test in Europe Conference and Exhibition (DATE'05)