Note

This is example is available as a Jupyter notebook. Download it and all necessary data files here.

L1 Norm Regression

This example implements a regression model which minimizes mean absolute error as a fitting criteria. This metric cannot be handled by the typical Ordinary Least Squares (OLS) regression implementation, but is well suited to linear programming (and therefore, to Gurobi!).

L1 norm regressionaims to choose weights \(w\) in order to minimize the following loss function:

\[\min_w \lvert Xw - y \rvert\]

To model the L1 regression loss function using linear programming, we need to introduce a number of auxiliary variables. Here \(I\) is the set of data points and \(J\) the set of fields. Response values \(y_i\) are predicted from predictor values \(x_{ij}\) by fitting coefficients \(w_j\). To handle the absolute value, non-negative continuous variables \(u_i\) and \(v_i\) are introduced.

\[\begin{split}\begin{aligned} \min \quad & \sum_i u_i + v_i \\ \mbox{s.t.} \quad & \sum_j w_j x_{ij} + u_i - v_i = y_i \quad & \forall i \in I \\ & u_i, v_i \ge 0 \quad & \forall i \in I \\ & w_j \,\, \text{free} \quad & \forall j \in J \end{aligned}\end{split}\]

We start with the usual imports for gurobipy-pandas models. Additionally, we import some sklearn modules, both to load an example dataset and to help perform model validation.

[1]:

import gurobipy as gp
from gurobipy import GRB
import gurobipy_pandas as gppd
import pandas as pd

from sklearn import datasets
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split

We first load the diabetes regression dataset into pandas dataframes, then perform scaling and a train-test split. Additionally, we add a constant column to capture an intercept term.

[3]:

diabetes_X, diabetes_y = datasets.load_diabetes(return_X_y=True, as_frame=True)
scaler = StandardScaler()
diabetes_X_scaled = pd.DataFrame(
    data=scaler.fit_transform(diabetes_X),
    columns=diabetes_X.columns,
    index=diabetes_X.index,
)
diabetes_X_scaled["const"] = 1.0
X_train, X_test, y_train, y_test = train_test_split(
    diabetes_X_scaled, diabetes_y, random_state=42
)
X_train

[3]:
age sex bmi bp s1 s2 s3 s4 s5 s6 const
16 -0.115937 -0.938537 0.889214 1.038895 0.516642 -0.501640 1.564414 -0.830301 1.099061 0.586922 1.0
408 1.335088 -0.938537 -1.059520 2.269385 0.661281 0.406893 -0.370637 0.496320 1.220180 0.848171 1.0
432 0.189542 -0.938537 1.161130 -0.119214 1.210908 0.940162 -0.061029 0.488562 1.170736 2.241496 1.0
316 0.342282 1.065488 0.300062 0.025550 0.024870 -0.448971 -0.680245 0.721302 1.576064 0.848171 1.0
3 -1.872441 -0.938537 -0.243771 -0.770650 0.256292 0.525397 -0.757647 0.721302 0.476983 -0.196823 1.0
... ... ... ... ... ... ... ... ... ... ... ...
106 -2.025181 -0.938537 -1.603353 -0.915414 -0.958674 -0.732065 0.171178 -0.830301 -1.250310 -1.764314 1.0
270 0.113172 1.065488 0.639957 1.762713 -0.785107 -0.995407 0.325982 -0.830301 0.181658 0.325674 1.0
348 0.647761 -0.938537 -0.425049 -0.119214 -0.090841 -0.620144 1.641816 -0.830301 -0.229228 -0.022657 1.0
435 -0.268676 -0.938537 -0.493028 -0.843032 -0.351191 0.097465 -0.370637 -0.054499 -0.808569 -0.806403 1.0
102 -1.948811 -0.938537 0.594638 -0.336359 0.776992 0.525397 1.177404 -0.830301 -0.108108 -0.022657 1.0

331 rows × 11 columns

Model Formulation

Create the Gurobi model:

[4]:

model = gp.Model()

Then we create an unbounded variable for each column coefficient. Here we use the free function gppd.add_vars to create these variables. Note that what we want here is one variable per column, and in pandas, the column labels are also an index. So, we can pass this index directly to the gurobipy-pandas function to create a Series of variables aligned to this column index.

[5]:

coeffs = gppd.add_vars(model, X_train.columns, name="coeff", lb=-GRB.INFINITY)
model.update()
coeffs

[5]:

age        <gurobi.Var coeff[age]>
sex        <gurobi.Var coeff[sex]>
bmi        <gurobi.Var coeff[bmi]>
bp          <gurobi.Var coeff[bp]>
s1          <gurobi.Var coeff[s1]>
s2          <gurobi.Var coeff[s2]>
s3          <gurobi.Var coeff[s3]>
s4          <gurobi.Var coeff[s4]>
s5          <gurobi.Var coeff[s5]>
s6          <gurobi.Var coeff[s6]>
const    <gurobi.Var coeff[const]>
Name: coeff, dtype: object

Given this column-oriented series and our training dataset, we can formulate the linear relationships representing the regression formula using native pandas operations:

[6]:

relation = (X_train * coeffs).sum(axis="columns")
relation

[6]:

16     -0.11593687913484971 coeff[age] + -0.938536660...
408    1.3350883235096633 coeff[age] + -0.93853666088...
432    0.18954211089557407 coeff[age] + -0.9385366608...
316    0.34228160591078594 coeff[age] + 1.06548847975...
3      -1.8724410718097866 coeff[age] + -0.9385366608...
                             ...
