Skip to content

Commit f5bcd54

Browse files
idc9idc9
committed
scipy lin prog, cvxpy quad prog
1 parent 46ace1f commit f5bcd54

13 files changed

Lines changed: 466 additions & 312 deletions

‎ya_glm/add_init_params.py‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,7 @@ def init_wrapper(init):
2828
# start with init's current parameters
2929
init_params = list(signature(init).parameters.values())
3030
init_params = init_params[1:] # ignore self
31-
init_param_names = set(p.name for p in init_params)
31+
current_param_names = set(p.name for p in init_params)
3232

3333
empty_init_params = set(['self', 'args', 'kwargs'])
3434

@@ -48,7 +48,8 @@ def init_wrapper(init):
4848
cls_params = cls_params[1:] # ignore self
4949
# ignore parameter if it was already in init
5050
cls_params = [p for p in cls_params
51-
if p.name not in init_param_names]
51+
if p.name not in current_param_names]
52+
current_param_names.update([p.name for p in cls_params])
5253

5354
params.extend(cls_params)
5455

‎ya_glm/backends/cvxpy/glm_solver.py‎

Lines changed: 29 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
import cvxpy as cp
22
from functools import partial
3+
from time import time
34

45
from ya_glm.utils import clip_zero
56
from ya_glm.cvxpy.penalty import lasso, ridge
@@ -26,6 +27,7 @@ def solve_glm(X, y,
2627
zero_tol=1e-8,
2728
cp_kws={}):
2829

