在建立模型时，特征选择是一个重要环节，它指通过保留一部分特征子集来拟合模型，而舍弃其余特征。进行特征选择有多重原因:
如果特征数量N较小，可使用穷举搜索尝试所有可能的特征组合，保留使成本/目标函数最小的那个。但当N较大时，穷举搜索就行不通了，因为需尝试的组合数为2^N，这是指数级增长，N超过几十个就变得极其耗时。
此时需采用启发式算法，以有效方式探索搜索空间，寻找能使目标函数最小化的特征组合。具体来说，需寻找一个长度为N的0/1向量[1,0,0,1,0,1,1,0,...]，其中1表示选择该特征，0表示舍弃。目标是找到一个能最小化目标函数的这样一个向量。搜索空间的维度等于特征数量N，每一维只有0/1两种取值可能。
寻找一个好的启发式算法并非易事。R语言中regsubsets()函数提供了一些选项。Scikit-learn库也提供了几种启发式特征选择方法，前提是问题能很好地符合它们的技术假设。然而，要找到一种通用的、高效的启发式算法仍是一个挑战。在本系列文章中，我们将探讨几种即使在特征数量N很大、目标函数可为任意可计算函数(只要不过于缓慢)的情况下，也能给出合理结果的协方差矩阵适应进化算法方法。
特征选择是机器学习中一个重要的预处理步骤,它通过选择出对预测目标贡献最大的特征子集,从而提高模型的性能和简洁性。常见的特征选择方法包括Filter(过滤式)、Wrapper(包装式)和Embedded(嵌入式)等。其中,协方差矩阵适应进化算法(Covariance Matrix Adaptation Evolution Strategy, CMA-ES)是一种高效的Wrapper式特征选择算法。
在本文中，我们将使用著名的房价预测数据集(来自Kaggle[1] ，共213个特征，1453个样本)对CMA-ES算法的特征选择能力进行说明。我们所使用的模型是线性回归模型，目标是最小化贝叶斯信息准则(BIC)，它是一种评估模型质量的指标，值越小表示模型越好。与之类似的指标还有AIC(Akaike信息准则)，两者都能有效避免过拟合。当然，我们也可以使用R平方或调整R平方作为目标函数，只是需要注意此时目标是最大化而非最小化。
图片
import numpy as npimport pandas as pdimport matplotlib.pyplot as pltimport seaborn as snsimport osimport randomimport copyfrom itertools import repeatimport shutilimport timeimport mathimport arraydf = pd.read_csv('train.csv')# drop columns with NaNdf.dropna(axis=1, inplace=True)# drop irrelevant columnsdf.drop(columns=['Id'], inplace=True)# drop major outliersdf = df[df['BsmtFinSF1'] <= 5000]df = df[df['MiscVal'] <= 5000]df = df[df['LotArea'] <= 100000]columns_non_categorical = df.select_dtypes(exclude='object').columns.to_list()columns_non_categorical.sort()columns_non_categorical.remove('SalePrice')columns_non_categorical = ['SalePrice'] + columns_non_categorical# alphabetically sort columns, keep target firstdf_temp = df[['SalePrice']]df.drop(columns=['SalePrice'], inplace=True)df.sort_index(axis=1, inplace=True)df = pd.concat([df_temp, df], axis=1)del df_temp无论选择何种评估指标作为目标函数，CMA-ES算法都能通过迭代搜索的方式，找到一个使目标函数值最优的特征子集，从而实现自动高效的特征选择。下面将对算法的原理和应用进行详细阐述。
我们将尝试通过特征选择来最小化 BIC，因此这里是在启用所有特征选择之前，从 statsmodels.api.OLS() 中得到的 BIC 基准值：
X = df.drop(columns=['SalePrice']).copy()y = df['SalePrice'].copy()X = add_constant(X)linear_model = sm.OLS(y, X).fit()print(f'baseline BIC: {linear_model.bic}')baseline BIC: 34570.166173470934现在，让我们来看看一种著名的、久经考验的特征选择技术，我们将把它与后面介绍的更复杂的技术进行比较。
顺序特征搜索是一种贪婪的特征选择算法，它包含前向搜索和后向搜索两种策略。以前向搜索为例，算法流程如下:
SFS是一种贪婪算法，它每一步的选择都是基于当前最优解的局部决策，无法回头修正之前的决策。但它的搜索复杂度只有O(N^2)，相比暴力搜索指数级的复杂度，计算效率大幅提高。当然，这种高效是以潜在的全局最优解被忽略为代价的。
SFS的后向策略则是从全量特征集合出发，每轮迭代移除一个使目标函数值改善最小的特征，直至所有特征被遍历过为止。
使用 mlxtend 库[2] 编写的 SFS 代码在 Python 中是什么样的：
import statsmodels.api as smfrom mlxtend.feature_selection import SequentialFeatureSelector as SFSfrom sklearn.base import BaseEstimatorclass DummyEstimator(BaseEstimator): # mlxtend 希望使用 sklearn 估计器，但这里不需要（使用 statsmodels OLS 代替）。 def fit(self, X, y=None, **kwargs): return selfdef neg_bic(m, X, y): # objective function lin_mod_res = sm.OLS(y, X, hascnotallow=True).fit() return -lin_mod_res.bicseq_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] - 1objective_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_seqseq_metrics = seq_res.get_metric_dict()它可以快速运行各种组合，这些就是汇总结果：
best k: 36best objective: 33708.98602877906R2 @ best k: 0.9075677543649224objective runs: 22791total run time: 42.448 sec在 213 个特征中，最佳特征数为 36 个。最佳 BIC 为 33708.986（特征选择前的基线值为 34570.166），在我的系统上完成这一过程用时不到 1 分钟。它调用目标函数 22.8k 次。
这些是最佳 BIC 值和 R 方值与所选特征数量的函数关系：
best_objective_seq = -np.infr2_of_best_k = 0r2_list = []best_k = 1best_features_seq = []# also extract R2 from the feature selection searchfor k in tqdm(seq_metrics.keys()): r2_eval_mod = sm.OLS( y, X[list(seq_metrics[k]['feature_names'])], hascnotallow=True ) r2_eval_mod_res = r2_eval_mod.fit() r2 = r2_eval_mod_res.rsquared r2_list.append(r2) score_k = seq_metrics[k]['avg_score'] if score_k > best_objective_seq: best_objective_seq = score_k best_k = k best_features_seq = list(seq_metrics[k]['feature_names']) r2_of_best_k = r2print(f'best k: {best_k}')print(f'best objective: {-best_objective_seq}')print(f'R2 @ best k: {r2_of_best_k}')print(f'objective runs: {objective_runs_sfs}')print(f'total run time: {run_time_seq:.3f} sec')print(f'best features: {best_features_seq}')sfs_avg = [-seq_metrics[k]['avg_score'] for k in sorted(seq_metrics.keys())]fig, ax = plt.subplots(1, 2, figsize=(10, 4), layout='constrained')ax[0].scatter(sorted(seq_metrics.keys()), sfs_avg, s=1)ax[0].set_xticks(sorted(seq_metrics.keys()), minor=True)ax[0].set_title('BIC')ax[0].set_xlabel('number of features')ax[0].set_ylabel('BIC')ax[1].scatter(sorted(seq_metrics.keys()), r2_list, s=1)ax[1].set_title('R2')ax[1].set_xlabel('number of features')ax[1].set_ylabel('R2')fig.suptitle('sequential forward search')fig.savefig('sfs-performance.png')fig.show()best k: 36best objective: 33708.98602877906R2 @ best k: 0.9075677543649224objective runs: 22791total run time: 42.448 secbest features: ['const', 'BedroomAbvGr', 'BldgType_1Fam', 'BsmtFinSF1', 'Condition1_Norm', 'Condition2_PosN', 'ExterQual_Ex', 'Exterior1st_BrkFace', 'Functional_Sev', 'Functional_Typ', 'GarageArea', 'GrLivArea', 'HeatingQC_Ex', 'KitchenQual_Ex', 'LandContour_HLS', 'LotArea', 'LotConfig_CulDSac', 'LowQualFinSF', 'MSSubClass', 'Neighborhood_BrkSide', 'Neighborhood_Crawfor', 'Neighborhood_NWAmes', 'Neighborhood_NoRidge', 'Neighborhood_NridgHt', 'Neighborhood_Somerst', 'Neighborhood_StoneBr', 'OverallCond', 'OverallQual', 'PoolArea', 'RoofMatl_WdShngl', 'SaleCondition_Normal', 'SaleType_New', 'ScreenPorch', 'TotalBsmtSF', 'WoodDeckSF', 'YearBuilt']
BIC and R-squared for SFS
SFS 的 BIC 和 R 平方
有关实际选择的特征名称：
'const', 'BedroomAbvGr', 'BldgType_1Fam', 'BsmtFinSF1', 'Condition1_Norm', 'Condition2_PosN', 'ExterQual_Ex', 'Exterior1st_BrkFace', 'Functional_Sev', 'Functional_Typ', 'GarageArea', 'GrLivArea', 'HeatingQC_Ex', 'KitchenQual_Ex', 'LandContour_HLS', 'LotArea', 'LotConfig_CulDSac', 'LowQualFinSF', 'MSSubClass', 'Neighborhood_BrkSide', 'Neighborhood_Crawfor', 'Neighborhood_NWAmes', 'Neighborhood_NoRidge', 'Neighborhood_NridgHt', 'Neighborhood_Somerst', 'Neighborhood_StoneBr', 'OverallCond', 'OverallQual', 'PoolArea', 'RoofMatl_WdShngl', 'SaleCondition_Normal', 'SaleType_New', 'ScreenPorch', 'TotalBsmtSF', 'WoodDeckSF', 'YearBuilt'CMA-ES是一种高效的无约束/约束黑箱连续优化的随机搜索算法。它属于进化计算的一种,但与传统的遗传算法有着明显区别。
与遗传算法直接对解个体进行变异和交叉操作不同，CMA-ES在连续域上对多元正态分布模型的参数(均值和协方差矩阵)进行更新迭代，间接实现对潜在解集群的适应性搜索。
算法不需要目标函数的梯度信息，即无需计算导数，因此具有很强的鲁棒性，可应用于非线性、非凸、多峰、多模态等各种复杂优化问题。同时，CMA-ES通过自适应策略有效利用样本信息，在保证全局性的同时提高了收敛速度。
CMA-ES已被广泛应用于机器学习、计算机视觉、计算生物学等诸多领域，并成为优选的优化算法之一。在Optuna、PyGMO等知名数值优化库中都有CMA-ES的实现版本。由于算法理论较为复杂，这里只是简要介绍，可参考文末的扩展阅读材料进一步学习。
考虑二维 Rastrigin 函数：
Rastrigin 函数
下面的热图显示了该函数的值--颜色越亮表示值越高。该函数的全局最小值位于原点（0，0），但其中有许多局部极值。我们需要通过 CMA-ES 找出全局最小值。
Rastrigin 函数热图
CMA-ES利用多元正态分布作为基础。它根据该分布在搜索空间中生成测试点。用户需要猜测分布的原始均值和标准偏差，然后算法会反复调整这些参数，并在搜索空间中扫描分布，以寻找最佳的目标函数值。以下是测试点的初始分布：
CMA-ES 分布
xi 表示算法在搜索空间中产生的每一步的点集。lambda 是产生的点数。该分布的平均值在每一步中都会更新，最终希望收敛到真正的解决方案。sigma 是分布的标准偏差，即测试点的分布。C 是协方差矩阵，它定义了分布的形状。根据 C 的值，分布可能是“圆形”，也可能是拉长的椭圆形。改变 C 可以让CMA-ES进入搜索空间的某些区域，或者避开其他区域。
第一代测试点
在图中创建了一个包含6个点的群体，这是优化器选择的默认群体大小。这是第一步。接下来，算法需要：
仅仅更新分布的平均值是非常简单的。工作原理如下：计算每个测试点的目标函数后，给这些点分配权重，目标值较高的点权重较大，然后根据它们的位置计算出加权和，这就是新的平均值。实际上，CMA-ES（协方差矩阵自适应演化策略）将分布均值向目标值较好的点移动。
更新 CMA-ES 分布均值
如果算法达到真实解决方案，分布的平均值将趋于该解决方案。协方差矩阵将导致分布的形状发生变化（圆形或椭圆形），这取决于目标函数的地理位置，会向有利的区域扩展，而回避不利的区域。
二维Rastrigin函数相对简单，因为它只有2个维度。然而，对于我们的特征选择问题，我们面临N=213个维度，并且空间并不是连续的。每个测试点都是一个N维向量，分量的值为"{0, 1}"。换句话说，每个测试点看起来像这样：1，0，0，1，1，0，......一个二进制向量。除此之外，问题是相同的：我们需要找到使目标函数（即OLS模型的BIC参数）最小化的点或向量。
对于具有 213 个维度和离散取值（0 或 1）的特征选择问题，您可以考虑使用遗传算法或者模拟退火算法来寻找最优解。这些算法对于高维度和离散空间的优化问题非常有效。
可以将特征选择问题建模为一个多维的优化问题，然后使用上述算法来寻找最小化目标函数的解。这些算法可以帮助你在高维度和离散空间中寻找较优的特征子集。
下面是使用 cmaes 库[3] 进行特征选择的 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], hascnotallow=True).fit() return lin_mod.bicX_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 = 100best_objective_cmaes_small = np.infbest_sol_raw_cmaes_small = Nonefor 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 优化器的定义涉及对均值和 sigma（标准偏差）进行一些初始猜测。然后，优化器会循环运行多代，创建测试点 x_for_eval ，并根据目标评估其，然后修改分布（均值、sigma、协方差矩阵）等。每个x_for_eval点都是一个二进制向量[1, 1, 1, 0, 0, 1, ...]，用于从数据集中选择特征。
请注意，这里使用的是 CMAwM() 优化器（带边际的 CMA），而不是默认的 CMA()。默认优化器在处理常规连续问题时效果很好，但这里的搜索空间是高维的，只允许两个离散值（0 和 1）。默认优化器就会卡在这个空间里。CMAwM()扩大了一点搜索空间（而它返回的解仍然是二进制向量），这似乎足以解除它的阻碍。
这段简单的代码确实有效，但还远未达到最佳状态。下面是一个更复杂、更优化的版本，它能更快地找到更好的解。
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], hascnotallow=True).fit() return lin_mod.bic# copy the original dataframes into local copies, onceX_cmaes = X.copy()y_cmaes = y.copy()rs = np.random.RandomState(seed=0)features_select = [f for f in X_cmaes.columns if f != 'const']n_features = len(features_select)cma_bounds = np.tile([0, 1], (n_features, 1))cma_steps = np.ones(n_features)n_max_resampling = 10 * n_featuresmean = np.full(n_features, 0.5)sigma = 1 / 6pop_size = 4 + math.floor(3 * math.log(n_features))margin = 1 / (n_features * pop_size)optimizer = CMAwM( mean=mean, sigma=sigma, bounds=cma_bounds, steps=cma_steps, n_max_resampling=n_max_resampling, seed=0, population_size=pop_size, margin=margin,)gen_max = 1000best_objective_cmaes = np.infbest_generation_cmaes = 0best_sol_raw_cmaes = Nonehistory_values_cmaes = np.full((gen_max,), np.nan)history_values_best_cmaes = np.full((gen_max,), np.nan)time_to_best_cmaes = np.infobjective_runs_cmaes = 0solutions_avg = np.full((gen_max, n_features), np.nan)n_jobs = os.cpu_count()iterator = tqdm(range(gen_max), desc='generation')t_start_cmaes = time.time()for generation in iterator: best_value_gen = np.inf # solutions fed back to the optimizer solutions_float = [] # binary-truncated solutions - the yes/no answers we're looking for solutions_binary = np.full((pop_size, n_features), np.nan) vals = np.full((pop_size,), np.nan) for i in range(pop_size): # re-seeding the RNG is very important # otherwise CMA-ES gets stuck in local extremes seed = rs.randint(1, 2**16) + generation optimizer._rng.seed(seed) fs_for_eval, fs_for_tell = optimizer.ask() solutions_binary[i,] = fs_for_eval value = cma_objective(fs_for_eval) objective_runs_cmaes += 1 vals[i] = value solutions_float.append((fs_for_tell, value)) optimizer.tell(solutions_float) solutions_avg[generation,] = solutions_binary.mean(axis=0) best_value_gen = vals.min() t_end_loop_cmaes = time.time() if best_value_gen < best_objective_cmaes: best_objective_cmaes = best_value_gen best_generation_cmaes = generation best_sol_raw_cmaes = solutions_binary[np.argmin(vals),] time_to_best_cmaes = t_end_loop_cmaes - t_start_cmaes print( f'gen: {best_generation_cmaes:5n}, new best objective: {best_objective_cmaes:.4f}' ) history_values_cmaes[generation] = best_value_gen history_values_best_cmaes[generation] = best_objective_cmaes if optimizer.should_stop(): print(f'Optimizer decided to stop.') iterator.close() break if os.path.isfile('break'): # to gracefully break the loop, manually create a file called 'break' print(f'Found break file, stopping now.') iterator.close() breakgen_completed_cmaes = generationbest_features_cmaes = [ features_select[i] for i, val in enumerate(best_sol_raw_cmaes.tolist()) if val == 1.0]print(f'best objective: {best_objective_cmaes}')print(f'best generation: {best_generation_cmaes}')print(f'objective runs: {objective_runs_cmaes}')print(f'time to best: {time_to_best_cmaes:.3f} sec')print(f'best features: {best_features_cmaes}')cma_mv_df = pd.DataFrame(solutions_avg, columns=features_select).iloc[ : gen_completed_cmaes + 1]if gen_completed_cmaes > 20: x_ticks = list( range(0, gen_completed_cmaes, round(gen_completed_cmaes / 20)) )else: x_ticks = list(range(0, gen_completed_cmaes))if x_ticks[-1] != gen_completed_cmaes: x_ticks.append(gen_completed_cmaes)fig, ax = plt.subplots( 2, 1, sharex=True, height_ratios=[20, 1], figsize=(12, 45), layout='constrained',)sns.heatmap( cma_mv_df.T, vmin=0.0, vmax=1.0, cmap='viridis', cbar=True, cbar_kws={'fraction': 0.01, 'anchor': (0.0, 1.0)}, ax=ax[0],)ax[0].axvline(x=best_generation_cmaes, color='C7')ax[0].tick_params(axis='both', reset=True)ax[0].set_xticks(x_ticks)ax[0].set_xticklabels(x_ticks)ax[0].set_title('Population average of feature selector values')ax[0].set_xlabel('generation')ax[1].scatter( range(gen_completed_cmaes + 1), history_values_cmaes[: gen_completed_cmaes + 1], s=1, label='current value',)ax[1].plot( range(gen_completed_cmaes + 1), history_values_best_cmaes[: gen_completed_cmaes + 1], color='C1', label='best value',)ax[1].vlines( x=best_generation_cmaes, ymin=history_values_cmaes[: gen_completed_cmaes + 1].min(), ymax=history_values_cmaes[: gen_completed_cmaes + 1].max(), colors='C7',)ax[1].tick_params(axis='both', reset=True)ax[1].legend()ax[1].set_title('Objective value')ax[1].set_xlabel('generation')fig.suptitle('CMA-ES')fig.savefig('cmaes-performance.png')fig.show()best objective: 33703.070530508514best generation: 921objective runs: 20000time to best: 48.326 secbest features: ['BedroomAbvGr', 'BldgType_1Fam', 'BsmtFinSF1', 'Condition1_Norm', 'Condition2_PosN', 'ExterQual_Ex', 'Exterior1st_BrkFace', 'Functional_Sev', 'Functional_Typ', 'GarageCars', 'GrLivArea', 'HeatingQC_Ex', 'KitchenQual_Ex', 'LandContour_HLS', 'LotArea', 'LowQualFinSF', 'Neighborhood_BrkSide', 'Neighborhood_Crawfor', 'Neighborhood_NWAmes', 'Neighborhood_NoRidge', 'Neighborhood_NridgHt', 'Neighborhood_Somerst', 'Neighborhood_StoneBr', 'OverallCond', 'OverallQual', 'PoolArea', 'RoofMatl_WdShngl', 'SaleCondition_Normal', 'SaleType_New', 'ScreenPorch', 'TotalBsmtSF', 'WoodDeckSF', 'YearBuilt']下图显示了复杂、优化的 CMA-ES 代码搜索最佳解决方案的历史。热图显示了每一代中每种功能的流行程度（更亮 = 更流行）。可以看到，有些特征总是非常流行，而另一些则很快过时，还有一些则是后来才 “发现 ”的。根据这个问题的参数，优化器选择的群体大小为 20 个点（个体），因此特征流行度是这 20 个点的平均值。
上下滑动查看更多CMA-ES 优化历史这些是经过优化的 CMA-ES 代码的主要统计信息：
best objective: 33703.070530508514best generation: 921objective runs: 20000time to best: 48.326 sec与 SFS 相比，它能找到更好（更小）的目标值，调用目标函数的次数更少（20k），所用时间也差不多。从所有指标来看，这都是比 SFS 的净收益。
同样，特征选择前的基准 BIC 为
baseline BIC: 34570.166173470934顺便提一句：在研究了传统优化算法（遗传算法、模拟退火等）之后，CMA-ES 给我带来了惊喜。它几乎没有超参数，计算量小，只需要少量个体（点）来探索搜索空间，而且性能相当不错。如果你需要解决优化问题，值得把它放在你的工具箱里。
[1]Kaggle: https://www.kaggle.com/c/house-prices-advanced-regression-techniques/data
[2]mlxtend 库: https://rasbt.github.io/mlxtend/
[3]cmaes 库: https://github.com/CyberAgentAILab/cmaes
本文链接：http://www.28at.com/showinfo-26-101117-0.html协方差矩阵适应进化算法实现高效特征选择
声明：本网页内容旨在传播知识，若有侵权等问题请及时与本网联系，我们将在第一时间删除处理。邮件：2376512515@qq.com
Copyright © 2016-2023 天津谷骐科技有限公司 版权所有 sitemap.xml
违法及侵权请联系：2376512515@qq.com 津ICP备18001702号
津公网安备 12010102000574号
