175 lines
7.3 KiB
Python
175 lines
7.3 KiB
Python
"""
|
|
============================================================
|
|
Sparse recovery: feature selection for sparse linear models
|
|
============================================================
|
|
|
|
Given a small number of observations, we want to recover which features
|
|
of X are relevant to explain y. For this :ref:`sparse linear models
|
|
<l1_feature_selection>` can outperform standard statistical tests if the
|
|
true model is sparse, i.e. if a small fraction of the features are
|
|
relevant.
|
|
|
|
As detailed in :ref:`the compressive sensing notes
|
|
<compressive_sensing>`, the ability of L1-based approach to identify the
|
|
relevant variables depends on the sparsity of the ground truth, the
|
|
number of samples, the number of features, the conditioning of the
|
|
design matrix on the signal subspace, the amount of noise, and the
|
|
absolute value of the smallest non-zero coefficient [Wainwright2006]
|
|
(http://statistics.berkeley.edu/sites/default/files/tech-reports/709.pdf).
|
|
|
|
Here we keep all parameters constant and vary the conditioning of the
|
|
design matrix. For a well-conditioned design matrix (small mutual
|
|
incoherence) we are exactly in compressive sensing conditions (i.i.d
|
|
Gaussian sensing matrix), and L1-recovery with the Lasso performs very
|
|
well. For an ill-conditioned matrix (high mutual incoherence),
|
|
regressors are very correlated, and the Lasso randomly selects one.
|
|
However, randomized-Lasso can recover the ground truth well.
|
|
|
|
In each situation, we first vary the alpha parameter setting the sparsity
|
|
of the estimated model and look at the stability scores of the randomized
|
|
Lasso. This analysis, knowing the ground truth, shows an optimal regime
|
|
in which relevant features stand out from the irrelevant ones. If alpha
|
|
is chosen too small, non-relevant variables enter the model. On the
|
|
opposite, if alpha is selected too large, the Lasso is equivalent to
|
|
stepwise regression, and thus brings no advantage over a univariate
|
|
F-test.
|
|
|
|
In a second time, we set alpha and compare the performance of different
|
|
feature selection methods, using the area under curve (AUC) of the
|
|
precision-recall.
|
|
"""
|
|
print(__doc__)
|
|
|
|
# Author: Alexandre Gramfort and Gael Varoquaux
|
|
# License: BSD 3 clause
|
|
|
|
import warnings
|
|
|
|
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
from scipy import linalg
|
|
|
|
from sklearn.linear_model import (RandomizedLasso, lasso_stability_path,
|
|
LassoLarsCV)
|
|
from sklearn.feature_selection import f_regression
|
|
from sklearn.preprocessing import StandardScaler
|
|
from sklearn.metrics import auc, precision_recall_curve
|
|
from sklearn.ensemble import ExtraTreesRegressor
|
|
from sklearn.utils.extmath import pinvh
|
|
from sklearn.exceptions import ConvergenceWarning
|
|
|
|
|
|
def mutual_incoherence(X_relevant, X_irelevant):
|
|
"""Mutual incoherence, as defined by formula (26a) of [Wainwright2006].
|
|
"""
|
|
projector = np.dot(np.dot(X_irelevant.T, X_relevant),
|
|
pinvh(np.dot(X_relevant.T, X_relevant)))
|
|
return np.max(np.abs(projector).sum(axis=1))
|
|
|
|
|
|
for conditioning in (1, 1e-4):
|
|
###########################################################################
|
|
# Simulate regression data with a correlated design
|
|
n_features = 501
|
|
n_relevant_features = 3
|
|
noise_level = .2
|
|
coef_min = .2
|
|
# The Donoho-Tanner phase transition is around n_samples=25: below we
|
|
# will completely fail to recover in the well-conditioned case
|
|
n_samples = 25
|
|
block_size = n_relevant_features
|
|
|
|
rng = np.random.RandomState(42)
|
|
|
|
# The coefficients of our model
|
|
coef = np.zeros(n_features)
|
|
coef[:n_relevant_features] = coef_min + rng.rand(n_relevant_features)
|
|
|
|
# The correlation of our design: variables correlated by blocs of 3
|
|
corr = np.zeros((n_features, n_features))
|
|
for i in range(0, n_features, block_size):
|
|
corr[i:i + block_size, i:i + block_size] = 1 - conditioning
|
|
corr.flat[::n_features + 1] = 1
|
|
corr = linalg.cholesky(corr)
|
|
|
|
# Our design
|
|
X = rng.normal(size=(n_samples, n_features))
|
|
X = np.dot(X, corr)
|
|
# Keep [Wainwright2006] (26c) constant
|
|
X[:n_relevant_features] /= np.abs(
|
|
linalg.svdvals(X[:n_relevant_features])).max()
|
|
X = StandardScaler().fit_transform(X.copy())
|
|
|
|
# The output variable
|
|
y = np.dot(X, coef)
|
|
y /= np.std(y)
|
|
# We scale the added noise as a function of the average correlation
|
|
# between the design and the output variable
|
|
y += noise_level * rng.normal(size=n_samples)
|
|
mi = mutual_incoherence(X[:, :n_relevant_features],
|
|
X[:, n_relevant_features:])
|
|
|
|
###########################################################################
|
|
# Plot stability selection path, using a high eps for early stopping
|
|
# of the path, to save computation time
|
|
alpha_grid, scores_path = lasso_stability_path(X, y, random_state=42,
|
|
eps=0.05)
|
|
|
|
plt.figure()
|
|
# We plot the path as a function of alpha/alpha_max to the power 1/3: the
|
|
# power 1/3 scales the path less brutally than the log, and enables to
|
|
# see the progression along the path
|
|
hg = plt.plot(alpha_grid[1:] ** .333, scores_path[coef != 0].T[1:], 'r')
|
|
hb = plt.plot(alpha_grid[1:] ** .333, scores_path[coef == 0].T[1:], 'k')
|
|
ymin, ymax = plt.ylim()
|
|
plt.xlabel(r'$(\alpha / \alpha_{max})^{1/3}$')
|
|
plt.ylabel('Stability score: proportion of times selected')
|
|
plt.title('Stability Scores Path - Mutual incoherence: %.1f' % mi)
|
|
plt.axis('tight')
|
|
plt.legend((hg[0], hb[0]), ('relevant features', 'irrelevant features'),
|
|
loc='best')
|
|
|
|
###########################################################################
|
|
# Plot the estimated stability scores for a given alpha
|
|
|
|
# Use 6-fold cross-validation rather than the default 3-fold: it leads to
|
|
# a better choice of alpha:
|
|
# Stop the user warnings outputs- they are not necessary for the example
|
|
# as it is specifically set up to be challenging.
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter('ignore', UserWarning)
|
|
warnings.simplefilter('ignore', ConvergenceWarning)
|
|
lars_cv = LassoLarsCV(cv=6).fit(X, y)
|
|
|
|
# Run the RandomizedLasso: we use a paths going down to .1*alpha_max
|
|
# to avoid exploring the regime in which very noisy variables enter
|
|
# the model
|
|
alphas = np.linspace(lars_cv.alphas_[0], .1 * lars_cv.alphas_[0], 6)
|
|
clf = RandomizedLasso(alpha=alphas, random_state=42).fit(X, y)
|
|
trees = ExtraTreesRegressor(100).fit(X, y)
|
|
# Compare with F-score
|
|
F, _ = f_regression(X, y)
|
|
|
|
plt.figure()
|
|
for name, score in [('F-test', F),
|
|
('Stability selection', clf.scores_),
|
|
('Lasso coefs', np.abs(lars_cv.coef_)),
|
|
('Trees', trees.feature_importances_),
|
|
]:
|
|
precision, recall, thresholds = precision_recall_curve(coef != 0,
|
|
score)
|
|
plt.semilogy(np.maximum(score / np.max(score), 1e-4),
|
|
label="%s. AUC: %.3f" % (name, auc(recall, precision)))
|
|
|
|
plt.plot(np.where(coef != 0)[0], [2e-4] * n_relevant_features, 'mo',
|
|
label="Ground truth")
|
|
plt.xlabel("Features")
|
|
plt.ylabel("Score")
|
|
# Plot only the 100 first coefficients
|
|
plt.xlim(0, 100)
|
|
plt.legend(loc='best')
|
|
plt.title('Feature selection scores - Mutual incoherence: %.1f'
|
|
% mi)
|
|
|
|
plt.show()
|