<jupyter_start><jupyter_text>Best Subset Selection: L0-Regression Objective and PrerequisitesIn this example, you will learn how to perform linear regression with feature selection using mathematical programming. We’ll show you how to construct a mixed-integer quadratic programming (MIQP) model of this linear regression problem, implement this model in the Gurobi Python API, and generate an optimal solution.This modeling example is at the intermediate level, where we assume that you know Python and are familiar with the Gurobi Python API. In addition, you should have some knowledge about building mathematical optimization models.**Download the Repository** You can download the repository containing this and other examples by clicking [here](https://github.com/Gurobi/modeling-examples/archive/master.zip). --- MotivationLinear regression was invented at the beginning of the 19th century and today, after more than 200 years, it is still used extensively in practical applications for description and prediction purposes:- In econometrics, it is useful to estimate the price elasticity of a particular product by regressing sales revenue on price and possibly other features such as demographics and competitor and retail information.- In health sciences, it can be applied to predict how long a patient will remain (i.e. length of stay) in the ER of a hospital based on patient information, triage assessment, medical test results, and the date/time of arrival.- In social sciences, it may shed light on future academic performance of students, so proactive measures can be taken to improve their learning outcomes.In general, linear regression is used to model the relationship between a continuous variable and other explanatory variables, which can be either continuous or categorical. When applying this technique, finding the subset of features that maximizes its perfomance is often of interest.--- Problem DescriptionLinear regression is a supervised learning algorithm used to predict a quantitative response. It assumes that there is a linear relationship between the feature vector $x_i \in \mathbb{R}^d$ and the response $y_i \in \mathbb{R}$. Mathematically speaking, for sample $i$ we have $y_i = \beta^T x_i + \epsilon_i$, where $\beta \in \mathbb{R}^d$ is the vector of feature weights, including the intercept, and  $\epsilon_i$ is a normally-distributed random variable with zero mean and constant variance representing the error term. We can learn the weights from a training dataset with $n$ observations $\{X \in \mathbb{M}^{nxd},y \in \mathbb{R}^n\}$ by minimizing the Residual Sum of Squares (RSS): $e^Te =(y-X\beta)^T (y-X\beta)=\beta^T X^T X\beta- 2y^TX\beta+y^T y$. The Ordinary Least Squares (OLS) method achieves this by taking the derivative of this quadratic and convex function and then finding the stationary point: $\beta_{OLS}=(X^T X)^{-1} X^T y$.In practice, some of the features are in fact not associated with the response. By including them, we only add unnecessary complexity to the model and increase variance to the weight estimates. However, finding the best performing model is no simple task as there is an exponential number of candidates, as one has to test $\sum_{s=1}^{d-1}{{d-1} \choose s}$ models. Since OLS rarely yields estimates that are exactly zero, thus discarding the features related to them, we need to resort to feature selection methods. Popular methods include:- Subset selection, e.g. stepwise selection.- Dimensionality reduction, e.g. principal component regression.- Shrinkage, e.g. the Lasso.The Lasso has undoubtedly been the method of choice for the last decade. Basically it fits a model containing all $d$ predictors, while incorporating a budget constraint based on the L1-norm of $\beta$, disregarding the intercept component. In fact, this method minimizes the RSS, subject to $\sum_{l=1}^{d-1}\mathopen|\beta_l\mathclose| \leq s$, where $s$ is a hyper-parameter representing the budget that is usually tuned via cross-validation. This constraint has the effect of shrinking all weight estimates, allowing some of them to be exactly zero when $s$ is small enough. Finally, it is worth noting that the unconstrained version of the Lasso is more frequently used. This version solves an unconstrained optimization problem where $RSS + \lambda \sum_{l=1}^{d-1}\mathopen|\beta_l\mathclose|$ is minimized, for a given value of the —modified— lagrangian multiplier $\lambda \in \mathbb{R}^+$.  A similar formulation is now presented, where the L0-norm is used instead. We now seek to minimize the RSS, subject to $\sum_{l=1}^{d-1}I(\beta_l \neq 0) \leq s$, where $I(\beta_l \neq 0)$ is an indicator function taking on the value of 1 if $\beta_j \neq 0$ and 0 otherwise. In this setting, $s$ represents the number of features to consider in the model. This optimization problem may be casted as a Mixed Integer Quadratic Program (MIQP). Traditionally, the feature selection problem has not been tackled this way because of the common belief in the statistics community that large-scale problems are intractable. But this is no longer the case, considering the computing power currently available and the performance of modern optimization solvers such as Gurobi.--- Solution ApproachMathematical programming is a declarative approach where the modeler formulates a mathematical optimization model that captures the key aspects of a complex decision problem. The Gurobi Optimizer solves such models using state-of-the-art mathematics and computer science. A mathematical optimization model has five components, namely:- Sets and indices.- Parameters.- Decision variables.- Objective function(s).- Constraints.We now present a MIQP formulation that finds the weight estimates for a linear regression problem, where exactly $s$ of those weights can be nonzero: Sets and Indices$i \in I=\{1,2,\dots,n\}$: Set of observations.$j \in J=\{1,2,\dots,d\}$: Set of features, where the last ID corresponds to the intercept.$l \in L = J \backslash \{d\}$: Set of features, where the intercept is excluded. Parameters$s \in \mathbb{N}$: Number of features to include in the model, ignoring the intercept.$Q = X^T X \in \mathbb{M}^{|J|x|J|}$: Quadratic component of the objective function.$c = y^TX \in \mathbb{R}^{|J|}$: Linear component of the objective function.$k = y^T y \in \mathbb{R}$: Constant component of the objective function.***Note:** Recall that the RSS is defined as $\beta^T X^T X\beta- 2y^TX\beta+y^T y$. Decision Variables$\beta_j \in \mathbb{R}$: Weight of feature $j$, representing the change in the response variable per unit-change of feature $j$.$iszero_l \in \{0,1\}$: 1 if weight of feature $l$ is exactly equal to zero, and 0 otherwise. Auxiliary variable used to manage the budget constraint. Objective Function- **Training error**: Minimize the Residual Sum of Squares (RSS):\begin{equation}\text{Min} \quad Z = \frac{1}{2}\sum_{j \in J}\sum_{j' \in J}Q_{j,j'}\beta_j\beta_{j'} - \sum_{j \in J}c_j\beta_j + \frac{1}{2}k\tag{0}\end{equation}Note that we use the fact that if $x^*$ is a minimizer of $f(x)$, it is also a minimizer of $a\cdot f(x)$, as long as $a > 0$. Constraints- **Mutual exclusivity**: For each feature $l$, if $iszero_l=1$, then $\beta_l=0$ (defined as a Special Ordered Set of type one):\begin{equation}(\beta_l, iszero_l): \text{SOS-1} \quad \forall l \in L\tag{1}\end{equation}- **Budget constraint**: Exactly $|L| - s$ feature weights must be equal to zero:\begin{equation}\sum_{l \in L}iszero_l = |L| - s\tag{2}\end{equation}This model, by means of constraint 2, implicitly considers all ${{d-1} \choose s}$ feature subsets at once. However, we also need to find the value for $s$ that maximizes the performance of the regression on unseen observations. Notice that the training RSS decreases monotonically as more features are considered (which translates to relaxing the MIQP), so it is not advisable to use it as the performance metric. Instead, we will estimate the Mean Squared Error (MSE) via cross-validation. This metric is defined as $\text{MSE}=\frac{1}{n}\sum_{i=1}^{n}{(y_i-\hat{y}_i)^2}$, where $y_i$ and $\hat{y}_i$ are the observed and predicted values for the ith observation, respectively. Then, we will fine-tune $s$ using grid search, provided that the set of possible values is quite small.--- Python ImplementationIn the following implementation, we use four main libraries:- **Numpy** for scientific computing.- **Scikit learn** for machine learning algorithms.- **Gurobi** for mathematical optimization.- **Matplotlib** for visualization.<jupyter_code>%pip install gurobipy
# Let's import all the necessary libraries

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from itertools import product
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import cross_val_score
from sklearn import linear_model
from sklearn.metrics import mean_squared_error as mse

import gurobipy as gp
from gurobipy import GRB

from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
# Create and deploy the optimization model

# NOTE: This function assumes the design matrix features does not contain
#a column for the intercept
def miqp(features, response, non_zero, warm_up=None, verbose=False):
    """
    Deploy and optimize the MIQP formulation of L0-Regression.
    """
    assert isinstance(non_zero, (int, np.integer))
    regressor = gp.Model()
    samples, dim = features.shape
    assert samples == response.shape[0]
    assert non_zero <= dim

    # Append a column of ones to the feature matrix to account for the y-intercept
    X = np.concatenate([features, np.ones((samples, 1))], axis=1)  
    
    # Decision variables
    beta = regressor.addMVar(dim + 1, lb=-GRB.INFINITY, name="beta") # Weights
    intercept = beta[dim] # Last decision variable captures the y-intercept
    intercept.varname = 'intercept'
    normvar = regressor.addVar(ub=non_zero, name='normvar')
    
    # Objective Function (OF): minimize 1/2 * RSS using the fact that
    # if x* is a minimizer of f(x), it is also a minimizer of k*f(x) iff k > 0
    Quad = np.dot(X.T, X)
    lin = np.dot(response.T, X)
    obj = 0.5 * (beta @ Quad @ beta)
    obj -= lin @ beta
    obj += 0.5 * np.dot(response, response)
    regressor.setObjective(obj, GRB.MINIMIZE)
    
    # Constraint sets
    regressor.addGenConstrNorm(normvar, beta[:dim], 0.0, "normconstr")
    
    if not verbose:
        regressor.params.OutputFlag = 0
    regressor.params.timelimit = 60
    regressor.params.mipgap = 0.001
    regressor.optimize()

    coeff = np.array([beta[i].X for i in range(dim)])
    return intercept.X, coeff     
# Define functions necessary to perform hyper-parameter tuning via cross-validation

def split_folds(features, response, train_mask):
    """
    Assign folds to either train or test partitions based on train_mask.
    """
    xtrain = features[train_mask,:]
    xtest = features[~train_mask,:]
    ytrain = response[train_mask]
    ytest = response[~train_mask]
    return xtrain, xtest, ytrain, ytest

def cross_validate(features, response, non_zero, folds, standardize, seed):
    """
    Train an L0-Regression for each fold and report the cross-validated MSE.
    """
    if seed is not None:
        np.random.seed(seed)
    samples, dim = features.shape
    assert samples == response.shape[0]
    fold_size = int(np.ceil(samples / folds))
    # Randomly assign each sample to a fold
    shuffled = np.random.choice(samples, samples, replace=False)
    mse_cv = 0
    # Exclude folds from training, one at a time, 
    #to get out-of-sample estimates of the MSE
    for fold in range(folds):
        idx = shuffled[fold * fold_size : min((fold + 1) * fold_size, samples)]
        train_mask = np.ones(samples, dtype=bool)
        train_mask[idx] = False
        xtrain, xtest, ytrain, ytest = split_folds(features, response, train_mask)
        if standardize:
            scaler = StandardScaler()
            scaler.fit(xtrain)
            xtrain = scaler.transform(xtrain)
            xtest = scaler.transform(xtest)
        intercept, beta = miqp(xtrain, ytrain, non_zero)
        ypred = np.dot(xtest, beta) + intercept
        mse_cv += mse(ytest, ypred) / folds
    # Report the average out-of-sample MSE
    return mse_cv

def L0_regression(features, response, folds=5, standardize=False, seed=None):
    """
    Select the best L0-Regression model by performing grid search on the budget.
    """
    dim = features.shape[1]
    best_mse = np.inf
    best = 0
    # Grid search to find best number of features to consider
    for i in range(1, dim + 1):
        val = cross_validate(features, response, i, folds=folds,
                             standardize=standardize, seed=seed)
        if val < best_mse:
            best_mse = val
            best = i
    if standardize:
        scaler = StandardScaler()
        scaler.fit(features)
        features = scaler.transform(features)
    intercept, beta = miqp(features, response, best)
    return intercept, beta<jupyter_output><empty_output><jupyter_text>--- BenchmarkWe now compare the performance of the aforementioned approach w.r.t. OLS regression on all features and the Lasso. The Boston dataset is used for this purpose. This dataset measures the prices of 506 houses, along with 13 features that provide insights about their neighbourhoods. We will use the original feature terminology, so the interested reader can visit [this website](https://www.cs.toronto.edu/~delve/data/boston/bostonDetail.html) for more information.Note that 20% of the samples are reserved for computing the out-of-sample MSE. The resulting metrics are displayed in a bar chart (shown below) to facilitate the comparison between models.<jupyter_code># Define how the bar chart should be displayed

def plot_bar_chart(performance):
    """
    Display the performance of all three models in a bar chart.
    """
    bar = plt.bar([1, 2, 3], performance, color=['r', 'g', 'y'],
                  tick_label=['OLS', 'Lasso', 'L0-Regression'])
    plt.title('Out-of-Sample MSE')
    x1, x2, y1, y2 = plt.axis()
    plt.axis((x1, x2, np.floor(np.min(performance)),
              np.ceil(np.max(performance))))
    plt.show()
# Load data and split into train (80%) and test (20%)
boston = load_boston()
X = boston['data']
y = boston['target']
Xtrain, Xtest, ytrain, ytest = train_test_split(X, y, test_size=0.20,
                                                random_state=10101)

# OLS regression using all features
lr = linear_model.LinearRegression()
lr.fit(Xtrain, ytrain)
# Lasso with cross-validated penalization (lambda)
lasso = linear_model.LassoCV(cv=5)
lasso.fit(Xtrain, ytrain)
# L0-regression where the best feature subset is selected via cross-validation
intercept, beta = L0_regression(Xtrain, ytrain, seed=10101)

# Compare their performance using a bar chart
performance = []
performance.append(mse(ytest, lr.predict(Xtest)))
performance.append(mse(ytest, lasso.predict(Xtest)))
performance.append(mse(ytest, np.dot(Xtest, beta) + intercept))
plot_bar_chart(performance)<jupyter_output><empty_output><jupyter_text>Notice that the Lasso performs poorly, as we did not standardize the features to be expressed in the same units (with an average of zero and variance of one). Contrary to OLS and L0-Regression, the Lasso is not scale-invariant because the budget constraint is based on the L1-norm. Remember that $\beta_l$ is interpreted as the change in the response per unit-change of feature $l$. Since the L1-norm takes the sum of absolute values, how much of the budget $\beta_l$ consumes depends on the units of measurement of the feature associated to it. Such preprocessing entails three steps, namely:For each feature $x_l$:1. Compute its sample average $\mu_l$ and sample standard deviation $\sigma_l$.2. Center by subtracting $\mu_l$ from $x_l$.3. Scale by dividing the resulting difference by $\sigma_l$.In order to report the performance of the Lasso after applying standardization, we need to perform hyper-parameter tuning on the L1-norm penalty via cross-validation. Unfortunately, we must not use the model class `LassoCV`. This is due to the fact that standardization is not supported, and doing that beforehand over the whole dataset would contaminate the folds. In order to prevent that from happening, we will perform random search as follows:<jupyter_code>np.random.seed(10101)
num_tries = 500
best_alpha = None
best_score = -np.inf
for i in range(num_tries):
    # log-linear search for alpha in the domain [0.001, 1000]
    exponent = np.random.uniform(-3, 3)
    alpha = np.power(10, exponent)
    pipeline = make_pipeline(StandardScaler(), linear_model.Lasso(alpha=alpha))
    scores = cross_val_score(pipeline, Xtrain, ytrain, cv=5, scoring='neg_mean_squared_error')
    avg_score = np.mean(scores)
    if avg_score > best_score:
        best_score = avg_score
        best_alpha = alpha<jupyter_output><empty_output><jupyter_text>Let's now compare the performance of the models considered when the features are preprocessed. Notice that our user-defined function `L0-regression` does support standardization of the features:<jupyter_code># Standardize the features so they have an avg of 0 and a sample var of 1
scaler = StandardScaler()
scaler.fit(Xtrain)
Xtrain_std = scaler.transform(Xtrain)
Xtest_std = scaler.transform(Xtest)

# OLS regression using all features
lr = linear_model.LinearRegression()
lr.fit(Xtrain_std, ytrain)

# Lasso with cross-validated penalization (lambda)
lasso = linear_model.Lasso(alpha=best_alpha)
lasso.fit(Xtrain_std, ytrain)
# L0-regression where the best feature subset is selected via cross-validation
intercept, beta = L0_regression(Xtrain, ytrain, standardize=True, seed=10101)

# Compare their performance using a Bar chart
performance = []
performance.append(mse(ytest, lr.predict(Xtest_std)))
performance.append(mse(ytest, lasso.predict(Xtest_std)))
performance.append(mse(ytest, np.dot(Xtest_std, beta) + intercept))
plot_bar_chart(performance)<jupyter_output><empty_output><jupyter_text>As expected, the Lasso is better (although marginally) than OLS. This is due to the fact that the Lasso can retrieve the $\beta_{OLS}$ estimate when the budget $s$ is big enough (alternatively, when $\lambda$ is small enough). On the other hand, it is marginally worse than L0-Regression, mainly because by shrinking $\beta$ we add bias to the estimates. Furthermore, observe that L0-Regression achieved the best performance with the fewest number of features. This is convenient, as it leads to a more interpretable model.<jupyter_code>ols_features = np.sum(np.abs(lr.coef_) >= 1e-8)
lasso_features = np.sum(np.abs(lasso.coef_) >= 1e-8)
l0_features = np.sum(np.abs(beta) >= 1e-8)
print("OLS regression kept {0} features.".format(ols_features))
print("The Lasso kept {0} features.".format(lasso_features))
print("L0-Regression kept {0} features.".format(l0_features))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08975031472933277, 'fraction_numerical': 0.008917331095258078, 'mean_word_length': 4.31502648452746, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 27, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '148060', 'license': 'Apache-2.0', 'n_tokens_mistral': 5540, 'n_tokens_neox': 5271, 'n_words': 2288, 'path': 'linear_regression/l0_regression_gcl.ipynb', 'repo_name': 'pobonomo/modeling-examples', 'size': '19064'}