30+
start_time = time()
2931
######################
3032
# objective function #
3133
######################
@@ -50,10 +52,15 @@ def solve_glm(X, y,
5052
if coef.value is None:
5153
raise RuntimeError("cvxpy solvers failed")
5254

53-
return process_output(problem=problem, coef=coef,
54-
intercept=intercept,
55-
fit_intercept=fit_intercept,
56-
zero_tol=zero_tol)
55+
coef, intercept, opt_data = \
56+
process_output(problem=problem, coef=coef,
57+
intercept=intercept,
58+
fit_intercept=fit_intercept,
59+
zero_tol=zero_tol)
60+
61+
opt_data['runtime'] = time() - start_time
62+
63+
return coef, intercept, opt_data
5764

5865

5966
def solve_glm_path(fit_intercept=True, cp_kws={}, zero_tol=1e-8,
@@ -65,14 +72,18 @@ def solve_glm_path(fit_intercept=True, cp_kws={}, zero_tol=1e-8,
6572
check_decr=check_decr)
6673

6774
# make sure we setup the right penalty
68-
if 'lasso_pen' in param_path:
69-
kws['lasso_pen'] = param_path['lasso_pen'][0]
70-
if 'ridge_pen' in param_path:
71-
kws['ridge_pen'] = param_path['ridge_pen'][0]
75+
if 'lasso_pen' in param_path[0]:
76+
kws['lasso_pen'] = param_path[0]['lasso_pen']
77+
if 'ridge_pen' in param_path[0]:
78+
kws['ridge_pen'] = param_path[0]['ridge_pen']
7279

80+
start_time = time()
7381
problem, coef, intercept, lasso_pen, ridge_pen = setup_problem(**kws)
82+
pre_setup_runtime = start_time - time()
7483

7584
for params in param_path:
85+
start_time = time()
86+
7687
if 'lasso_pen' in params:
7788
lasso_pen.value = params['lasso_pen']
7889

@@ -90,7 +101,13 @@ def solve_glm_path(fit_intercept=True, cp_kws={}, zero_tol=1e-8,
90101
fit_intercept=fit_intercept,
91102
zero_tol=zero_tol)
92103

93-
# if generator:
104+
fit_out = {'coef': fit_out[0],
105+
'intercept': fit_out[1],
106+
'opt_data': fit_out[2]}
107+
108+
fit_out['opt_data']['runtime'] = start_time - time()
109+
fit_out['opt_data']['pre_setup_runtime'] = pre_setup_runtime
110+
94111
yield fit_out, params
95112

96113

@@ -174,19 +191,16 @@ def objective(coef, intercept):
174191
def objective(coef, intercept):
175192
return glm_loss(X=X, y=y, coef=coef, intercept=intercept)
176193

177-
###################
178-
# setup variables #
179-
###################
194+
###############################
195+
# setup variables and problem #
196+
###############################
180197

181198
coef = cp.Variable(shape=X.shape[1], value=coef_init)
182199
if fit_intercept:
183200
intercept = cp.Variable(value=intercept_init)
184201
else:
185202
intercept = None
186203

187-
###########################
188-
# setup and solve problem #
189-
###########################
190204
problem = cp.Problem(cp.Minimize(objective(coef, intercept)))
191205

192206
return problem, coef, intercept, lasso_pen, ridge_pen

ya_glm/backends/scipy/__init__.py renamed to ya_glm/backends/quantile_lp/__init__.py

File renamed without changes.

‎ya_glm/backends/quantile_lp/cvxpy_quad_prog.py‎

Lines changed: 173 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,173 @@
1+
import cvxpy as cp
2+
from time import time
3+
4+
from ya_glm.backends.fista.glm_solver import process_param_path
5+
from ya_glm.backends.quantile_lp.utils import get_lin_prog_data, \
6+
get_coef_inter, get_quad_mat
7+
8+
9+
def solve(X, y, fit_intercept=True, quantile=0.5, sample_weights=None,
10+
lasso_pen=1, ridge_pen=None,
11+
lasso_weights=None, ridge_weights=None, tikhonov=None,
12+
coef_init=None, intercept_init=None,
13+
solver=None,
14+
cp_kws={}):
15+
"""
16+
Solves the L1 + L2 penalized quantile regression problem by formulating it as a linear quadratic program then appealing to cvxpy.
17+
"""
18+
19+
if lasso_weights is not None and lasso_pen is None:
20+
lasso_pen = 1
21+
22+
if (ridge_weights is not None or tikhonov is not None) \
23+
and ridge_pen is None:
24+
ridge_pen = 1
25+
26+
start_time = time()
27+
28+
problem, var, lasso_pen, ridge_pen = \
29+
setup_problem(X=X, y=y,
30+
fit_intercept=fit_intercept,
31+
quantile=quantile,
32+
sample_weights=sample_weights,
33+
lasso_pen=lasso_pen,
34+
ridge_pen=ridge_pen,
35+
lasso_weights=lasso_weights,
36+
ridge_weights=ridge_weights,
37+
tikhonov=tikhonov,
38+
coef_init=coef_init,
39+
intercept_init=intercept_init)
40+
41+
problem.solve(solver=solver, **cp_kws)
42+
# solve_with_backups(problem=problem, variable=var, **cp_kws)
43+
44+
opt_data = {**problem.solver_stats.__dict__,
45+
'status': problem.status,
46+
'runtime': time() - start_time}
47+
48+
if fit_intercept:
49+
n_params = X.shape[1] + 1
50+
else:
51+
n_params = X.shape[1]
52+
53+
coef, intercept = get_coef_inter(solution=var.value,
54+
n_params=n_params,
55+
fit_intercept=fit_intercept)
56+
57+
# coef = clip_zero(coef, zero_tol=zero_tol)
58+
# if fit_intercept:
59+
# intercept = clip_zero(intercept, zero_tol=zero_tol)
60+
# else:
61+
# intercept = None
62+
63+
return coef, intercept, opt_data
64+
65+
66+
def solve_path(fit_intercept=True, cp_kws={}, zero_tol=1e-8,
67+
lasso_pen_seq=None, ridge_pen_seq=None,
68+
check_decr=True, **kws):
69+
70+
param_path = process_param_path(lasso_pen_seq=lasso_pen_seq,
71+
ridge_pen_seq=ridge_pen_seq,
72+
check_decr=check_decr)
73+
74+
# make sure we setup the right penalty
75+
if 'lasso_pen' in param_path[0]:
76+
kws['lasso_pen'] = param_path[0]['lasso_pen']
77+
if 'ridge_pen' in param_path[0]:
78+
kws['ridge_pen'] = param_path[0]['ridge_pen']
79+
80+
start_time = time()
81+
problem, var, lasso_pen, ridge_pen = setup_problem(**kws)
82+
pre_setup_runtime = time() - start_time
83+
84+
for params in param_path:
85+
start_time = time()
86+
87+
if 'lasso_pen' in params:
88+
lasso_pen.value = params['lasso_pen']
89+
90+
if 'ridge_pen' in params:
91+
ridge_pen.value = params['ridge_pen']
92+
93+
problem.solve(**cp_kws)
94+
# solve_with_backups(problem=problem, variable=var, **cp_kws)
95+
96+
if var.value is None:
97+
raise RuntimeError("cvxpy solvers failed")
98+
99+
opt_data = {**problem.solver_stats.__dict__,
100+
'status': problem.status,
101+
'runtime': time() - start_time,
102+
'pre_setup_runtime': pre_setup_runtime}
103+
104+
if fit_intercept:
105+
n_params = kws['X'].shape[1] + 1
106+
else:
107+
n_params = kws['X'].shape[1]
108+
109+
coef, intercept = get_coef_inter(solution=var.value,
110+
n_params=n_params,
111+
fit_intercept=fit_intercept)
112+
113+
# coef = clip_zero(coef, zero_tol=zero_tol)
114+
# if fit_intercept:
115+
# intercept = clip_zero(intercept, zero_tol=zero_tol)
116+
# else:
117+
# intercept = None
118+
119+
fit_out = {'coef': coef, 'intercept': intercept, 'opt_data': opt_data}
120+
yield fit_out, params
121+
122+
123+
def setup_problem(X, y, fit_intercept=True, quantile=0.5, sample_weights=None,
124+
lasso_pen=1, ridge_pen=None,
125+
lasso_weights=None, ridge_weights=None, tikhonov=None,
126+
coef_init=None, intercept_init=None):
127+
128+
if lasso_pen is not None:
129+
lasso_pen = cp.Parameter(pos=True, value=lasso_pen)
130+
131+
if ridge_pen is not None:
132+
ridge_pen = cp.Parameter(pos=True, value=ridge_pen)
133+
134+
if coef_init is not None or intercept_init is not None:
135+
raise NotImplementedError("I do not think initialization works for these solvers")
136+
137+
######################
138+
# setup problem data #
139+
######################
140+
A_eq, b_eq, lin_coef, n_params = \
141+
get_lin_prog_data(X, y,
142+
fit_intercept=fit_intercept,
143+
quantile=quantile,
144+
lasso_pen=lasso_pen,
145+
sample_weights=sample_weights,
146+
lasso_weights=lasso_weights)
147+
148+
lin_coef = cp.hstack(lin_coef)
149+
150+
if ridge_pen is not None:
151+
quad_mat = get_quad_mat(X=X,
152+
fit_intercept=fit_intercept,
153+
weights=ridge_weights,
154+
tikhonov=tikhonov)
155+
156+
n_dim = A_eq.shape[1]
157+
var = cp.Variable(shape=n_dim)
158+
159+
####################
160+
# setup cp problem #
161+
####################
162+
if ridge_pen is None:
163+
objective = cp.Minimize(var.T @ lin_coef)
164+
else:
165+
objective = cp.Minimize(var.T @ lin_coef +
166+
0.5 * ridge_pen * cp.quad_form(var, quad_mat))
167+
168+
constraints = [var >= 0,
169+
A_eq @ var == b_eq]
170+
171+
problem = cp.Problem(objective, constraints)
172+
173+
return problem, var, lasso_pen, ridge_pen

‎ya_glm/backends/quantile_lp/glm_solver.py‎

Lines changed: 83 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,83 @@
1+
from ya_glm.backends.quantile_lp.scipy_lin_prog import solve as solve_lin_prog
2+
from ya_glm.backends.quantile_lp.cvxpy_quad_prog import solve as solve_quad_prog
3+
from ya_glm.backends.quantile_lp.cvxpy_quad_prog import solve_path
4+
5+
6+
def solve_glm(X, y,
7+
loss_func='quantile',
8+
loss_kws={},
9+
fit_intercept=True,
10+
11+
lasso_pen=None,
12+
lasso_weights=None,
13+
14+
ridge_pen=None,
15+
ridge_weights=None,
16+
tikhonov=None,
17+
18+
coef_init=None,
19+
intercept_init=None,
20+
21+
solver='default',
22+
solver_kws={}
23+
):
24+
"""
25+
Solves quantile regression with either a Linear Programming (LP) formulation (for unpenalizer or Lasso penalties) or a Quadratic Programming (QP) formulation (for ridge penalties). For LPs we uses scipy's linprog solver. For QPs we use cvxpy.
26+
27+
28+
29+
"""
30+
if loss_func != 'quantile':
31+
raise NotImplementedError("This solver only works for quantile regression")
32+
33+
if coef_init is not None or intercept_init is not None:
34+
raise NotImplementedError("I do not think initialization works for these solvers")
35+
36+
if lasso_weights is not None and lasso_pen is None:
37+
lasso_pen = 1
38+
39+
if (ridge_weights is not None or tikhonov is not None) \
40+
and ridge_pen is None:
41+
ridge_pen = 1
42+
43+
quantile = loss_kws['quantile']
44+
45+
kws = {'X': X,
46+
'y': y,
47+
'fit_intercept': fit_intercept,
48+
'quantile': quantile,
49+
'lasso_pen': lasso_pen,
50+
'lasso_weights': lasso_weights,
51+
# 'sample_weights': None, # TODO: add
52+
**solver_kws}
53+
54+
if ridge_pen is None:
55+
if solver == 'default':
56+
solver = 'highs'
57+
58+
return solve_lin_prog(solver=solver,
59+
**kws)
60+
61+
else:
62+
if solver == 'default':
63+
solver = 'ECOS'
64+
65+
return solve_quad_prog(**kws,
66+
ridge_pen=ridge_pen,
67+
ridge_weights=ridge_weights,
68+
tikhonov=tikhonov,
69+
solver=solver,
70+
cp_kws=solver_kws)
71+
72+
73+
def solve_glm_path(loss_func='quantile',
74+
loss_kws={}, **kws):
75+
"""
76+
Path algorithm for the Linear and Quadratic Program formulations of quantile regression solved using cvxpy. This is not a true path algorithm in the sense that (I believe) the solution is not reused. However this does save time by resuing the solver setups.
77+
"""
78+
79+
if loss_func != 'quantile':
80+
raise NotImplementedError("This solver only works for quantile regression")
81+
82+
quantile = loss_kws.pop('quantile', 0.5)
83+
return solve_path(loss_func=loss_func, quantile=quantile, **kws)

0 commit comments

Comments
 (0)