test pure glmnet cd python implementation against cd_fast.enet_coordinate_descent

cd_fast2.py

import numpy as np


def fsign( f):
    if f == 0:
        return 0
    elif f > 0:
        return 1.0
    else:
        return -1.0


def enet_coordinate_descent2(w, alpha, beta, X, y, max_iter):

    n_samples = X.shape[0]
    n_features = X.shape[1]

    norm_cols_X = (X ** 2).sum(axis=0)

    R = y - np.dot(X,w)

    for n_iter in range(max_iter):

        for ii in xrange(n_features):

            w_ii = w[ii]

            # Naive Updates
            partial_residual = y - np.inner(X, w) + X[:, ii] * w_ii
            tmp = np.dot(X[:, ii], partial_residual)

            # Covariance Updates
            #problem_sum = np.array([np.dot(X[:,ii], X[:,k])*w[k] for k in xrange(n_features)]).sum()
            #tmp = np.dot(X[:,ii], y) - problem_sum +  w_ii * n_samples

            w[ii] = fsign(tmp) * max(abs(tmp) - alpha, 0) \
                / (norm_cols_X[ii] + beta)

    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
    alpha = alpha * rho * n_samples
    beta = alpha * (1.0 - rho) * n_samples
    w = np.zeros(n_features)
    X = np.asfortranarray(X) 
    result_org = cd_fast.enet_coordinate_descent(w, alpha, beta,
                                X, y, max_iter=1000, tol=1e-4, positive=False)
    w = np.zeros(n_features)
    #print result_org
    my_result = enet_coordinate_descent2(w, alpha, beta, X, y, max_iter=1000)
    assert_array_almost_equal(my_result,
            np.array([0.52323384, -0.00908868]),4)


def test_big_data():
    X, y = make_regression(n_samples=10, n_features=20, n_informative=10,
                    random_state=0)
    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.zeros(n_features)
    X = np.asfortranarray(X) 
    result_org = cd_fast.enet_coordinate_descent(w, alpha, beta,
                                X, y, max_iter=1000, tol=1e-4, 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,
    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]),2)