Together with other scikit-learn developers we've created an umbrella organization for scikit-learn-related projects named scikit-learn-contrib. The idea is for this organization to host projects that are deemed too specific or too experimental to be included in the scikit-learn codebase but still offer an API which is compatible with scikit-learn and would like to benefit of the visibility of being labeled as scikit-learn-compatible.

We've set two requirements to being under this umbrella right now (this might evolve in the future). The first requirement is to have a scikit-learn compatible API, i.e., to follow the guide on the scikit-learn documentation so that objects can be used by scikit-learn meta-estimators (such as GridSearchCV). The second condition is that projects should be actively maintaned and have a high-quality codebase. Judging the quality of a codebase is difficult and subjective, but we agreed that at the bare minimum, the source code should be tested using continuous integration tools such as travis and reach a good test coverage (above 80%). More information is not available on the scikit-learn-contrib repository.

The first project to be hosted by this organization is lightning, but we hope that others will follow. If you would like to submit a new project, open an issue at the main project and we will look into it. There is also a project template for new and old projects.

Recently I've implemented, together with Arnaud Rachez, the SAGA[^{1}] algorithm in the lightning machine learning library (which by the way, has been recently moved to the new scikit-learn-contrib project). The lightning library uses the same API as scikit-learn but is particularly adapted to online learning. As for the SAGA algorithm, its performance is similar to other variance-reduced stochastic algorithms such as SAG[^{3}] or SVRG[^{2}] but it has the advantage with respect to SAG[^{3}] that it allows non-smooth penalty terms (such as $\ell_1$ regularization). It is implemented in lightning as SAGAClassifier and SAGARegressor.

We have taken care to make this implementation as efficient as possible. As for most stochastic gradient algorithms, a naive implementation takes 3 lines of code and is straightforward to implement. However, there are many tricks that are time-consuming and error-prone to implement but make a huge difference in efficiency.

A small example, more as a sanity check than to claim anything. The following plot shows the suboptimality as a function of time for three similar methods: SAG, SAGA and SVRG. The dataset used in the RCV1 dataset (test set, obtained from the libsvm webpage), consisting of 677.399 samples and 47.236 features. Interestingly, all methods can solve this rather large-scale problem within a few seconds. Within them, SAG and SAGA have a very similar performance and SVRG seems to be reasonably faster.

A note about the benchmarks: it is difficult to compare fairly stochastic gradient methods because at the end it usually boils down to how you choose the step size. In this plot I set the step size of all methods to 1/(3L) , where L is the Lipschitz constant of the objective function, as I think this is a popular choice. I would have prefered 1/L but SVRG was not converging for this step size. The code for the benchmarks can be found here.
Cross-validation iterators in scikit-learn are simply generator objects, that is, Python objects that implement the `__iter__`

method and that for each call to this method return (or more precisely, `yield`

) the indices or a boolean mask for the train and test set. Hence, implementing new cross-validation iterators that behave as the ones in scikit-learn is easy with this in mind. Here goes a small code snippet that implements a holdout cross-validator generator following the scikit-learn API.

import numpy as np
from sklearn.utils import check_random_state
class HoldOut:
"""
Hold-out cross-validator generator. In the hold-out, the
data is split only once into a train set and a test set.
Unlike in other cross-validation schemes, the hold-out
consists of only one iteration.
Parameters
----------
n : total number of samples
test_size : 0 < float < 1
Fraction of samples to use as test set. Must be a
number between 0 and 1.
random_state : int
Seed for the random number generator.
"""
def __init__(self, n, test_size=0.2, random_state=0):
self.n = n
self.test_size = test_size
self.random_state = random_state
def __iter__(self):
n_test = int(np.ceil(self.test_size * self.n))
n_train = self.n - n_test
rng = check_random_state(self.random_state)
permutation = rng.permutation(self.n)
ind_test = permutation[:n_test]
ind_train = permutation[n_test:n_test + n_train]
yield ind_train, ind_test

