{'fraction_non_alphanumeric': 0.04869575816674793, 'fraction_numerical': 0.02029497805948318, 'mean_word_length': 4.028807845540913, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Best subset selection in linear regression via bi-objective mixed integer linear programming

We study the problem of choosing the best subset of p features in linear regression given n observations. This problem naturally contains two objective functions including minimizing the amount of bias and minimizing the number of predictors. The existing approaches transform the problem into a single-objective optimization problem. We explain the main weaknesses of existing approaches, and to overcome their drawbacks, we propose a bi-objective mixed integer linear programming approach. A computational study shows the efficacy of the proposed approach.


Introduction
Linear regression models should have two important characteristics in practice including prediction accuracy and interpretability (Tibshirani, 1996). The traditional approach of constructing regression models is to minimize the sum of squared residuals. It is evident that models obtained in this approach have low biases. However, their prediction accuracy can be low due to their large variances. Furthermore, models constructed by this approach may contain a large number of predictors and so data analysts struggle in interpreting them.
In general, reducing the number of predictors in a regression model can improve not only the interpretability but also, sometimes, the prediction accuracy by reducing the variance (Tibshirani, 1996). Hence, there is often a trade-off between the amount of bias and the practical characteristics of a regression model. In other words, finding a desirable regression model is naturally a biobjective optimization problem that minimizes the amount of bias and the number of predictors, simultaneously. This classical Statistics problem, that is finding the best subset of p potential predictors and estimating their coefficients in a linear regression model given n observations is called the Best Subset Selection Problem (BSSP). Due to the intensive computation of solving bi-objective optimization problems, BSSP has always been transformed into a single-objective optimization problem (e.g., see Ren and Zhang (2010) and Bertsimas et al. (2016)). Such transformation may cause practical issues that will be discussed in this paper.
Albeit, recent algorithmic and theoretical advances have made bi-objective optimization problems computationally tractable in practice. More precisely, under some mild conditions, we are now able to solve them reasonably fast. To the best of our knowledge, this paper is the first work to tackle BSSP utilizing a bi-objective optimization approach.
The structure of the paper is organized as follows: In Section 2, two main single-objective approaches for BSSP are presented and their practical drawbacks are discussed. In Section 3, a new bi-objective optimization formulation for BSSP is developed. In Section 4, the computational results are reported. Finally, in Section 5, some concluding remarks are provided.
2 Single-objective models for BSSP As discussed in Section 1, BSSP is naturally a bi-objective optimization problem which can be stated as min β∈F {z 1 ( β), z 2 ( β)} where F is the feasible set of parameter estimator vector β, z 1 ( β) is the total bias and z 2 ( β) is the number of predictors. In the literature, the following two approaches have widely been used to convert BSSP to a single-objective optimization problem: (i) The weighted sum approach: Given some λ > 0, BSSP has been reformulated as min β∈F z 1 ( β) + λz 2 ( β).
(ii) The goal programming approach: Given some k ∈ Z ≥ , BSSP has been reformulated as min β∈F : z 2 ( β)≤k For further details, interested readers are referred to Chen et al. (1998), Meinshausen and Bhlmann (2006), Zhang and Huang (2008), Bickel et al. (2009), Candés andPlan (2009), Ren and Zhang (2010), Ciuperca (2014) and Wang and Tian (2017) for the weighted sum approach, and Miller (2002) and Bertsimas et al. (2016) for the goal programming approach. Although, those two optimization programs (i) and (ii) can be solved significantly faster than a bi-objective optimization problem, their drawbacks are explained and illustrated here. However, we first need to present some notation and definitions from bi-objective optimization.

Preliminaries
A Bi-Objective Mixed Integer Linear Program (BOMILP) can be stated as follows: where X := {(x 1 , x 2 ) ∈ Z n 1 ≥ × R n 2 ≥ : A 1 x 1 + A 2 x 2 ≤ b} represents the feasible set in the decision space, Z n 1 ≥ := {s ∈ Z n 1 : s ≥ 0}, R n 2 ≥ := {s ∈ R n 2 : s ≥ 0}, A 1 ∈ R m×n 1 , A 2 ∈ R m×n 2 , and b ∈ R m . It is assumed that X is bounded and z i (x 1 , x 2 ) = c ⊺ i,1 x 1 + c ⊺ i,2 x 2 where c i,1 ∈ R n 1 and c i,2 ∈ R n 2 for i = 1, 2 represents a linear objective function. The image Y of X under vectorvalued function z := (z 1 , z 2 ) ⊺ represents the feasible set in the objective/criterion space, that is x 2 ) for all (x 1 , x 2 ) ∈ X }. Note that BOMILP is called Bi-objective Linear Program (BOLP) and Bi-Objective Integer Linear Program (BOILP) for the special cases of n 1 = 0 and n 2 = 0, respectively.
Definition 1. A feasible solution (x 1 , x 2 ) ∈ X is called efficient or Pareto optimal, if there is no is called a nondominated point. The set of all efficient solutions is denoted by X E . The set of all nondominated points z(x 1 , x 2 ) for (x 1 , x 2 ) ∈ X E is denoted by Y N and referred to as the nondominated frontier.
Definition 3. Let Y e be the set of extreme points of the convex hull of Y, that is the smallest convex The convex hull of Y  In summary, based on Definition 1, the elements of Y can be divided into two groups including dominated and nondominated points. Furthermore, based on Definitions 2 and 3, the nondominated points can be divided into unsupported nondominated points, non-extreme supported nondominated points and extreme supported nondominated points. Overall, bi-objective optimization problems are concerned with finding all elements of Y N , that is all nondominated points, including supported and unsupported nondominated points. An illustration of the set Y and its corresponding categories are shown in Figure 1. It is well-known that in a BOLP, both the set of efficient solutions X E and the set of nondominated points Y N are supported and connected (Sayin, 1996). Consequently, to describe all nondominated points in a BOLP, it suffices to find all extreme supported nondominated points. A typical illustration of the nondominated frontier of a BOLP is displayed in Figure 2a.
Since we assume that X is bounded, the set of nondominated points of a BOILP is finite.
However, due to the existence of unsupported nondominated points in a BOILP, finding all nondominated points is more challenging than in a BOLP. A typical nondominated frontier of a BOILP is shown in Figure 2b where the rectangles are unsupported nondominated points.
Finding all nondominated points of a BOMILP is even more challenging. Nonetheless, if at most one of the objective functions of a BOMILP contains continuous decision variables, then the set of nondominated points is finite and BOILP solution approaches can be utilized to solve it (Stidsen et al., 2014). However, in all other cases that more than one objective function contain continuous decision variables, the nondominated frontier of a BOMILP may contain connected parts as well as supported and unsupported nondominated points. Therefore, in these cases, the set of nondominated points may not be finite and BOILP algorithms cannot be applied to solve them anymore. A typical nondominated frontier of a BOMILP is illustrated in Figure 2c for further discussions on the properties of BOILPs and BOMILPs and algorithms to solve them.
The wighted sum approach fails to find rectangles The goal programming approach may incorrectly find rectangles

