如果特征N的数量很小，那么穷举搜索可能是可行的:比如说尝试所有可能的特征组合，只保留成本/目标函数最小的那一个。但是如果N很大，那么穷举搜索肯定是不可能的。因为对于N的组合是一个指数函数，所以在这种情况下，必须使用启发式方法:以一种有效的方式探索搜索空间，寻找能够最小化用于执行搜索的目标函数的特征组合。
找到一个好的启发式算法并非易事。R中的regsubsets()函数有一些可以使用的选项。此外，scikit-learn提供了几种执行启发式特征选择的方法，但是找到一个好的、通用的启发式——以最通用的形式——是一个难题。所以本文将探讨一些即使在N相当大的情况下也能很好地工作的方法。
数据集
我们这里使用Kaggle上非常流行的House Prices数据集(MIT许可)——然后经过一些简单的特征转换后，最终得到一个213个特征(N=213)和1453个观察值的数据集。我使用的模型是线性回归，statsmodels.api.OLS()，我们试图最小化的目标函数是BIC，贝叶斯信息标准，一种信息损失的度量，所以BIC越低越好。它与AIC相似，只不过BIC倾向于生成更简洁的模型(它更喜欢特征较少的模型)。最小化AIC或BIC可以减少过拟合。但也可以使用其他目标函数，例如r方(目标中已解释的方差)或调整后的r方——只要记住r平方越大越好，所以这是一个最大化问题。
目标函数的选择在这里是无关紧要的。重要的是我们要有一个目标函数，因为需要尝试使用各种技术来优化它。
但首先，让我们看看一些众所周知的、久经考验的特征选择技术，我们将使用它与后面描述的更复杂的技术进行基线比较。
顺序特征搜索 SFS
SFS相当简单，它从选择一个特征开始，然后选择使目标函数最小的特征。一旦一个特性被选中，它就会永远被选中。然后它会尝试添加另一个特征(这时就有2个了)，然后再次最小化目标。当所有特征一起被尝试时，搜索就结束了，使目标最小化的组合最为特征选择的结果。
SFS是一种贪婪算法——每个选择都是局部最优的——而且它永远不会回去纠正自己的错误。但他有一个优点就是即使N很大，它还是很快的。它尝试的组合总数是N(N+1)/2，这是一个二次多项式。
在Python中使用mlextend库的SFS代码是这样的
import statsmodels.api as sm
from mlxtend.feature_selection import SequentialFeatureSelector as SFS
from sklearn.base import BaseEstimator
class DummyEstimator(BaseEstimator):
# mlxtend wants to use an sklearn estimator, which is not needed here
# (statsmodels OLS is used instead)
# create a dummy estimator to pacify mlxtend
def fit(self, X, y=None, **kwargs):
return self
def neg_bic(m, X, y):
# objective function
lin_mod_res = sm.OLS(y, X, hasconst=True).fit()
return -lin_mod_res.bic
seq_selector = SFS(
DummyEstimator(),
k_features=(1, X.shape[1]),
forward=True,
floating=False,
scoring=neg_bic,
cv=None,
n_jobs=-1,
verbose=0,
# make sure the intercept is not dropped
fixed_features=['const'],
)
n_features = X.shape[1] - 1
objective_runs_sfs = round(n_features * (n_features + 1) / 2)
t_start_seq = time.time()
# mlxtend will mess with your dataframes if you don't .copy()
seq_res = seq_selector.fit(X.copy(), y.copy())
t_end_seq = time.time()
run_time_seq = t_end_seq - t_start_seq
seq_metrics = seq_res.get_metric_dict()
它快速找到了这些组合:
best k: 36
best objective: 33708.98602877906
R2 @ best k: 0.9075677543649224
objective runs: 22791
total run time: 44.850 sec
213个特征，最好的是36个。最好的BIC是33708.986，在我的系统上完成不到1分钟。它调用目标函数22.8k次。
以下是最佳BIC值和r平方值，作为所选特征数量的函数:
现在我们来试试更复杂的。
协方差矩阵自适应演化 CMA-ES
这是一个数值优化算法。它与遗传算法属于同一类(它们都是进化的)，但CMA-ES与遗传算法截然不同。它是一个随机算法，没有导数，不需要计算目标函数的导数(不像梯度下降，它依赖于偏导数)。它的计算效率很高，被用于各种数值优化库，如Optuna。在这里我们只简要介绍一下CMA-ES，有关更详细的解释，请参阅最后链接部分的文献。
考虑二维Rastrigin函数:
下面的热图显示了这个函数的值——颜色越亮意味着值越高。该函数在原点(0,0)处具有全局最小值，但它夹杂着许多局部极值。我们想通过CMA-ES找到全局最小值。
CMA-ES基于多元正态分布。它从这个分布中生成搜索空间中的测试点。但在开始之前，必须必须猜测分布的原始均值，以及它的标准差，但之后算法会迭代地修改所有这些参数，在搜索空间中扫描分布，寻找最佳的目标函数值。
Xi是算法在搜索空间中每一步产生的点的集合。Lambda是生成的点的个数。分布的平均值在每一步都会被更新，并且最终会收敛于真正的解。Sigma是分布的标准差——测试点的分布。C是协方差矩阵，它定义了分布的形状。根据C值的不同，分布可能呈“圆形”或更细长的椭圆形。对C的修改允许CMA-ES“潜入”搜索空间的某些区域，或避开其他区域。
上图中生成了6个点的种群，这是优化器为这个问题选择的默认种群大小。这是第一步。然后算法进行下面的步骤:
1、计算每个点的目标函数(Rastrigin)
2、更新均值、标准差和协方差矩阵，根据从目标函数中学到的信息，有效地创建一个新的多元正态分布
3、从新的分布中生成一组新的测试点
4、重复，直到满足某个准则(要么收敛于某个平均值，要么超过最大步数等)。
对于更新分布的均值是很简单的:计算每个测试点的目标函数后，给这些点分配权重，目标值越好的点权重越大，并从它们的位置计算加权和，这就成为新的均值。CMA-ES有效地将分布均值向目标值较好的点移动:
如果算法收敛于真解，那么分布的平均值也将收敛于该解。标准差会收敛于0。协方差矩阵将根据目标函数的位置改变分布的形状(圆形或椭圆形)，扩展到有希望的区域，并避开不好的区域。
下面是一个显示了CMA-ES解决拉斯特里金问题时测试点的时间演变的GIF动画:
将CMA-ES用于特征选择
2D Rastrigin函数相对简单，因为它只有2个维度。但对于我们的特征选择问题，有N=213个维度。而且空间是不连续的。每个测试点是一个n维向量，其分量值为{0,1}。换句话说，每个测试点看起来像这样:[1,0,0,1,1,1,0，…]-一个二进制向量。
下面是使用cmaes库进行特征选择的CMA-ES代码的简单版本:
def cma_objective(fs):
features_use = ['const'] + [
f for i, f in enumerate(features_select) if fs[i,] == 1
]
lin_mod = sm.OLS(y_cmaes, X_cmaes[features_use], hasconst=True).fit()
return lin_mod.bic
X_cmaes = X.copy()
y_cmaes = y.copy()
features_select = [f for f in X_cmaes.columns if f != 'const']
dim = len(features_select)
bounds = np.tile([0, 1], (dim, 1))
steps = np.ones(dim)
optimizer = CMAwM(
mean=np.full(dim, 0.5),
sigma=1 / 6,
bounds=bounds,
steps=steps,
n_max_resampling=10 * dim,
seed=0,
)
max_gen = 100
best_objective_cmaes_small = np.inf
best_sol_raw_cmaes_small = None
for gen in tqdm(range(max_gen)):
solutions = []
for _ in range(optimizer.population_size):
x_for_eval, x_for_tell = optimizer.ask()
value = cma_objective(x_for_eval)
solutions.append((x_for_tell, value))
if value < best_objective_cmaes_small:
best_objective_cmaes_small = value
best_sol_raw_cmaes_small = x_for_eval
optimizer.tell(solutions)
best_features_cmaes_small = [
features_select[i]
for i, val in enumerate(best_sol_raw_cmaes_small.tolist())
if val == 1.0
]
print(f'best objective: {best_objective_cmaes_small}')
print(f'best features: {best_features_cmaes_small}')
CMA-ES优化器定义了对平均值和标准差的一些初始猜测。然后它循环许多代，创建测试点x_for_eval，用目标评估，修改分布(mean, sigma, covariance matrix)等等。每个x_for_eval点都是一个二进制向量[1,1,1,0,0,1，…]]用于从数据集中选择特征。
这里使用的是CMAwM()优化器(带边距的CMA)而不是默认的CMA()。默认的优化器可以很好地处理规则的、连续的问题，但是这里的搜索空间是高维的，并且只允许两个离散值(0和1)。默认优化器会卡在这个空间中。CMAwM()稍微扩大了搜索空间(虽然它返回的解仍然是二进制向量)，可以以解除阻塞。
下图显示了CMA-ES代码寻找最佳解决方案的运行记录。热图显示了每一代每个特征的流行/流行程度(越亮=越受欢迎)。你可以看到有些特征总是很受欢迎，而另一些特征很快就过时了，还有一些特征后来才被“发现”。给定该问题的参数，优化器选择的总体大小为20个点(个体)，因此特征受欢迎程度是在这20个点上平均的。
运行结果如下：
best objective: 33703.070530508514
best generation: 921
objective runs: 20000
time to best: 46.912 sec
它能够找到比SFS更好(更小)的目标值，它调用目标函数的次数更少(20k)，并且花费了大约相同的时间。
在研究了传统的优化算法(遗传算法、模拟退火等)之后，CMA-ES是一个非常好的解决方案，它几乎没有超参数，计算量很轻，它只需要少量的个体(点)来探索搜索空间，但它的性能却很好。如果需要解决优化问题，用它来测试对比是非常有帮助的。
遗传算法
最后我们再看看常用的遗传算法
遗传算法受到生物进化和自然选择的启发。在自然界中，生物(粗略地说)是根据它们所处环境中促进生存和繁殖成功的基因(特征)而被选择的。
对于特征选择，有N个特征并试图找到n个长度的二进制向量[1,0,0,1,1,1，…]，选择特征(1 =特征选择，0 =特征拒绝)，以最小化成本/目标函数。
每个这样的向量可以被认为是一个“个体”。每个向量分量(值0或值1)成为一个“基因”。通过应用进化和选择，有可能进化出一个个体群体，使其接近于我们感兴趣的目标函数的最佳值。
以下是GA的简要介绍。首先生成一群个体(向量)，每个向量的长度为n。向量分量值(基因)是从{0,1}中随机选择的。在下面的图表中，N=12，总体规模为8。
在群体创建后，通过目标函数评估每个个体。
保留客观价值最好的个体，抛弃客观价值最差的个体。这里有许多可能的策略，从单纯的排名(这与直觉相反，并不是很有效)到随机对比选择(从长远来看，这是非常有效的)。还记得“explore-exploit ”困境吗?使用GA，我们很容易陷入这个问题，所以随机对比会比排名好好很多。
一旦最优秀的个体被选择出来，不太适合的个体被抛弃，就可以通过两种技术在基因库中引入变异了:交叉和突变。
交叉的工作原理与自然界完全一样，当两个生物交配并产生后代时:来自父母双方的遗传物质在后代中“混合”，具有一定程度的随机性。
突变也是自然中发生的事情，当遗传物质发生随机突变时，新的价值被引入基因库，增加了它的多样性。
然后个体再次通过目标函数进行评估，选择发生，然后交叉，突变等等。
如果目标函数在某些代数内停止改进，则循环可能会被中断。或者可以为评估的总代数设置一个值，或者运行时间等等，在停止后具有最佳客观价值的个人应该被认为是问题的“解决方案”。
我们可以将GA封装在Optuna中，让Optuna自己寻找最佳的超参数——但这在计算上是非常非常耗时的。
将遗传算法用于特征选择
为了演示，我们使用了非常强大的deep库，如果你不熟悉这个库可能看起来有点困难，我们将尽量保持代码的简单。
# to maximize the objective
# fitness_weights = 1.0
# to minimize the objective
fitness_weights = -1.0
# copy the original dataframes into local copies, once
X_ga = X.copy()
y_ga = y.copy()
# 'const' (the first column) is not an actual feature, do not include it
X_features = X_ga.columns.to_list()[1:]
try:
del creator.FitnessMax
del creator.Individual
except Exception as e:
pass
creator.create("FitnessMax", base.Fitness, weights=(fitness_weights,))
creator.create(
"Individual", array.array, typecode='b', fitness=creator.FitnessMax
)
try:
del toolbox
except Exception as e:
pass
toolbox = base.Toolbox()
# Attribute generator
toolbox.register("attr_bool", random.randint, 0, 1)
# Structure initializers
toolbox.register(
"individual",
tools.initRepeat,
creator.Individual,
toolbox.attr_bool,
len(X_features),
)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
def evalOneMax(individual):
# objective function
# create True/False selector list for features
# and add True at the start for 'const'
cols_select = [True] + [i == 1 for i in list(individual)]
# fit model using the features selected from the individual
lin_mod = sm.OLS(y_ga, X_ga.loc[:, cols_select], hasconst=True).fit()
return (lin_mod.bic,)
toolbox.register("evaluate", evalOneMax)
toolbox.register("mate", tools.cxTwoPoint)
toolbox.register("mutate", tools.mutFlipBit, indpb=0.05)
toolbox.register("select", tools.selTournament, tournsize=3)
random.seed(0)
pop = toolbox.population(n=300)
hof = tools.HallOfFame(1)
pop, log = algorithms.eaSimple(
pop, toolbox, cxpb=0.5, mutpb=0.2, ngen=10, halloffame=hof, verbose=True
)
best_individual_ga_small = list(hof[0])
best_features_ga_small = [
X_features[i] for i, val in enumerate(best_individual_ga_small) if val == 1
]
best_objective_ga_small = (
sm.OLS(y_ga, X_ga[['const'] + best_features_ga_small], hasconst=True)
.fit()
.bic
)
print(f'best objective: {best_objective_ga_small}')
print(f'best features: {best_features_ga_small}')
代码创建了定义个体和整个种群的对象，以及用于评估(目标函数)、交叉、突变和选择的策略。它从300个个体的种群开始，然后调用eaSimple()(一个交叉、突变、选择的固定序列)，我们这里只运行10代，其中绝对最好的个体被保留下来，以免在选择过程中意外突变或被跳过等。
这段简单的代码很容易理解，但效率不高，运行的结果如下：
best objective: 33705.569572544795
best generation: 787
objective runs: 600525
time to best: 157.604 sec
下面的热图显示了各代中每个特征的受欢迎程度(颜色越亮=越受欢迎)。可以看到有些特征总是很受欢迎，有些特征很快就被拒绝了，而另一些特征可能随着时间的推移变得更受欢迎或不那么受欢迎。
方法对比总结
我们尝试了三种不同的技术:SFS、CMA-ES和GA。
这些测试是在AMD Ryzen 7 5800X3D(8/16核)机器上进行的，运行Ubuntu 22.04和Python 3.11.7。SFS和GA是使用有16个线程来运行目标函数。CMA-ES是单进程的——在多线程中运行它似乎并没有提供显著的改进，这可能是算法没有支持多线程，下面是运行时间
SFS: 44.850 sec
GA: 157.604 sec
CMA-ES: 46.912 sec
目标函数调用的次数:
SFS: 22791
GA: 600525
CMA-ES: 20000
目标函数的最佳值-越少越好:
SFS: 33708.9860
GA: 33705.5696
CMA-ES: 33703.0705
CMA-ES找到了最佳目标函数。它的运行时间与SFS相当。它只调用目标函数20k次，是所有方法中调用次数最少的。别忘了，他还是单线程的
GA能够在目标函数上超过SFS，但它是最慢的。它调用目标函数的次数比其他方法多一个数量级。这是因为我们可以认为他是一个半随机的过程，因为遗传突变这个东西没有算法可解释。
SFS速度很快(可以在所有CPU内核上运行)，但性能一般。但它是目前最简单的算法。
如果你只是想用一个简单的算法快速估计出最佳的特征集，那么SFS还不错。如果你想要绝对最好的客观价值，CMA-ES似乎是首选，并且它也不慢。
最后本文的代码：
https://avoid.overfit.cn/post/3bd329a18abe4ab4af36e6b7ceaef469
作者：Florin Andrei
