Recursive Feature Elimination
Recursive feature elimination (RFE) repeatedly fits a model, removes its lowest-ranked inputs, and fits again. The repeated fitting matters because a feature's apparent importance can change when another feature disappears.
Unlike an individual feature score, RFE judges features through a model trained with the current subset. This makes its choices depend on the learner, preprocessing, data and removal schedule. It is a search procedure, not a certificate that the surviving columns are the truly relevant ones.
Follow the fit–rank–remove loop
Choose the estimator, how to read its importance, the number of features to retain, and how many to remove per round. With one-at-a-time removal from $P$ features down to $K$, the procedure fits $P-K+1$ models, including the final fit.
For a single-target linear model on standardized inputs, absolute coefficient magnitude supplies a simple ranking. A negative coefficient is not automatically weak: its magnitude can be large. Scikit-learn squares these coefficients for ranking, which gives the same order in this single-target example.
Refitting can reverse the next removal decision
Use a small synthetic regression task with four inputs. A and B generate the target signal; C is a noisy copy of A; D is independent noise. There are only 30 training rows, so chance relationships can affect the fit. The reproducible generator is below.
Standardize using training statistics and fit ridge regression with penalty parameter 1. We will keep two columns. The bars show the absolute coefficients after each actual fit; a removed feature is marked as absent rather than assigned a new coefficient.
Keeping the two largest coefficients from the initial fit would choose A and B. RFE chooses A and C because it recomputes coefficients after removing D. The comparison demonstrates a different procedure; it does not establish that the recursive choice is better.
Correlated inputs can share or trade influence across fits. A low coefficient in one fitted subset does not prove that an input is intrinsically useless. Scaling also matters: coefficient magnitudes measured against different input units are not directly comparable. See standardization for the training-time transformation.
Check whether the selected model predicts well
Our generator uses B in the target, yet this run removes it. Evaluate on 4,000 fresh synthetic rows from the same generating rule. Mean squared error averages squared prediction mistakes; lower is better.
The lesson is to validate the complete selection procedure. Increasing sophistication does not remove sampling noise or make a greedy path globally optimal. Compare against the full eligible input set and sensible simpler baselines.
Python: reproduce the coefficient trace, RFE outputs and evaluation
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Ridge
from sklearn.feature_selection import RFE
from sklearn.metrics import mean_squared_error
def make_data(seed, n):
rng = np.random.default_rng(seed)
signal = rng.normal(size=(n, 2))
X = np.column_stack([
signal[:, 0], signal[:, 1],
signal[:, 0] + .3 * rng.normal(size=n),
rng.normal(size=n),
])
y = 2 * signal[:, 0] + .5 * signal[:, 1] + .8 * rng.normal(size=n)
return X, y
X, y = make_data(1, 30)
scaler = StandardScaler().fit(X) # Training rows only.
Z = scaler.transform(X)
names = np.array(["A", "B", "C", "D"])
keep, trace = list(range(4)), []
while True:
model = Ridge(alpha=1).fit(Z[:, keep], y)
trace.append((keep.copy(), model.coef_.copy()))
print(names[keep], np.round(model.coef_, 4))
if len(keep) == 2:
break
weakest = int(np.argmin(np.abs(model.coef_)))
keep.pop(weakest) # Remove a local position; retain original column indices.
selector = RFE(Ridge(alpha=1), n_features_to_select=2, step=1).fit(Z, y)
print(selector.support_) # [True, False, True, False]
print(selector.ranking_) # [1, 2, 1, 3]
test_X, test_y = make_data(100, 4000) # Separate synthetic evaluation sample.
test_Z = scaler.transform(test_X)
initial_top_two = sorted(np.argsort(np.abs(trace[0][1]))[-2:])
for label, columns in [("All four", list(range(4))),
("Initial top two", initial_top_two), ("RFE", keep)]:
fitted = Ridge(alpha=1).fit(Z[:, columns], y)
mse = mean_squared_error(test_y, fitted.predict(test_Z[:, columns]))
print(label, mse)
# All four: 0.874248; initial top two (A,B): 0.682276; RFE (A,C): 1.001822.
The generator uses independent standard-normal sources. The target is twice A plus half B plus independent noise with standard deviation 0.8; C adds independent noise with standard deviation 0.3 to A. The fitted scaler is reused for the evaluation sample.
Read the outputs without inventing extra meaning
support_ identifies the retained original columns. Here it is true for A and C. ranking_ is 1 for every retained feature, 2 for B, removed last, and 3 for D, removed first.
These ranks describe the elimination history. They do not say A and C are equally predictive or that rank 2 has half the value of rank 1. When several features are removed together, they can share a rank. Save the support mask with the preprocessing and final estimator so prediction inputs remain aligned.
Choose the removal schedule and feature count
A larger removal step reduces the number of model fits, but it skips intermediate subsets. Several features can be removed using scores from the same fit, before the remaining model has a chance to adjust. Ties and estimator randomness can also affect the path; set reproducible seeds where needed and inspect stability across training folds.
The standard estimator must expose a usable importance, such as coefficients or tree importances. A fitted pipeline can expose it through an importance_getter path. The selector inherits the weaknesses of that importance measure; tree feature importance discusses its common biases.
RFE uses a feature count you choose. RFECV uses cross-validation to choose that count. RFECV evaluates elimination paths inside its training folds, then selects a count and refits using its full input training set. Different folds can retain different feature identities even at the same count.
Keep preprocessing and selection inside the evaluation boundary
For fixed-count RFE, place training-fitted preprocessing and selection inside the pipeline passed to cross-validation. If you choose the count with RFECV, its inner folds also need their own preprocessing fit. Putting a scaler before RFECV in an outer pipeline fits that scaler before RFECV makes its inner splits.
Python: give each RFECV fit its own scaler
from sklearn.feature_selection import RFECV
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import KFold
estimator = make_pipeline(StandardScaler(), Ridge(alpha=1))
selector_cv = RFECV(
estimator,
importance_getter="named_steps.ridge.coef_",
cv=KFold(5, shuffle=True, random_state=12),
scoring="neg_mean_squared_error",
min_features_to_select=1,
step=1,
).fit(X, y) # Only the development training data.
Use an untouched evaluation set or an outer cross-validation loop to assess this count-selection procedure. The maximum inner-CV score was used to choose a count, so it is not an independent final evaluation.
Preserve group or time boundaries when observations require them. If you use a different downstream learner from the one driving elimination, validate that combination too. The nested-CV chapter explains how to evaluate a procedure that tunes its choices internally.
Implement actual refitting in the elimination loop
The C++ version below fits ridge coefficients after every removal, using coordinate updates. Its input contract is a dense training matrix already standardized with training statistics and a target centered by its training mean. Call eliminate(Z, centered_y, 2) to retain two columns.
C++17: ridge fitting and one-at-a-time elimination
#include <algorithm>
#include <cmath>
#include <numeric>
#include <stdexcept>
#include <vector>
using Matrix = std::vector<std::vector<double>>;
struct Selection { std::vector<std::size_t> columns; std::vector<double> weights; };
// Internal ridge solver: Z columns and y are already centered.
std::vector<double> ridge_fit(const Matrix& Z, const std::vector<double>& y,
const std::vector<std::size_t>& columns, double alpha) {
std::vector<double> weights(columns.size(), 0), residual = y;
for (int sweep = 0; sweep < 10000; ++sweep) {
double change = 0;
for (std::size_t j = 0; j < columns.size(); ++j) {
double numerator = 0, denominator = alpha;
for (std::size_t i = 0; i < Z.size(); ++i) {
double x = Z[i][columns[j]];
numerator += x * (residual[i] + x * weights[j]);
denominator += x * x;
}
double next = numerator / denominator, delta = next - weights[j];
if (!std::isfinite(next)) throw std::runtime_error("nonfinite ridge update");
for (std::size_t i = 0; i < Z.size(); ++i)
residual[i] -= Z[i][columns[j]] * delta;
change = std::max(change, std::abs(delta)); weights[j] = next;
}
if (change < 1e-10) return weights;
}
throw std::runtime_error("ridge solver did not converge");
}
Selection eliminate(const Matrix& Z, const std::vector<double>& centered_y,
std::size_t target, double alpha = 1) {
if (Z.empty() || Z[0].empty() || centered_y.size() != Z.size())
throw std::invalid_argument("invalid training shape");
const auto P = Z[0].size();
if (target == 0 || target > P || !std::isfinite(alpha) || alpha <= 0)
throw std::invalid_argument("invalid target count or ridge penalty");
for (const auto& row : Z) {
if (row.size() != P) throw std::invalid_argument("ragged matrix");
for (double x : row) if (!std::isfinite(x))
throw std::invalid_argument("finite features required");
}
for (double y : centered_y) if (!std::isfinite(y))
throw std::invalid_argument("finite labels required");
std::vector<std::size_t> keep(P); std::iota(keep.begin(), keep.end(), 0);
while (true) {
auto weights = ridge_fit(Z, centered_y, keep, alpha); // Refit after every removal.
if (keep.size() == target) return {keep, weights};
auto weakest = std::min_element(weights.begin(), weights.end(),
[](double a, double b) { return std::abs(a) < std::abs(b); });
keep.erase(keep.begin() + (weakest - weights.begin()));
}
}
The ridge helper minimizes squared residuals plus alpha × sum of squared weights, matching the Python ridge objective. It holds the other weights fixed while updating one coefficient and adjusts the residual incrementally. The public function validates shapes and finite inputs; the internal solver assumes that validation. It stops on small coefficient changes and reports failure if it does not converge.
The result contains original column indices and coefficients from the final refit. For prediction, apply the saved feature standardization, select those columns, and add back the saved target mean. This teaching implementation removes one column per round and uses the first minimum to resolve an exact tie.
References and further reading
- Scikit-learn: RFE — the fit–rank–remove procedure, importance access and output ranks.
- Scikit-learn: RFECV — cross-validated feature-count selection and fold-specific outputs.
- Scikit-learn: Ridge — the regression objective used for the numerical example.
- Guyon et al.: Gene Selection for Cancer Classification using Support Vector Machines — the original SVM-based RFE method; this chapter uses ridge regression for its small illustration.