Contrary to other cross-validation schemes, holdout relies on a single split of the data. It is well known than in practice holdout performs much worse than KFold or LeaveOneOut schemes. However, holdout has the advantage that its theoretical properties are easier to derive. For examples of this see e.g. Section 8.7 of Theory of classification: a survey of some recent advances and the very recent The reusable holdout.

My latest contribution for scikit-learn is
an implementation of the isotonic regression model that I coded with
Nelle Varoquaux and
Alexandre Gramfort. This model
finds the best least squares fit to a set of points, given the
constraint that the fit must be a non-decreasing
function. The example
on the scikit-learn website gives an intuition on this model.

The original points are in red, and the estimated ones are in
green. As you can see, there is one estimation (green point) for each
data sample (red point). Calling $y \in \mathbb{R}^n$ the input data,
the model can be written concisely as an optimization problem over $x$

$$
\text{argmin}_x |y - x |^2 \\
\text{subject to } x_0 \leq x_1 \leq \cdots \leq x_n
$$

The algorithm implemented in scikit-learn ^{3} is the pool adjacent
violators algorithm ^{1}, which is an efficient linear time
$\mathcal{O}(n)$ algorithm. The algorithm sweeps through the data
looking for violations of the monotonicity constraint. When it finds
one, it adjusts the estimate to the best possible fit with
constraints. Sometimes it also needs to modify previous points to make
sure the new estimate does not violate the constraints. The following
picture shows how it proceeds at each iteration

This tutorial introduces the concept of pairwise preference used in most ranking problems. I'll use scikit-learn and for learning and matplotlib for visualization.

In the ranking setting, training data consists of lists of items with some order specified between items in each list. This order is typically induced by giving a numerical or ordinal score or a binary judgment (e.g. "relevant" or "not relevant") for each item, so that for any two samples `a`

and `b`

, either `a < b`

, `b > a`

or `b`

and `a`

are not comparable.

For example, in the case of a search engine, our dataset consists of results that belong to different queries and we would like to only compare the relevance for results coming from the same query.

This order relation is usually domain-specific. For instance, in information retrieval the set of comparable samples is referred to as a "query id". The goal behind this is to compare only documents that belong to the same query (Joachims 2002). In medical imaging on the other hand, the order of the labels usually depend on the subject so the comparable samples is given by the different subjects in the study (Pedregosa et al 2012).

import itertools
import numpy as np
from scipy import stats
import pylab as pl
from sklearn import svm, linear_model, cross_validation

To start with, we'll create a dataset in which the target values consists of three graded measurements Y = {0, 1, 2} and the input data is a collection of 30 samples, each one with two features.

The set of comparable elements (queries in information retrieval) will consist of two equally sized blocks, $X = X_1 \cup X_2$, where each block is generated using a normal distribution with different mean and covariance. In the pictures, we represent $X_1$ with round markers and $X_2$ with triangular markers.

np.random.seed(0)
theta = np.deg2rad(60)
w = np.array([np.sin(theta), np.cos(theta)])
K = 20
X = np.random.randn(K, 2)
y = [0] * K
for i in range(1, 3):
X = np.concatenate((X, np.random.randn(K, 2) + i * 4 * w))
y = np.concatenate((y, [i] * K))
# slightly displace data corresponding to our second partition
X[::2] -= np.array([3, 7])
blocks = np.array([0, 1] * (X.shape[0] / 2))
# split into train and test set
cv = cross_validation.StratifiedShuffleSplit(y, test_size=.5)
train, test = iter(cv).next()
X_train, y_train, b_train = X[train], y[train], blocks[train]
X_test, y_test, b_test = X[test], y[test], blocks[test]
# plot the result
idx = (b_train == 0)
pl.scatter(X_train[idx, 0], X_train[idx, 1], c=y_train[idx],
marker='^', cmap=pl.cm.Blues, s=100)
pl.scatter(X_train[~idx, 0], X_train[~idx, 1], c=y_train[~idx],
marker='o', cmap=pl.cm.Blues, s=100)
pl.arrow(0, 0, 8 * w[0], 8 * w[1], fc='gray', ec='gray',
head_width=0.5, head_length=0.5)
pl.text(0, 1, '$w$', fontsize=20)
pl.arrow(-3, -8, 8 * w[0], 8 * w[1], fc='gray', ec='gray',
head_width=0.5, head_length=0.5)
pl.text(-2.6, -7, '$w$', fontsize=20)
pl.axis('equal')
pl.show()

