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.
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.
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]]