refactor joint L1-L2 optimization in CV code

This commit is contained in:
Brian Quistorff 2018-10-01 13:56:00 -07:00
Родитель bc069a70a5
Коммит 583c7bb86d
3 изменённых файлов: 103 добавлений и 92 удалений

Просмотреть файл

@ -4,7 +4,7 @@ from RidgeSC.fit_loo import loo_v_matrix, loo_weights, loo_score
from RidgeSC.fit_ct import ct_v_matrix, ct_weights, ct_score
# Public API
from RidgeSC.cross_validation import score_train_test, score_train_test_sorted_lambdas, CV_score
from RidgeSC.cross_validation import score_train_test, score_train_test_sorted_lambdas, CV_score, joint_penalty_optimzation
from RidgeSC.tensor import tensor
from RidgeSC.weights import weights
from RidgeSC.lambda_utils import get_max_lambda, L2_pen_guestimate

Просмотреть файл

@ -337,6 +337,42 @@ def CV_score(X,Y,
return total_score
def joint_penalty_optimzation(X, Y, L1_pen_start, L2_pen_start, bounds, X_treat = None, Y_treat = None):
from scipy.optimize import fmin_l_bfgs_b, differential_evolution
import time
# -----------------------------------------------------------------
# Optimization of the L2 and L1 Penalties Simultaneously, keeping their
# product constant. Heuristically, this has been most efficient means of
# optimizing the L2 Parameter.
# -----------------------------------------------------------------
# build the objective function to be minimized
# cache for L2_obj_func
n_calls = [0,]
temp_results =[]
def L1_L2_obj_func (x):
n_calls[0] += 1
t1 = time.time();
score = CV_score(X = X, Y = Y,
X_treat = X_treat, Y_treat = Y_treat,
# if LAMBDA is a single value, we get a single score, If it's an array of values, we get an array of scores.
LAMBDA = L1_pen_start * np.exp(x[0]),
L2_PEN_W = L2_pen_start / np.exp(x[0]),
# suppress the analysis type message
quiet = True)
t2 = time.time();
temp_results.append((n_calls[0],x,score))
print("calls: %s, time: %0.4f, x0: %0.4f, Cross Validation Error: %s" % (n_calls[0], t2 - t1, x[0], score))
#print("calls: %s, time: %0.4f, x0: %0.4f, x1: %0.4f, Cross Validation Error: %s, R-Squared: %s" % (n_calls[0], t2 - t1, x[0], x[1], score, 1 - score / SS ))
return score
# the actual optimization
diff_results = differential_evolution(L1_L2_obj_func, bounds = bounds)
return diff_results
# ------------------------------------------------------------
# utilities for maintaining a worker pool

Просмотреть файл

@ -288,122 +288,97 @@ if __name__ == "__main__":
# product constant. Heuristically, this has been most efficient means of
# optimizing the L2 Parameter.
# -----------------------------------------------------------------
if False:
# build the objective function to be minimized
# cache for L2_obj_func
n_calls = [0,]
temp_results =[]
SS = np.power(Y_pre_treated - np.mean(Y_pre_treated),2).sum()
best_L1_penalty_ct = np.float(1908.9329)
# build the objective function to be minimized
def L1_L2_const_obj_func (x):
n_calls[0] += 1
t1 = time.time();
score = SC.CV_score(X = X_control,
Y = Y_pre_control,
X_treat = X_treated,
Y_treat = Y_pre_treated,
# if LAMBDA is a single value, we get a single score, If it's an array of values, we get an array of scores.
LAMBDA = best_L1_penalty_ct * np.exp(x[0]),
L2_PEN_W = L2_pen_start_ct / np.exp(x[0]),
# suppress the analysis type message
quiet = True)
t2 = time.time();
temp_results.append((n_calls[0],x,score))
print("calls: %s, time: %0.4f, x0: %0.4f, Cross Validation Error: %s, out-of-sample R-Squared: %s" % (n_calls[0], t2 - t1, x[0], score, 1 - score / SS ))
#print("calls: %s, time: %0.4f, x0: %0.4f, x1: %0.4f, Cross Validation Error: %s, R-Squared: %s" % (n_calls[0], t2 - t1, x[0], x[1], score, 1 - score / SS ))
return score
# cache for L2_obj_func
n_calls = [0,]
temp_results =[]
SS = np.power(Y_pre_treated - np.mean(Y_pre_treated),2).sum()
best_L1_penalty_ct = np.float(1908.9329)
# the actual optimization
print("Starting L1-L2 Penalty (product constant) optimization")
def L1_L2_obj_func (x):
n_calls[0] += 1
t1 = time.time();
score = SC.CV_score(X = X_control,
Y = Y_pre_control,
X_treat = X_treated,
Y_treat = Y_pre_treated,
# if LAMBDA is a single value, we get a single score, If it's an array of values, we get an array of scores.
LAMBDA = best_L1_penalty_ct * np.exp(x[0]),
L2_PEN_W = L2_pen_start_ct / np.exp(x[0]),
# suppress the analysis type message
quiet = True)
t2 = time.time();
temp_results.append((n_calls[0],x,score))
print("calls: %s, time: %0.4f, x0: %0.4f, Cross Validation Error: %s, out-of-sample R-Squared: %s" % (n_calls[0], t2 - t1, x[0], score, 1 - score / SS ))
#print("calls: %s, time: %0.4f, x0: %0.4f, x1: %0.4f, Cross Validation Error: %s, R-Squared: %s" % (n_calls[0], t2 - t1, x[0], x[1], score, 1 - score / SS ))
return score
#results = differential_evolution(L1_L2_const_obj_func, bounds = ((-6,6,),) )
#results = differential_evolution(L1_L2_obj_func, bounds = ((-6,6,),)*2 )
# the actual optimization
print("Starting L2 Penalty optimization")
import pdb; pdb.set_trace()
NEW_best_L1_penalty_ct = best_L1_penalty_ct * np.exp(results.x[0])
best_L2_penalty = L2_pen_start_ct * np.exp(results.x[1])
results = differential_evolution(L1_L2_obj_func, bounds = ((-6,6,),) )
#results = differential_evolution(L1_L2_obj_func, bounds = ((-6,6,),)*2 )
import pdb; pdb.set_trace()
NEW_best_L1_penalty_ct = best_L1_penalty_ct * np.exp(results.x[0])
best_L2_penalty = L2_pen_start_ct * np.exp(results.x[1])
print("DE optimized L2 Penalty: %s, DE optimized L1 penalty: %s" % (NEW_best_L1_penalty_ct, best_L2_penalty,) )
print("DE optimized L2 Penalty: %s, DE optimized L1 penalty: %s" % (NEW_best_L1_penalty_ct, best_L2_penalty,) )
# -----------------------------------------------------------------
# Optimization of the L2 Parameter alone
# -----------------------------------------------------------------
if False:
# -----------------------
# build the objective function to be minimized
# -----------------------
# -----------------------
# build the objective function to be minimized
# -----------------------
# OBJECTIVE FUNCTION(S) TO MINIMIZE USING DE
# OBJECTIVE FUNCTION(S) TO MINIMIZE USING DE
# cache for L2_obj_func
n_calls = [0,]
temp_results =[]
SS = np.power(Y_pre_treated - np.mean(Y_pre_treated),2).sum()
best_L1_penalty_ct = np.float(1908.9329)
# cache for L2_obj_func
n_calls = [0,]
temp_results =[]
SS = np.power(Y_pre_treated - np.mean(Y_pre_treated),2).sum()
best_L1_penalty_ct = np.float(1908.9329)
def L2_obj_func (x):
n_calls[0] += 1
t1 = time.time();
def L2_obj_func (x):
n_calls[0] += 1
t1 = time.time();
score = SC.CV_score(X = X_control,
Y = Y_pre_control,
X_treat = X_treated,
Y_treat = Y_pre_treated,
# if LAMBDA is a single value, we get a single score, If it's an array of values, we get an array of scores.
LAMBDA = best_L1_penalty_ct,
L2_PEN_W = L2_pen_start_ct * np.exp(x[0]),
# suppress the analysis type message
quiet = True)
score = SC.CV_score(X = X_control,
Y = Y_pre_control,
X_treat = X_treated,
Y_treat = Y_pre_treated,
# if LAMBDA is a single value, we get a single score, If it's an array of values, we get an array of scores.
LAMBDA = best_L1_penalty_ct,
L2_PEN_W = L2_pen_start_ct * np.exp(x[0]),
# suppress the analysis type message
quiet = True)
t2 = time.time();
temp_results.append((n_calls[0],x,score))
print("calls: %s, time: %0.4f, x0: %0.4f, Cross Validation Error: %s, out-of-sample R-Squared: %s" %
(n_calls[0], t2 - t1, x[0], score, 1 - score / SS ))
return score
t2 = time.time();
temp_results.append((n_calls[0],x,score))
print("calls: %s, time: %0.4f, x0: %0.4f, Cross Validation Error: %s, out-of-sample R-Squared: %s" %
(n_calls[0], t2 - t1, x[0], score, 1 - score / SS ))
return score
# the actual optimization
print("Starting L2 Penalty optimization")
# the actual optimization
print("Starting L2 Penalty optimization")
results = differential_evolution(L2_obj_func, bounds = ((-6,6,),))
NEW_L2_pen_start_ct = L2_pen_start_ct * np.exp(results.x[1])
print("DE optimized L2 Penalty: %s, using fixed L1 penalty: %s" % (NEW_L2_pen_start_ct, best_L1_penalty_ct,) )
results = differential_evolution(L2_obj_func, bounds = ((-6,6,),))
NEW_L2_pen_start_ct = L2_pen_start_ct * np.exp(results.x[0])
print("DE optimized L2 Penalty: %s, using fixed L1 penalty: %s" % (NEW_L2_pen_start_ct, best_L1_penalty_ct,) )
# -----------------------------------------------------------------
# Optimization of the L2 and L1 Penalties Simultaneously
# -----------------------------------------------------------------
# build the objective function to be minimized
# cache for L2_obj_func
n_calls = [0,]
temp_results =[]
SS = np.power(Y_pre_treated - np.mean(Y_pre_treated),2).sum()
best_L1_penalty_ct = np.float(1908.9329)
def L1_L2_obj_func (x):
n_calls[0] += 1
t1 = time.time();
score = SC.CV_score(X = X_control,
Y = Y_pre_control,
X_treat = X_treated,
Y_treat = Y_pre_treated,
# if LAMBDA is a single value, we get a single score, If it's an array of values, we get an array of scores.
LAMBDA = best_L1_penalty_ct * np.exp(x[0]),
L2_PEN_W = L2_pen_start_ct * np.exp(x[1]),
# suppress the analysis type message
quiet = True)
t2 = time.time();
temp_results.append((n_calls[0],x,score))
print("calls: %s, time: %0.4f, x0: %0.4f, x1: %0.4f, Cross Validation Error: %s, out-of-sample R-Squared: %s" % (n_calls[0], t2 - t1, x[0], x[1], score, 1 - score / SS ))
return score
# the actual optimization
print("Starting L2 Penalty optimization")
print("Starting L1-L2 Joint Penalty optimization")
results = differential_evolution(L1_L2_obj_func, bounds = ((-6,6,),)*2 )
results = SC.joint_penalty_optimzation(X = X_control, Y = Y_pre_control, L1_pen_start = best_L1_penalty_ct, L2_pen_start = L2_pen_start_ct, bounds = ((-6,6,),)*2, X_treat = X_treated, Y_treat = Y_pre_treated)
import pdb; pdb.set_trace()
NEW_best_L1_penalty_ct = best_L1_penalty_ct * np.exp(results.x[0])