Feature Selection: Letting the Label Choose, Not Just Variance

Lecture 13, and the topic closes a loop opened back in the PCA post: reducing dimensionality without using the label can throw away important information. Feature selection fixes that in the most direct way possible, using the label to decide what to drop.

The dataset: digits, pixel by pixel

MNIST again, 60 thousand training images, 784 pixels each. The professor prints a digit as ASCII art (X where there's ink, blank where there isn't), a pretty literal reminder that each of these 784 "variables" is just one pixel, and most of them (the image borders, for instance) probably never have any ink at all.

class Normalizer255(BaseEstimator, TransformerMixin):
    def fit(self, X, y=None):
        return self
    def transform(self, X, y=None):
        return X / 255.0

pipeline = make_pipeline(Normalizer255(), RidgeClassifier())

Output (all 784 pixels): 0.8604 test accuracy.

Notice the normalization: instead of StandardScaler (which needs to compute the training set's mean and spread), the professor uses a fixed constant, 255, because he already knows ahead of time that a grayscale image pixel ranges from 0 to 255. Domain knowledge replacing a computation that would otherwise be unnecessary.

Dropping by variance: doesn't always help

from sklearn.feature_selection import VarianceThreshold
X_transformed = VarianceThreshold(0.1).fit_transform(X_train)

Output: from 784 pixels down to 695 (89 dropped, probably the borders that almost never have ink). Accuracy after the cut: 0.8449, worse than the original 784.

Worth stating that honestly: even a cut that looks obvious (drop a pixel that almost never changes) made the result worse here, instead of simplifying for free. VarianceThreshold only looks at X, never at the label y, so it has no way of knowing whether that small bit of variance it dropped was carrying signal relevant to telling the digits apart.

The homemade ruler: within-class versus between-class

The professor builds a vectorized Euclidean distance function from scratch (the same idea you can derive from $\|a-b\|^2 = \|a\|^2 - 2a\cdot b + \|b\|^2$, avoiding a loop over every pair of points) and uses it to define a quality ruler per variable:

def distance_score(X, y):
    in_dist, out_dist = 0, 0
    for label in np.unique(y):
        in_class = X[y == label]
        out_class = X[y != label]
        in_dist += pairwise_distances(in_class, in_class).mean()
        out_dist += pairwise_distances(in_class, out_class).mean()
    return in_dist / (out_dist + 1e-8)

The idea: for every variable, measure the average distance within each class (do all the "3"s look alike on this variable?) and the average distance between different classes (are a "3" and a "7" pretty different on this variable?). A good variable has low within-class distance and high between-class distance, so the within/between ratio comes out small. That's essentially the same idea as the Fisher discriminant Bishop describes (chapter 4.1.4, which I already cited back in the classification post in a different form): separation between classes divided by variation within each class, just applied variable by variable here, instead of along a learned projection direction.

The cost of computing this over everything

Computing pairwise distance across 60 thousand points is expensive (the cost grows with the square of the point count), so the professor cuts down to a sample of 1200 training examples before ranking:

Output: ranking the 784 variables on that smaller sample took around 20 seconds. Still much faster than trying this on all 60 thousand points.

A direct reminder that "how much this costs to compute" is also part of a technique's design, not just "what it measures."

Actually selecting

class UnivariatedRanking(BaseEstimator, TransformerMixin):
    def __init__(self, n_features):
        self.n_features = n_features
    def fit(self, X, y=None):
        self.scores_ = univariate_ranking(X, y)
        return self
    def transform(self, X, y=None):
        X_sorted = X[:, self.scores_.argsort()]
        return X_sorted[:, :self.n_features]

With only 1200 training examples (much less than the 60 thousand from before, to fit the time budget), the baseline with all 784 variables already drops quite a bit, to 0.7678, the price of training on less data. Selecting only the 180 best variables (23% of the total):

Output: 0.7681. Practically tied with using everything, with less than a quarter of the variables.

That's the core difference from PCA: here the label participates in the choice from the start, so the technique doesn't risk throwing away exactly the information that separates the classes, the same risk that dropped accuracy from 100% to 47% on that artificial dataset in the PCA post. Supervised variable selection and unsupervised dimensionality reduction solve similar-looking problems, but with very different guarantees.

Automating "how many variables to keep"

There was still the question of picking 180 somehow, and the professor solves that too: HybridRanking uses the univariate ranking to order the variables, then tests incrementally (1 variable, 2 variables, 3...) on a separate validation set, keeping whichever size gave the best accuracy:

for i in range(1, len(score_idxs)):
    X_selected = X_tr[:, score_idxs[0:i]]
    self.estimator.fit(X_selected, y_tr)
    acc = self.estimator.score(X_val[:, score_idxs[0:i]], y_val)
    if acc > best_score:
        best_score, best_set = acc, i

Output: the algorithm picked 199 variables on its own, with accuracy 0.7702551020408164, exactly matching, down to the last decimal, the result of manually choosing the ranking's top 199 variables.

The combination (a cheap ranking first, then testing sizes incrementally) is a classic feature-selection pattern: a cheap filter (distance_score, computed once, sorted) shrinks the search space, and a more expensive method (actually training and validating) decides the cutoff, without needing to train a model for every one of the $2^{784}$ possible variable combinations.

Wrapping up

Practical application

I reproduced the same ruler (distance_score, univariate ranking) on a smaller digits dataset, the same Digits (1797 images, 8×8 pixels) that already showed up in the K-means post, fast enough to run on all 64 pixels without needing to shrink the sample.

scores = univariate_ranking(X_train, y_train)
sorted_idxs = np.argsort(scores)
for n in [10, 20, 32]:
    model = RidgeClassifier().fit(X_train[:, sorted_idxs[:n]], y_train)
    acc = accuracy_score(y_test, model.predict(X_test[:, sorted_idxs[:n]]))

With half the variables (32 of 64), accuracy drops by just one percentage point (0.9389 to 0.9278). The curve isn't linear: cutting from 64 to 32 barely hurts, but cutting from 32 to 20 and then to 10 starts hurting fast, a sign that most of the useful signal really is concentrated in a fraction of the variables, exactly the assumption that makes feature selection a technique worth using.