test pure glmnet cd python implementation against cd_fast.enet_coordinate_descent

cd_fast2.py

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

test_enet.py

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

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.

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