-
-
Save ibayer/2984543 to your computer and use it in GitHub Desktop.
test pure glmnet cd python implementation against cd_fast.enet_coordinate_descent
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import numpy as np | |
from cd_regression import enet_f | |
def fsign(f): | |
if f == 0: | |
return 0 | |
elif f > 0: | |
return 1.0 | |
else: | |
return -1.0 | |
def check_convergence(y, X, w, value_enet_f): | |
if value_enet_f < enet_f(y, X, w, rho=0.3, alpha=0.5): | |
print "+" | |
else: | |
print "-" | |
value_enet_f = enet_f(y, X, w, rho=0.3, alpha=0.5) | |
print value_enet_f | |
return value_enet_f | |
def enet_coordinate_descent2(w, l2_reg, l1_reg, X, y, max_iter): | |
n_samples = X.shape[0] | |
n_features = X.shape[1] | |
norm_cols_X = (X ** 2).sum(axis=0) | |
Xy = np.dot(X.T,y) | |
gradient = np.zeros(n_features) | |
feature_inner_product = np.zeros(shape=(n_features, n_features)) | |
active_set = set(range(n_features)) | |
#debug | |
value_enet_f = 0 | |
for n_iter in range(max_iter): | |
for ii in active_set: | |
w_ii = w[ii] | |
# initial calculation | |
if n_iter == 0: | |
feature_inner_product[:, ii] = np.dot(X[:, ii], X) | |
gradient[ii] = Xy[ii] - np.dot(feature_inner_product[:, ii], w) | |
tmp = gradient[ii] + w_ii * norm_cols_X[ii] | |
w[ii] = fsign(tmp) * max(abs(tmp) - l2_reg, 0) \ | |
/ (norm_cols_X[ii] + l1_reg) | |
# update gradients, if coef changed | |
if w_ii != w[ii]: | |
for j in active_set: | |
if n_iter >= 1 or j <= ii: | |
gradient[j] -= feature_inner_product[ii, j] * \ | |
(w[ii] - w_ii) | |
# debug | |
#value_enet_f = check_convergence(y, X, w, value_enet_f) | |
#print value_enet_f | |
#remove inactive features | |
tmp_s = set.copy(active_set) | |
for j in tmp_s: | |
if w[j] == 0: | |
active_set.remove(j) | |
return w |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import numpy as np | |
from sklearn.linear_model import cd_fast | |
from cd_fast2 import enet_coordinate_descent2 | |
from sklearn.linear_model.coordinate_descent import ElasticNet | |
from sklearn.linear_model.base import center_data | |
from numpy.testing import assert_array_almost_equal, assert_almost_equal, \ | |
assert_equal | |
from sklearn.datasets.samples_generator import make_regression | |
# ATTENTION does not pass with w = 0 as start value | |
def test_line(): | |
X = np.array([[-1], [0.], [1.]]) | |
y = np.array([-1.0, 0.0, 1.0]) # just a straight line | |
n_samples, n_features = X.shape | |
rho = 0.3 | |
alpha = 0.5 | |
alpha = alpha * rho * n_samples | |
beta = alpha * (1.0 - rho) * n_samples | |
w = np.array([0.2]) | |
my_result = enet_coordinate_descent2(w, alpha, beta, X, y, max_iter=100) | |
assert_array_almost_equal(my_result, | |
np.array([0.52631579])) | |
# cd_fast.enet_coordinate_descent(w, alpha, beta, | |
# X, y, max_iter=100, tol=1e-4, positive=False)[0]) | |
def test_2d(): | |
X = np.array([[-1, 0.0], [0., 1.0], [1., -1.]]) | |
y = np.array([-1.0, 0.0, 1.0]) # just a straight line | |
rho = 0.3 | |
alpha = 0.5 | |
n_samples, n_features = X.shape | |
l2_reg = alpha * rho * n_samples | |
l1_reg = alpha * (1.0 - rho) * n_samples | |
w = np.zeros(n_features) | |
X = np.asfortranarray(X) | |
result_org, gap, tol = cd_fast.enet_coordinate_descent(w, l2_reg, l1_reg, | |
X, y, max_iter=100, tol=1e-7, positive=False) | |
w = np.zeros(n_features) | |
#print result_org | |
my_result = enet_coordinate_descent2(w, l2_reg, l1_reg, X, y, max_iter=100) | |
assert_array_almost_equal(my_result, result_org, 9) | |
# assert_array_almost_equal(my_result, | |
# np.array([0.52323384, -0.00908868]),7) | |
def test_active_set(): | |
# test set with zeros in solution coef | |
X, y = make_regression(n_samples=40, n_features=20, n_informative=5, | |
random_state=0) | |
n_samples, n_features = X.shape | |
rho = 0.80 | |
alpha = 10 | |
alpha = alpha * rho * n_samples | |
beta = alpha * (1.0 - rho) * n_samples | |
w = np.zeros(n_features) | |
X = np.asfortranarray(X) | |
result_org, gap, tol = cd_fast.enet_coordinate_descent(w, alpha, beta, | |
X, y, max_iter=1000, tol=1e-9, positive=False) | |
w = np.zeros(n_features) | |
my_result = enet_coordinate_descent2(w, alpha, beta, X, y, max_iter=1000) | |
print result_org[0] | |
assert_array_almost_equal(my_result, result_org, 7) | |
# np.array([38.18037338, 18.4702112, 9.86198851, -1.46801215, 16.52490931 | |
# , 14.26861543, 18.15508878, 36.40871624, 0., 12.35964046 | |
# ,6.98213445, 30.17242224,7.07032768,4.42177579, -1.73831861 | |
# , -7.26278943,0.34912212, 48.84641316,8.05922053, 10.301779]) |
Yes I know, it's basically a place holder at the moment. It's not clear to me what the best solution is, re-computation result's in low memory but high cost.
On the other side only the values for active features need to stored. But I don't know how many to expect, so how much memory do I allocate? What if the allocated memory is not enough? If I dynamical allocate memory, do I set a limit?
Yes I know, it's basically a place holder at the moment. It's not clear to me what the best solution is, re-computation result's in low memory but high cost.
for Lasso at least you know that the number of active features is less
than min(n_samples, n_features)
for Elastic-Net is below n_features which is not very helpful…
On the other side only the values for active features need to stored. But I don't know how many to expect, so how much memory do I allocate? What if the allocated memory is not enough? If I dynamical allocate memory, do I set a limit?
you can investigate different dynamic allocation schemes
in lasso_lars we use a fixed size (user specified)
https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/linear_model/least_angle.py#L112
but it might not be the optimal thing to do for your case.
HTH
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
feature_inner_product = np.zeros(shape=(n_features, n_features))
is forbidden as the point is to work with huge feature space. You must avoid this.