sag_multinomial_benchmark.py

import gc
import time
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import fetch_rcv1
from sklearn.preprocessing import LabelBinarizer
from sklearn.linear_model import LogisticRegression
from sklearn.linear_model.logistic import _multinomial_loss
from sklearn.externals.joblib import Memory

m = Memory(cachedir='.', verbose=0)


# compute multinomial logistic loss
def get_loss(w, intercept, X, y, C):
    w = np.vstack((w.T, intercept.T)).T.ravel()
    lbin = LabelBinarizer()
    Y_multi = lbin.fit_transform(y)
    if Y_multi.shape[1] == 1:
        Y_multi = np.hstack([1 - Y_multi, Y_multi])
    loss, _, _ = _multinomial_loss(w, X, Y_multi, 1 / C, np.ones(X.shape[0]))
    return loss


# We use joblib to cache individual fits.
@m.cache(ignore=('X', 'y'))
def bench_one(X_shape, X, y, name, clf_type, clf_params, n_iter):
    clf = clf_type(**clf_params)
    clf.set_params(max_iter=n_iter, random_state=42)

    st = time.time()
    clf.fit(X, y)
    end = time.time()

    train_loss = get_loss(clf.coef_, clf.intercept_, X, y, clf.C)
    train_score = clf.score(X, y)
    test_score = clf.score(X_test, y_test)
    duration = end - st

    return train_loss, train_score, test_score, duration


def bench(clf, results):
    for (name, solver, iter_range, train_losses, train_scores,
         test_scores, durations) in results:
        print("training %s" % name)
        clf.set_params(solver=solver)
        clf_type = type(clf)
        clf_params = clf.get_params()

        for n_iter in iter_range:
            gc.collect()

            train_loss, train_score, test_score, duration = bench_one(
                X.shape, X, y, name, clf_type, clf_params, n_iter)

            train_losses.append(train_loss)
            train_scores.append(train_score)
            test_scores.append(test_score)
            durations.append(duration)
            print("classifier: %s" % name)
            print("train_loss: %.8f" % train_loss)
            print("train_score: %.8f" % train_score)
            print("test_score: %.8f" % test_score)
            print("time for fit: %.8f seconds" % duration)
            print("")

        print("")
    return results


def plot_train_losses(results):
    plt.figure()
    for (name, _, _, train_losses, _, _, durations) in results:
        plt.plot(durations, train_losses, '-o', label=name)
        plt.legend(loc=0)
        plt.xlabel("seconds")
        plt.ylabel("train loss")


def plot_train_scores(results):
    plt.figure()
    for (name, _, _, _, train_scores, _, durations) in results:
        plt.plot(durations, train_scores, '-o', label=name)
        plt.legend(loc=0)
        plt.xlabel("seconds")
        plt.ylabel("train score")
        plt.ylim((0.88, 1))


def plot_test_scores(results):
    plt.figure()
    for (name, _, _, _, _, test_scores, durations) in results:
        plt.plot(durations, test_scores, '-o', label=name)
        plt.legend(loc=0)
        plt.xlabel("seconds")
        plt.ylabel("test score")
        plt.ylim((0.88, 1))


def plot_dloss(results):
    plt.figure()
    pobj_final = []
    for (name, _, _, train_losses, _, _, durations) in results:
        pobj_final.append(train_losses[-1])

    indices = np.argsort(pobj_final)
    pobj_best = pobj_final[indices[0]]

    for (name, _, _, train_losses, _, _, durations) in results:
        log_pobj = np.log(abs(np.array(train_losses) - pobj_best)) / np.log(10)

        plt.plot(durations, log_pobj, '-o', label=name)
        plt.legend(loc=0)
        plt.xlabel("seconds")
        plt.ylabel("log(best - train_loss)")

# Load RCV1 dataset (large, sparse and multinomial)
rcv1 = fetch_rcv1()
X = rcv1.data
y = rcv1.target
# take only 4 categories (CCAT, ECAT, GCAT, MCAT)
y = y[:, [1, 4, 6, 10]].astype(np.float64)
# remove samples that have more than one category
mask = np.asarray(y.sum(axis=1) == 1).ravel()
y = y[mask, :].indices
X = X[mask, :]

# Split in training and testing sets
n = 23149
X_test = X[:n, :]
y_test = y[:n]
X = X[n:, :]
y = y[n:]

clf = LogisticRegression(C=1., tol=1e-15, multi_class='multinomial',
                         fit_intercept=True)

iter_range = list(range(1, 15, 1))
lbfgs_iter_range = list(range(1, 100, 10))

results = [("SAG", "sag", iter_range, [], [], [], []),
           ("Newton-CG", "newton-cg", iter_range, [], [], [], []),
           ("L-BFGS", "lbfgs", lbfgs_iter_range, [], [], [], [])]

results = bench(clf, results)

plot_train_scores(results)
plot_test_scores(results)
plot_train_losses(results)
plot_dloss(results)
plt.show()