In the plot we clearly see that for both blocks there's a common vector w such that the projection onto w gives a list with the correct ordering.

However, because linear considers that output labels live in a metric space it will consider that all pairs are comparable. Thus if we fit this model to the problem above it will fit both blocks at the same time, yielding a result that is clearly not optimal. In the following plot we estimate $\hat{w}$ using an l2-regularized linear model.

ridge = linear_model.Ridge(1.)
ridge.fit(X_train, y_train)
coef = ridge.coef_ / linalg.norm(ridge.coef_)
pl.scatter(X_train[idx, 0], X_train[idx, 1], c=y_train[idx],
marker='^', cmap=pl.cm.Blues, s=100)
pl.scatter(X_train[~idx, 0], X_train[~idx, 1], c=y_train[~idx],
marker='o', cmap=pl.cm.Blues, s=100)
pl.arrow(0, 0, 7 * coef[0], 7 * coef[1], fc='gray', ec='gray',
head_width=0.5, head_length=0.5)
pl.text(2, 0, '$\hat{w}$', fontsize=20)
pl.axis('equal')
pl.title('Estimation by Ridge regression')
pl.show()

To assess the quality of our model we need to define a ranking score. Since we are interesting in a model that *orders* the data, it is natural to look for a metric that compares the ordering of our model to the given ordering. For this, we use Kendall's tau correlation coefficient, which is defined as (P - Q)/(P + Q), being P the number of concordant pairs and Q is the number of discordant pairs. This measure is used extensively in the ranking literature (e.g Optimizing Search Engines using Clickthrough Data).

We thus evaluate this metric on the test set for each block separately.

for i in range(2):
tau, _ = stats.kendalltau(
ridge.predict(X_test[b_test == i]), y_test[b_test == i])
print('Kendall correlation coefficient for block %s: %.5f' % (i, tau))

Kendall correlation coefficient for block 0: 0.71122
Kendall correlation coefficient for block 1: 0.84387

## The pairwise transform

As proved in (Herbrich 1999), if we consider linear ranking functions, the ranking problem can be transformed into a two-class classification problem. For this, we form the difference of all comparable elements such that our data is transformed into $(x'_k, y'_k) = (x_i - x_j, sign(y_i - y_j))$ for all comparable pairs.

This way we transformed our ranking problem into a two-class classification problem. The following plot shows this transformed dataset, and color reflects the difference in labels, and our task is to separate positive samples from negative ones. The hyperplane {x^T w = 0} separates these two classes.

# form all pairwise combinations
comb = itertools.combinations(range(X_train.shape[0]), 2)
k = 0
Xp, yp, diff = [], [], []
for (i, j) in comb:
if y_train[i] == y_train[j] \
or blocks[train][i] != blocks[train][j]:
# skip if same target or different group
continue
Xp.append(X_train[i] - X_train[j])
diff.append(y_train[i] - y_train[j])
yp.append(np.sign(diff[-1]))
# output balanced classes
if yp[-1] != (-1) ** k:
yp[-1] *= -1
Xp[-1] *= -1
diff[-1] *= -1
k += 1
Xp, yp, diff = map(np.asanyarray, (Xp, yp, diff))
pl.scatter(Xp[:, 0], Xp[:, 1], c=diff, s=60, marker='o', cmap=pl.cm.Blues)
x_space = np.linspace(-10, 10)
pl.plot(x_space * w[1], - x_space * w[0], color='gray')
pl.text(3, -4, '$\{x^T w = 0\}$', fontsize=17)
pl.axis('equal')
pl.show()

As we see in the previous plot, this classification is separable. This will not always be the case, however, in our training set there are no order inversions, thus the respective classification problem is separable.

