How do I model this linear programming problem in Python?

Question

I have been tasked to program a Dantzig Selector using Python, but I was given no guidelines and do not have much experience in linear programming or data science. I cannot find the information I need in LP module manuals, or in other questions on this site.

This is the problem. I am looking for the column vector ˆβ. Sorry that I do not have any code for this part of my program, as I don't know how to approach this problem. I have tried several approaches, but none correctly reflected the problem, so I rejected and deleted them.

min||ˆβ||l1 s.t. ||xT(y-xˆβ(||l(inf) <= δ

It can be rewritten as a linear program.

ˆβ is a kx1 column vector and the Dantzig Selector I am looking for.

• y is a nx1 column vector of observations/responses
• X is a nxk sample matrix, where k >> n
• δ is a noise variable

Here are more details that may be useful.

Here is my working code so far. Data values are all just samples/placeholders. I have already prepared X, y, and some δ values. However, I cannot find the right LP function to give me ˆβ.

import numpy as np
import random
import math

#n = no. runs = 5
n = 5
#k = no. variables = 23
k = 23

#y = vector of observations/responses (nx1, binary decisions)
y = np.array([[1],
              [0],
              [0],
              [1],
              [0]])
#X = predictor/sample matrix (nxk)
X = np.array([[1.1, 0, 0.7, 0.8, 0.9, 0.2, 0.3, 0.5, 0.2, 0.2, 1.2, 1.1, 0.5, 0.5, 0.7, 1.2, 1.3, 0.8, 0.9, 1.7, 1.2, 1.9, 0.9],
              [0.3, 0.1, 0.7, 0.4, 0.9, 0.9, 0.1, 0.8, 0.1, 0.2, 1.1, 0, 0.9, 0.4, 1.4, 1.4, 0.1, 0.5, 1.8, 1.6, 1.2, 1.8, 0.3],
              [0.1, 0.1, 0.3, 0.9, 0.7, 0.8, 0, 0.7, 0.8, 0.2, 1.1, 1.1, 0.5, 0.5, 0.8, 1.5, 0.2, 0.5, 1.6, 1.5, 1.2, 1.7, 0.5],
              [1.2, 0.2, 0.9, 0.8, 0.6, 0.2, 0.3, 0.5, 0.3, 0.2, 1.2, 1.1, 0.5, 0, 0.7, 1.2, 1.3, 0.8, 0.9, 1.7, 1.2, 1.9, 0.9],
              [0.2, 0.1, 0.6, 0, 0.5, 1.1, 0.2, 0.5, 0.9, 0.2, 1.2, 1.1, 0.8, 1.6, 0.5, 1.3, 0.2, 0.5, 1.7, 1.2, 1.2, 1.9, 0.1]])


#estimate missing data (0)
X_row_minima = np.where(X>0,X,X.max()).min(1)
X[X==0] = X_row_minima/2

#unit length normalize X
X = X/np.linalg.norm(X, ord=2, axis=1, keepdims=True)

#standardize y to zero mean
y = y - np.mean(y) / np.std(y)

#transpose X (kxn)
Xt = np.transpose(X)

#solve d0
Xty = np.matmul(Xt,y)
d0 = max(abs(Xty))

#generate 100 evenly-spaced d values
d = np.linspace(0, d0, 100)

This is my first post on this site. I apologize for the lack of details in the post compared to others.

Answer 1

I am not sure if I have this right, as I haven't had to figure out a Danzig Selector before. If possible, I would suggest testing with a dataset where you know what to expect for an answer. From what I can tell, your problem is something like this.

minimize:   sum(u)
subject to: -u_i <= beta_i for all i in len(u)
            u_i >= beta_i for all i in len(u)
            [X^T (y - X beta)]_j <= delta_j for all j in len(y)
            [X^T (y - X beta)]_j >= -delta_j for all j in len(y)

Writing this in cvxpy is pretty straightforward, I think. If I'm off then hopefully it gets you in the right direction. Using your values for X, y, k:

import cvxpy as cp
beta = cp.Variable((k,1))
u = cp.Variable((k,1))
delta = 1

objective = cp.Minimize(cp.sum(u))

constraints = [beta >= -u,
               beta <= u,
               Xt @ (y - X @ beta) <= delta,
               Xt @ (y - X @ beta) >= -delta]
               

prob = cp.Problem(objective, constraints)
prob.solve(verbose=True)

print("Problem status: ", prob.status)
print("Optimal value", prob.value)
print("Beta values", beta.value)

My resulting output:

===============================================================================
                                     CVXPY
                                    v1.1.15
===============================================================================
(CVXPY) Oct 12 11:00:50 PM: Your problem has 46 variables, 4 constraints, and 0 parameters.
(CVXPY) Oct 12 11:00:50 PM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Oct 12 11:00:50 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Oct 12 11:00:50 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
-------------------------------------------------------------------------------
                                  Compilation
-------------------------------------------------------------------------------
(CVXPY) Oct 12 11:00:50 PM: Compiling problem (target solver=ECOS).
(CVXPY) Oct 12 11:00:50 PM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> ECOS
(CVXPY) Oct 12 11:00:50 PM: Applying reduction Dcp2Cone
(CVXPY) Oct 12 11:00:50 PM: Applying reduction CvxAttr2Constr
(CVXPY) Oct 12 11:00:50 PM: Applying reduction ConeMatrixStuffing
(CVXPY) Oct 12 11:00:50 PM: Applying reduction ECOS
(CVXPY) Oct 12 11:00:50 PM: Finished problem compilation (took 3.027e-03 seconds).
-------------------------------------------------------------------------------
                                Numerical solver
-------------------------------------------------------------------------------
(CVXPY) Oct 12 11:00:50 PM: Invoking solver ECOS  to obtain a solution.

ECOS 2.0.7 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS

It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
 0  +0.000e+00  -1.596e+02  +3e+02  2e-01  6e-01  1e+00  4e+00    ---    ---    1  1  - |  -  -
 1  +4.462e+00  -6.924e+00  +4e+01  1e-02  6e-02  2e-01  4e-01  0.9228  4e-02   0  0  0 |  0  0
 2  +4.980e-01  -8.219e-01  +4e+00  1e-03  8e-03  3e-02  5e-02  0.9010  3e-02   0  0  0 |  0  0
 3  +3.673e-02  -4.271e-02  +3e-01  9e-05  5e-04  1e-03  3e-03  0.9630  2e-02   0  0  0 |  0  0
 4  +4.054e-04  -4.749e-04  +3e-03  1e-06  5e-06  1e-05  3e-05  0.9890  1e-04   0  0  0 |  0  0
 5  +4.490e-06  -5.272e-06  +3e-05  1e-08  6e-08  2e-07  4e-07  0.9890  1e-04   0  0  0 |  0  0
 6  +4.972e-08  -5.852e-08  +4e-07  1e-10  6e-10  2e-09  4e-09  0.9890  1e-04   0  0  0 |  0  0
 7  +5.507e-10  -6.496e-10  +4e-09  1e-12  7e-12  2e-11  4e-11  0.9890  1e-04   0  0  0 |  0  0

OPTIMAL (within feastol=6.9e-12, reltol=-nan(ind), abstol=4.1e-09).
Runtime: 0.002074 seconds.

-------------------------------------------------------------------------------
                                    Summary
-------------------------------------------------------------------------------
(CVXPY) Oct 12 11:00:50 PM: Problem status: optimal
(CVXPY) Oct 12 11:00:50 PM: Optimal value: 5.507e-10
(CVXPY) Oct 12 11:00:50 PM: Compilation took 3.027e-03 seconds
(CVXPY) Oct 12 11:00:50 PM: Solver (including time spent in interface) took 2.358e-03 seconds
Problem status:  optimal
Optimal value 5.507030257315668e-10
Beta values [[ 1.43866286e-11]
 [ 2.52691994e-12]
 [ 1.26496527e-11]
 [ 1.25889902e-11]
 [ 1.32324154e-11]
 [ 9.34898130e-12]
 [ 4.19460380e-12]
 [ 1.07283774e-11]
 [ 7.54676427e-12]
 [ 3.94492714e-12]
 [ 1.68905798e-11]
 [ 1.61644206e-11]
 [ 1.09177395e-11]
 [ 9.11945636e-12]
 [ 1.38073390e-11]
 [ 1.51133406e-11]
 [ 1.59259478e-11]
 [ 1.24992988e-11]
 [ 1.10965990e-11]
 [ 1.50106084e-11]
 [ 1.66261028e-11]
 [-1.92542442e-13]
 [ 1.24499858e-11]]