Drawbacks of transforming BSSP into a single-objective model
Suppose that for each β ∈ F, the corresponding point (z 1 ( β), z 2 ( β)) ⊺ is plotted into the criterion space. Figures 3a and 3b show two typical plots of such pairs for all β ∈ F. In these two figures, all filled circles and rectangles are nondominated points of the problem and unfilled rectangles and circles are dominated points. In Figure 3a, the region defined by the dashed lines is the convex hull of all feasible points. In this case, it is impossible that the wighted sum approach finds the filled rectangles for any arbitrary weight as all filled rectangles are unsupported nondominated points (i.e., they are interior points of the convex hull). So, this illustrates that there may exist many nondominated points, but the weighted sum approach can fail to find most of them for any arbitrary weight. Figure 3b is helpful for understanding the main drawback of the goal programming approach. It is obvious that depending on the value of k, the goal programming approach may find one of the unfilled rectangles which are dominated points. So, the main drawback of the goal programming approach is that it may even fail to find a nondominated point.
The main contribution of our research presented here is to overcome both of these disadvantages by utilizing bi-objective optimization techniques. We note that in the literature of BSSP, z 1 ( β) is mainly defined as the sum of squared residuals. The reason lies in the fact that the sum of squared residuals is a smooth (convex) function. However, to be able to exploit existing bi-objective mixed integer linear programming solvers, we use the sum of absolute residuals for z 1 ( β). We conclude this section by providing two remarks.
Remark 4. If we incorporate additional linear constraints on the vector of parameter estimators of the regression model, β, it will be more likely that the goal programming approach fails to find a nondominated point.
Remark 5. Unlike the weighted sum and goal programming approaches that new parameters λ and k, respectively, should be employed and tuned by the user, the bi-objective optimization approach does not need any extra parameter.

A bi-objective formulation for BSSP
Let X = [x 1 , . . . , x p ] ∈ R n×p be the model matrix (it is assumed that x 1 = 1), β ∈ R p×1 be the vector of regression coefficients, and y ∈ R n×1 be the response vector. It is assumed that β is unknown and should be estimated. Let β ∈ R p×1 denote an estimate for β. To solve BSSP for this set of data, we construct the following BOMILP and denote it by BSSP-BOMILP: such that: r j l j ≤ β j ≤ r j u j for j = 1, . . . , p p j=1 x ij β j − y i ≤ γ i for i = 1, . . . , n r j ∈ {0, 1} for j = 1, . . . , p β j ∈ R for j = 1, . . . , p, where l j ∈ R and u j ∈ R are, respectively, a lower bound and an upper bound (known) for β j for j = 1, . . . , p, γ i for i = 1, . . . , n is a non-negative continuous variable that captures the value of |y i − p j=1 x ij β j | in any efficient solution, and r j for j = 1, . . . , p is a binary decision variable that takes the value of one if β j = 0, implying that the predictor j is active. By these definitions, for any efficient solution, the first objective function, z 1 ( β) = n i=1 γ i , takes the value of the sum of absolute residuals and the second objective function, z 2 ( β) = p j=1 r j , computes the number of predictors. Constraint (3) ensures that if β j = 0 then r j = 1 for j = 1, . . . , p. Constraints (4) and (5) guarantee that |y i − p j=1 x ij β j |≤ γ i for i = 1, . . . , n. Note that since we minimize the first objective function, we have |y i − p j=1 x ij β j |= γ i for i = 1, . . . , n in an efficient solution.
Remark 6. The BSSP-BOMILP can handle additional linear constraints and variables. Furthermore, by choosing tight bounds in Constraint (3), we can speed up the solution time of BSSP-BOMILP.
Hence, we should try to choose l j /u j as large/small as possible.
Remark 7. Since only one of the objective functions in BSSP-BOMILP contains continuous variables, based on our discussion in Section 2.1, the set of nondominated points of BSSP-BOMILP is finite.
is a trivial nondominated point of BSSP-BOMILP where there is no parameter selected in the estimated regression model. Hence, we exclude this trivial nondominated point by adding the constraint p j=1 r j ≥ 1 to BSSP-BOMILP.