106    -2.0251805668249987 coeff[age] + -0.9385366608...
270    0.11317236338796813 coeff[age] + 1.06548847975...
348    0.6477605959412097 coeff[age] + -0.93853666088...
435    -0.2686763741500616 coeff[age] + -0.9385366608...
102    -1.9488108193173925 coeff[age] + -0.9385366608...
Length: 331, dtype: object

To incorporate the absolute error component, we need to introduce the variables \(u_i\) and \(v_i\), and constrain them such that they measure the error for each data point. This can be done compactly using the gppd dataframe accessor and method chaining:

[7]:

fit = (
    relation.to_frame(name="MX")
    .gppd.add_vars(model, name="u")
    .gppd.add_vars(model, name="v")
    .join(y_train)
    .gppd.add_constrs(model, "target == MX + u - v", name="fit")
)
model.update()
fit

[7]:
MX u v target fit
16 -0.11593687913484971 coeff[age] + -0.938536660... <gurobi.Var u[16]> <gurobi.Var v[16]> 166.0 <gurobi.Constr fit[16]>
408 1.3350883235096633 coeff[age] + -0.93853666088... <gurobi.Var u[408]> <gurobi.Var v[408]> 189.0 <gurobi.Constr fit[408]>
432 0.18954211089557407 coeff[age] + -0.9385366608... <gurobi.Var u[432]> <gurobi.Var v[432]> 173.0 <gurobi.Constr fit[432]>
316 0.34228160591078594 coeff[age] + 1.06548847975... <gurobi.Var u[316]> <gurobi.Var v[316]> 220.0 <gurobi.Constr fit[316]>
3 -1.8724410718097866 coeff[age] + -0.9385366608... <gurobi.Var u[3]> <gurobi.Var v[3]> 206.0 <gurobi.Constr fit[3]>
... ... ... ... ... ...
106 -2.0251805668249987 coeff[age] + -0.9385366608... <gurobi.Var u[106]> <gurobi.Var v[106]> 134.0 <gurobi.Constr fit[106]>
270 0.11317236338796813 coeff[age] + 1.06548847975... <gurobi.Var u[270]> <gurobi.Var v[270]> 202.0 <gurobi.Constr fit[270]>
348 0.6477605959412097 coeff[age] + -0.93853666088... <gurobi.Var u[348]> <gurobi.Var v[348]> 148.0 <gurobi.Constr fit[348]>
435 -0.2686763741500616 coeff[age] + -0.9385366608... <gurobi.Var u[435]> <gurobi.Var v[435]> 64.0 <gurobi.Constr fit[435]>
102 -1.9488108193173925 coeff[age] + -0.9385366608... <gurobi.Var u[102]> <gurobi.Var v[102]> 302.0 <gurobi.Constr fit[102]>

331 rows × 5 columns

Each step of the above code appends a new column to the returned dataframe:

  1. Start with the column “MX”, representing the \(Xw\) term

  2. Append a column of new variables \(u_i\) (one for each data point)

  3. Append a column of new variables \(v_i\) (one for each data point)

  4. Join the true target data from y_train

  5. Add one constraint per row such that u_i - v_i measures the error between prediction and target

Now we are in a position to formulate the mean absolute error expression, and minimize it in the objective function:

[8]:

abs_error = fit["u"] + fit["v"]
mean_abs_error = abs_error.sum() / fit.shape[0]
model.setObjective(mean_abs_error, sense=GRB.MINIMIZE)

Extracting Solutions

With the model formulated, we solve it using the Gurobi Optimizer:

[9]:

model.optimize()

Gurobi Optimizer version 13.0.3 build v13.0.3rc0 (linux64 - "Ubuntu 24.04.5 LTS")

CPU model: AMD EPYC 9R14, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 4 logical processors, using up to 4 threads

Optimize a model with 331 rows, 673 columns and 4303 nonzeros (Min)
Model fingerprint: 0x99b4224e
Model has 662 linear objective coefficients
Coefficient statistics:
  Matrix range     [1e-03, 4e+00]
  Objective range  [3e-03, 3e-03]
  Bounds range     [0e+00, 0e+00]
  RHS range        [2e+01, 3e+02]

Presolve time: 0.00s
Presolved: 331 rows, 673 columns, 4303 nonzeros

Iteration    Objective       Primal Inf.    Dual Inf.      Time
       0      handle free variables                          0s
     364    4.3725902e+01   0.000000e+00   0.000000e+00      0s

Solved in 364 iterations and 0.01 seconds (0.02 work units)
Optimal objective  4.372590220e+01

Using the series accessor, we can extract the result coefficients as a series and display them on a plot:

[10]:

coeffs.gppd.X.plot.bar()

[10]:

<Axes: >
../_images/examples_regression_19_1.png

We can also extract the distribution of errors from the abs_error series we constructed, and plot the results as a histogram:

[11]:

# abs_error is a Series of linear expressions
abs_error.head()

[11]:

16       u[16] + v[16]
408    u[408] + v[408]
432    u[432] + v[432]
316    u[316] + v[316]
3          u[3] + v[3]
dtype: object

[12]:

# .gppd.get_value() evaluates these expressions at the current solution
abs_error.gppd.get_value().plot.hist();
../_images/examples_regression_22_0.png