/
nonneg_l2.py
116 lines (99 loc) · 3.78 KB
/
nonneg_l2.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
#!/usr/bin/env python
import scipy.sparse as spspa
import scipy as sp
import numpy as np
import quadprog.problem as qp
from quadprog.solvers.solvers import GUROBI, CPLEX, OSQP
import quadprog.solvers.osqp.osqp as osqp
# Generate and solve a nonnegative least-squares problem
class nonneg_l2(object):
"""
Nonnegative least-squares problem is defined as
minimize || Ax - b ||^2
subjec to x >= 0
Arguments
---------
m, n - Dimensions of matrix A <int>
osqp_opts - Parameters of OSQP solver
dens_lvl - Density level of matrix A <float>
version - QP reformulation ['dense', 'sparse']
"""
def __init__(self, m, n, dens_lvl=1.0, version='dense',
osqp_opts={}):
# Generate data
Ad = spspa.random(m, n, density=dens_lvl, format='csc')
x_true = np.ones(n) / n + np.random.randn(n) / np.sqrt(n)
bd = Ad.dot(x_true) + 0.5*np.random.randn(m)
# Construct the problem
if version == 'dense':
# minimize 1/2 x.T (A.T * A) x - (A.T * b).T * x
# subject to x >= 0
P = Ad.T.dot(Ad)
q = -Ad.T.dot(bd)
A = spspa.eye(n)
lA = np.zeros(n)
uA = np.inf * np.ones(n)
elif version == 'sparse':
# minimize 1/2 y.T*y
# subject to y = Ax - b
# x >= 0
Im = spspa.eye(m)
P = spspa.block_diag((spspa.csc_matrix((n, n)), Im), format='csc')
q = np.zeros(n + m)
A = spspa.vstack([
spspa.hstack([Ad, -Im]),
spspa.hstack([spspa.eye(n), spspa.csc_matrix((n, m))]),
]).tocsc()
lA = np.append(bd, np.zeros(n))
uA = np.append(bd, np.inf * np.ones(n))
else:
assert False, "Unhandled version"
# Create a quadprogProblem and store it in a private variable
self._prob = qp.quadprogProblem(P, q, A, lA, uA)
# Create an OSQP object and store it in a private variable
self._osqp = osqp.OSQP(**osqp_opts)
self._osqp.problem(P, q, A, lA, uA)
def solve(self, solver=OSQP):
"""
Solve the problem with a specificed solver.
"""
if solver == OSQP:
results = self._osqp.solve()
elif solver == CPLEX:
results = self._prob.solve(solver=CPLEX, verbose=1)
elif solver == GUROBI:
results = self._prob.solve(solver=GUROBI, OutputFlag=0)
else:
assert False, "Unhandled solver"
return results
# ==================================================================
# == Solve a small example
# ==================================================================
# Set the random number generator
sp.random.seed(1)
# Set problem
m = 3
n = 2*m
# Set options of the OSQP solver
options = {'eps_abs': 1e-4,
'eps_rel': 1e-4,
'alpha': 1.6,
'scale_problem': True,
'scale_steps': 4,
'polish': False}
# Create a lasso object
nonneg_l2_obj = nonneg_l2(m, n, dens_lvl=0.5, version='sparse',
osqp_opts=options)
# Solve with different solvers
resultsCPLEX = nonneg_l2_obj.solve(solver=CPLEX)
resultsGUROBI = nonneg_l2_obj.solve(solver=GUROBI)
resultsOSQP = nonneg_l2_obj.solve(solver=OSQP)
# Print objective values
print "CPLEX Objective Value: %.3f" % resultsCPLEX.objval
print "GUROBI Objective Value: %.3f" % resultsGUROBI.objval
print "OSQP Objective Value: %.3f" % resultsOSQP.objval
print "\n"
# Print timings
print "CPLEX CPU time: %.3f" % resultsCPLEX.cputime
print "GUROBI CPU time: %.3f" % resultsGUROBI.cputime
print "OSQP CPU time: %.3f" % resultsOSQP.cputime