Bounds for the regression coefficients
In this section, we develop a data-driven approach to find bounds l j and u j for j = 1, . . . , p such that l j ≤ β j ≤ u j , in the lack of any additional information. To achieve this, we firstly present the following proposition.
Motivated from Proposition 9, we can solve the following optimization problem to find u j for j = 1, . . . , p: There are several ways to transform (9) into a linear program (e.g., see Dielman (2005)). Here, we propose the following linear programing model: It should be noted that (10) is a relaxation of (9) since γ i over-calculates |y i − p j ′ =1 x ij ′ β j ′ | for i = 1, . . . , n. Analogously, l j for j = 1, . . . , p can be computed by changing 'max' into 'min' in (10).

Computational results
We conduct a computational study to show the performance of the ǫ-constraint method on BSSP-BOMILP, numerically. We use C++ to code the ǫ-constraint method. In this computational study, the algorithm uses CPLEX 12.7 as the single-objective integer programming solver. All computational experiments are carried out on a Dell PowerEdge R630 with two Intel Xeon E5-2650 2.2 GHz 12-Core Processors (30MB), 128GB RAM, and the RedHat Enterprise Linux 6.8 operating system.
We allow CPLEX to employ at most 10 threads at the same time.
We design six classes of instances, each denoted by C(p, n) where p ∈ {20, 40} and n ∈ {2p, 3p, 4p}. Based on this construction, we generate three regression models for each class as follows: • Set all x i1 = 1 and all other x ij for j > 1 are randomly drawn from the discrete uniform distribution on interval [−50, 50]; • Two steps are taken to construct y i for i = 1, . . . , n: (1) A vector β is generated such that the two third of its components are zero, and the others are randomly drawn from the uniform distribution on interval (0, 1); (2) Set y i = p j=1 x ij β j + ε i (with at most one decimal place) where ε i is randomly generated from the standard normal distribution; • Optimal values of l j and u j for j = 1, . . . , p are computed by solving (10). and p = 40 in about 5 seconds and 7 minutes in average, respectively.
To highlight the drawbacks of existing approaches including the weighted sum approach and the goal programming approach, the nondominated frontier of the Instance 1 from the Class C(20, 40) is illustrated in Figure 4. The filled rectangles and circles are unsupported and supported nondominated points, respectively. As we discussed previously, it is impossible to find any of the unsupported nondominated points using the weighted sum approach. Also, observe that many of the nondominated points lies on an almost vertical line. This implies that all these points are almost optimal for the goal programming approach when k = 7, . . . , 20.
We note that selecting a desirable nondominated point in the nondominated frontier depends on decision maker(s). Here, we introduce a heuristic algorithm to do so. Let T = (T 1 , T 2 ) ⊺ ∈ R 2 and B = (B 1 , B 2 ) ⊺ ∈ R 2 be the top and bottom endpoints of the nondominated frontier. One may simply choose the point that has the minimum Euclidean distance from the (imaginary) ideal point, that is (T 1 , B 2 ) ⊺ . Based on this algorithm, in Figure 4, the (imaginary) ideal point is (22.6, 0) ⊺ and the closest nondominated point to it is (27.2, 8) ⊺ . The generated instance that we discuss in Figure 4 is y = 0.42 + 0.86x 2 + 0.63x 3 + 0.68x 4 + 0.42x 5 + 0.50x 6 + 0.25x 7 and the estimated linear regression model corresponding to the selected nondominated point is y = 0.4 + 0.87x 2 + 0.65x 3 + 0.68x 4 + 0.43x 5 + 0.51x 6 + 0.24x 7 + 0.01x 14 , which are very close together.

Conclusion
In this paper, we discussed the practical drawbacks of transforming the best subset selection problem in linear regression modeling into a single-objective optimization problem. To resolve those disadvantages, we developed a new bi-objective optimization model for BSSP. The efficacy of this new model was shown through numerical results. Hopefully, the simplicity, versatility and per-formance of this new bi-objective optimization model encourage practitioners to consider it for constructing linear regression models.