在建立模型时，特征选择是一个重要环节，它指通过保留一部分特征子集来拟合模型，而舍弃其余特征。进行特征选择有多重原因:
- 保持模型的可解释性(过多特征会增加解释难度)
- 避免维数灾难
- 优化与模型相关的目标函数(如R平方、AIC等)
- 防止过拟合等
如果特征数量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 np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import os
import random
import copy
from itertools import repeat
import shutil
import time
import math
import array
df = pd.read_csv('train.csv')
# drop columns with NaN
df.dropna(axis=1, inplace=True)
# drop irrelevant columns
df.drop(columns=['Id'], inplace=True)
# drop major outliers
df = 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 first
df_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)
顺序特征搜索是一种贪婪的特征选择算法，它包含前向搜索和后向搜索两种策略。以前向搜索为例，算法流程如下:
- 首先从全部N个特征中选择一个使目标函数值最优的单特征子集。
- 在已选特征子集的基础上，再添加一个新特征，形成两个特征的子集，选择能使目标函数进一步最小化的那个组合。
- 重复步骤2，每轮迭代都只增加一个新特征，直到所有特征都被尝试加入子集。
- 从所有被尝试过的特征子集中，选择使目标函数值最小的那个作为最终输出。
SFS是一种贪婪算法，它每一步的选择都是基于当前最优解的局部决策，无法回头修正之前的决策。但它的搜索复杂度只有O(N^2)，相比暴力搜索指数级的复杂度，计算效率大幅提高。当然，这种高效是以潜在的全局最优解被忽略为代价的。
SFS的后向策略则是从全量特征集合出发，每轮迭代移除一个使目标函数值改善最小的特征，直至所有特征被遍历过为止。
使用 mlxtend 库[2] 编写的 SFS 代码在 Python 中是什么样的：
import statsmodels.api as sm
from mlxtend.feature_selection import SequentialFeatureSelector as SFS
from sklearn.base import BaseEstimator
class DummyEstimator(BaseEstimator):
# mlxtend 希望使用 sklearn 估计器，但这里不需要（使用 statsmodels OLS 代替）。
def fit(self, X, y=None, **kwargs):
return self
def neg_bic(m, X, y):
# objective function
lin_mod_res = sm.OLS(y, X, hascnotallow=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: 42.448 sec在 213 个特征中，最佳特征数为 36 个。最佳 BIC 为 33708.986（特征选择前的基线值为 34570.166），在我的系统上完成这一过程用时不到 1 分钟。它调用目标函数 22.8k 次。
这些是最佳 BIC 值和 R 方值与所选特征数量的函数关系：
best_objective_seq = -np.inf
r2_of_best_k = 0
r2_list = []
best_k = 1
best_features_seq = []
# also extract R2 from the feature selection search
for 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 = r2
print(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: 36
best objective: 33708.98602877906
R2 @ best k: 0.9075677543649224
objective runs: 22791
total run time: 42.448 sec
best 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通过自适应策略有效利用样本信息，在保证全局性的同时提高了收敛速度。
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个点的群体，这是优化器选择的默认群体大小。这是第一步。接下来，算法需要：
- 测算每个点的目标函数（Rastrigin）
- 根据从目标函数中获得的知识，更新均值、标准差和协方差矩阵，有效地创建一个新的多元正态分布
- 使用新的分布产生一组新的测试点
- 重复这个过程，直到达到某个标准（要么收敛到某个平均值，要么超过最大步数等）。
仅仅更新分布的平均值是非常简单的。工作原理如下：计算每个测试点的目标函数后，给这些点分配权重，目标值较高的点权重较大，然后根据它们的位置计算出加权和，这就是新的平均值。实际上，CMA-ES（协方差矩阵自适应演化策略）将分布均值向目标值较好的点移动。
更新 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.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 优化器的定义涉及对均值和 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, once
X_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_features
mean = np.full(n_features, 0.5)
sigma = 1 / 6
pop_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 = 1000
best_objective_cmaes = np.inf
best_generation_cmaes = 0
best_sol_raw_cmaes = None
history_values_cmaes = np.full((gen_max,), np.nan)
history_values_best_cmaes = np.full((gen_max,), np.nan)
time_to_best_cmaes = np.inf
objective_runs_cmaes = 0
solutions_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()
break
gen_completed_cmaes = generation
best_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.070530508514
best generation: 921
objective runs: 20000
time to best: 48.326 sec
best 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.070530508514
best generation: 921
objective runs: 20000
time 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
