linear_model.py

from sklearn.preprocessing import LabelEncoder
from sklearn.utils.validation import check_random_state, check_X_y, check_array
import numpy as np
import theano.tensor as T
import theano

class LogisticRegression(object):

    floatX = theano.config.floatX

    def __init__(self, l1=0.0, l2=0.0, random_state=0, tol=1.0e-6,
                 max_iter=100, learning_rate=0.5):

        # Parameters
        self.l1 = l1
        self.l2 = l2
        self.random_state = check_random_state(random_state)
        self.tol = tol
        self.max_iter = max_iter
        self.learning_rate = learning_rate

        # Attributes
        self.weights_ = None
        self.bias_ = None
        self.label_encoder = None
        self.n_features = None

    def _theano_init(self):
        n_features = self.n_features
        self.weights_ = theano.shared(
            value=self.random_state.rand(n_features).astype(self.floatX),
            name="weights"
        )

        self.bias_ = theano.shared(
            value=self.random_state.rand(1).astype(self.floatX),
            name="bias",
            broadcastable=(True,)
        )

        # X is of shape (n_samples, n_features)
        X = T.matrix()
        y = T.vector()

        n_samples, _ = X.shape

        # s is of shape (n_samples,)
        scores = self._scorer(X)

        regularization = self.l1 * self.weights_.norm(1) \
                         + 0.5 * self.l2 * self.weights_.norm(2) ** 2
        loss = -1 / n_samples * (y * scores - T.log(1 + T.exp(scores))).sum()
        loss += regularization

        # Gradient descent
        grad_weights = T.grad(cost=loss, wrt=self.weights_)
        grad_bias = T.grad(cost=loss, wrt=self.bias_)

        lr = self.learning_rate
        updates = [(self.weights_, self.weights_ -  lr * grad_weights),
                   (self.bias_, self.bias_ - lr * grad_bias)]

        predict = T.where(scores >= 0, 1, 0)

        # Storing the functions
        self._predict = theano.function([X], predict)
        self._grad_weights = theano.function([X, y], grad_weights)
        self._grad_bias = theano.function([X, y], grad_bias)

        self._gradient_step = theano.function(
            inputs=[X, y],
            outputs=loss,
            updates=updates
        )

    def _scorer(self, X):
        """Should work with either X being a theano variable or numpy variable.
        """
        return self.weights_.dot(X.T) + self.bias_

    def fit(self, X, y):
        X, y = check_X_y(X, y, dtype=self.floatX)

        self.n_features = X.shape[1]
        self._theano_init()

        if len(np.unique(y)) != 2:
            raise ValueError()

        # Preprocessing: make sure input is -1 or +1
        self.label_encoder = LabelEncoder()
        y = self.label_encoder.fit_transform(y)
        y[y == 0] = -1

        # Performs the gradient descents
        self._fit(X, y)
        return self

    def _fit(self, X, y):

        for _ in xrange(self.max_iter):
            g_weights = self._grad_weights(X, y)
            g_bias = self._grad_bias(X, y)
            grad_concat = np.hstack([g_weights, g_bias])

            loss = self._gradient_step(X, y)
            mag = np.linalg.norm(grad_concat, ord=1)
            if mag < self.tol * self.n_features:
                break

    def predict(self, X):
        X = check_array(X, dtype=self.floatX)
        if X.shape[1] != self.n_features:
            raise ValueError()

        # Binarize, and relabel the output
        predictions = self._predict(X)
        return self.label_encoder.inverse_transform(predictions)

sample.py

from sklearn.datasets import make_classification
from sklearn.metrics import accuracy_score
from linear_model import LogisticRegression
import sklearn

X, y = make_classifcation(n_samples=200, n_informative=20)
mylogreg = LogisticRegression(l1=0, l2=1).fit(X, y)
predictions = mylogreg.predict(X)
print accuracy_score(y, predictions)

# Now that you mention it, maybe I should make C = 1 / n_samples ?
sklearn_logreg = sklearn.linear_model.LogisticRegression(C=1.0)
sklearn_logreg.fit(X, y)
predictions = sklearn_logreg.predict(X)
print accuracy_score(y, predictions)