We will now finally train an Support Vector Machine model on the transformed data.
This model is known as RankSVM, although we note that the pairwise transform is more general and can be used together with any linear model. We will then plot the training data together with the estimated coefficient $\hat{w}$ by RankSVM.

clf = svm.SVC(kernel='linear', C=.1)
clf.fit(Xp, yp)
coef = clf.coef_.ravel() / linalg.norm(clf.coef_)
pl.scatter(X_train[idx, 0], X_train[idx, 1], c=y_train[idx],
marker='^', cmap=pl.cm.Blues, s=100)
pl.scatter(X_train[~idx, 0], X_train[~idx, 1], c=y_train[~idx],
marker='o', cmap=pl.cm.Blues, s=100)
pl.arrow(0, 0, 7 * coef[0], 7 * coef[1], fc='gray', ec='gray',
head_width=0.5, head_length=0.5)
pl.arrow(-3, -8, 7 * coef[0], 7 * coef[1], fc='gray', ec='gray',
head_width=0.5, head_length=0.5)
pl.text(1, .7, '$\hat{w}$', fontsize=20)
pl.text(-2.6, -7, '$\hat{w}$', fontsize=20)
pl.axis('equal')
pl.show()

Finally we will check that as expected, the ranking score (Kendall tau) increases with the RankSVM model respect to linear regression.

for i in range(2):
tau, _ = stats.kendalltau(
np.dot(X_test[b_test == i], coef), y_test[b_test == i])
print('Kendall correlation coefficient for block %s: %.5f' % (i, tau))

Kendall correlation coefficient for block 0: 0.83627
Kendall correlation coefficient for block 1: 0.84387

This is indeed higher than the values (0.71122, 0.84387) obtained in the case of linear regression.

**Original ipython notebook for this blog post can be found here**

Ridge coefficients for multiple values of the regularization parameter
can be elegantly computed by updating the *thin* SVD decomposition of
the design matrix:

import numpy as np
from scipy import linalg
def ridge(A, b, alphas):
"""
Return coefficients for regularized least squares
min ||A x - b||^2 + alpha ||x||^2
Parameters
----------
A : array, shape (n, p)
b : array, shape (n,)
alphas : array, shape (k,)
Returns
-------
coef: array, shape (p, k)
"""
U, s, Vt = linalg.svd(X, full_matrices=False)
d = s / (s[:, np.newaxis].T ** 2 + alphas[:, np.newaxis])
return np.dot(d * U.T.dot(y), Vt).T

This can be used to efficiently compute what it's *regularization
path*, that is, to plot the coefficients as a function of the
regularization parameter. Since the bottleneck of the algorithm is the
singular value decomposition, computing the coefficients for other
values of the regularization parameter basically comes for free.

A variant of this algorithm can then be used to compute the
optimal regularization parameter in the sense of leave-one-out
cross-validation and is implemented in scikit-learn's RidgeCV (for
which Mathieu Blondel has an excelent post by ). This optimal
parameter is denoted with a vertical dotted line in the following
picture, full code can be found here.

I've been working lately with Alexandre Gramfort coding the LARS
algorithm in scikits.learn. This algorithm computes the solution to
several general linear models used in machine learning: LAR, Lasso,
Elasticnet and Forward Stagewise. Unlike the implementation by
coordinate descent, the LARS algorithm gives the full coefficient path
along the regularization parameter, and thus it is specially well suited
for performing model selection.

The algorithm is coded
mostly in python, with some tiny parts in C (because I already had the
code for Cholesky deletes in C) and a cython interface for the BLAS
function dtrsv, which will be proposed to scipy once I stabilize this
code. The algorithm is mostly complete, allowing some optimizations,
like using a precomputed Gram matrix or specify maximum number of
features/iterations, but could still be extended to compute other
models, like ElasticNet or Forward Stagewise. I haven't done any
benchmarks yet, but preliminary ones by Alexandre Gramfort showed that
it is roughly equivalent to this Matlab implementation. Using
PyMVPA, it shouldn't be difficult to benchmark it against the R
implementation, though.