在建立模型時(shí),特征選擇是一個(gè)重要環(huán)節(jié),它指通過保留一部分特征子集來(lái)擬合模型,而舍棄其余特征。進(jìn)行特征選擇有多重原因:
如果特征數(shù)量N較小,可使用窮舉搜索嘗試所有可能的特征組合,保留使成本/目標(biāo)函數(shù)最小的那個(gè)。但當(dāng)N較大時(shí),窮舉搜索就行不通了,因?yàn)樾鑷L試的組合數(shù)為2^N,這是指數(shù)級(jí)增長(zhǎng),N超過幾十個(gè)就變得極其耗時(shí)。
此時(shí)需采用啟發(fā)式算法,以有效方式探索搜索空間,尋找能使目標(biāo)函數(shù)最小化的特征組合。具體來(lái)說(shuō),需尋找一個(gè)長(zhǎng)度為N的0/1向量[1,0,0,1,0,1,1,0,...],其中1表示選擇該特征,0表示舍棄。目標(biāo)是找到一個(gè)能最小化目標(biāo)函數(shù)的這樣一個(gè)向量。搜索空間的維度等于特征數(shù)量N,每一維只有0/1兩種取值可能。
尋找一個(gè)好的啟發(fā)式算法并非易事。R語(yǔ)言中regsubsets()函數(shù)提供了一些選項(xiàng)。Scikit-learn庫(kù)也提供了幾種啟發(fā)式特征選擇方法,前提是問題能很好地符合它們的技術(shù)假設(shè)。然而,要找到一種通用的、高效的啟發(fā)式算法仍是一個(gè)挑戰(zhàn)。在本系列文章中,我們將探討幾種即使在特征數(shù)量N很大、目標(biāo)函數(shù)可為任意可計(jì)算函數(shù)(只要不過于緩慢)的情況下,也能給出合理結(jié)果的協(xié)方差矩陣適應(yīng)進(jìn)化算法方法。
特征選擇是機(jī)器學(xué)習(xí)中一個(gè)重要的預(yù)處理步驟,它通過選擇出對(duì)預(yù)測(cè)目標(biāo)貢獻(xiàn)最大的特征子集,從而提高模型的性能和簡(jiǎn)潔性。常見的特征選擇方法包括Filter(過濾式)、Wrapper(包裝式)和Embedded(嵌入式)等。其中,協(xié)方差矩陣適應(yīng)進(jìn)化算法(Covariance Matrix Adaptation Evolution Strategy, CMA-ES)是一種高效的Wrapper式特征選擇算法。
在本文中,我們將使用著名的房?jī)r(jià)預(yù)測(cè)數(shù)據(jù)集(來(lái)自Kaggle[1] ,共213個(gè)特征,1453個(gè)樣本)對(duì)CMA-ES算法的特征選擇能力進(jìn)行說(shuō)明。我們所使用的模型是線性回歸模型,目標(biāo)是最小化貝葉斯信息準(zhǔn)則(BIC),它是一種評(píng)估模型質(zhì)量的指標(biāo),值越小表示模型越好。與之類似的指標(biāo)還有AIC(Akaike信息準(zhǔn)則),兩者都能有效避免過擬合。當(dāng)然,我們也可以使用R平方或調(diào)整R平方作為目標(biāo)函數(shù),只是需要注意此時(shí)目標(biāo)是最大化而非最小化。
圖片
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
無(wú)論選擇何種評(píng)估指標(biāo)作為目標(biāo)函數(shù),CMA-ES算法都能通過迭代搜索的方式,找到一個(gè)使目標(biāo)函數(shù)值最優(yōu)的特征子集,從而實(shí)現(xiàn)自動(dòng)高效的特征選擇。下面將對(duì)算法的原理和應(yīng)用進(jìn)行詳細(xì)闡述。
我們將嘗試通過特征選擇來(lái)最小化 BIC,因此這里是在啟用所有特征選擇之前,從 statsmodels.api.OLS() 中得到的 BIC 基準(zhǔn)值:
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
現(xiàn)在,讓我們來(lái)看看一種著名的、久經(jīng)考驗(yàn)的特征選擇技術(shù),我們將把它與后面介紹的更復(fù)雜的技術(shù)進(jìn)行比較。
順序特征搜索是一種貪婪的特征選擇算法,它包含前向搜索和后向搜索兩種策略。以前向搜索為例,算法流程如下:
SFS是一種貪婪算法,它每一步的選擇都是基于當(dāng)前最優(yōu)解的局部決策,無(wú)法回頭修正之前的決策。但它的搜索復(fù)雜度只有O(N^2),相比暴力搜索指數(shù)級(jí)的復(fù)雜度,計(jì)算效率大幅提高。當(dāng)然,這種高效是以潛在的全局最優(yōu)解被忽略為代價(jià)的。
SFS的后向策略則是從全量特征集合出發(fā),每輪迭代移除一個(gè)使目標(biāo)函數(shù)值改善最小的特征,直至所有特征被遍歷過為止。
使用 mlxtend 庫(kù)[2] 編寫的 SFS 代碼在 Python 中是什么樣的:
import statsmodels.api as smfrom mlxtend.feature_selection import SequentialFeatureSelector as SFSfrom sklearn.base import BaseEstimatorclass DummyEstimator(BaseEstimator): # mlxtend 希望使用 sklearn 估計(jì)器,但這里不需要(使用 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()
它可以快速運(yùn)行各種組合,這些就是匯總結(jié)果:
best k: 36best objective: 33708.98602877906R2 @ best k: 0.9075677543649224objective runs: 22791total run time: 42.448 sec
在 213 個(gè)特征中,最佳特征數(shù)為 36 個(gè)。最佳 BIC 為 33708.986(特征選擇前的基線值為 34570.166),在我的系統(tǒng)上完成這一過程用時(shí)不到 1 分鐘。它調(diào)用目標(biāo)函數(shù) 22.8k 次。
這些是最佳 BIC 值和 R 方值與所選特征數(shù)量的函數(shù)關(guān)系:
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 平方
有關(guān)實(shí)際選擇的特征名稱:
'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是一種高效的無(wú)約束/約束黑箱連續(xù)優(yōu)化的隨機(jī)搜索算法。它屬于進(jìn)化計(jì)算的一種,但與傳統(tǒng)的遺傳算法有著明顯區(qū)別。
與遺傳算法直接對(duì)解個(gè)體進(jìn)行變異和交叉操作不同,CMA-ES在連續(xù)域上對(duì)多元正態(tài)分布模型的參數(shù)(均值和協(xié)方差矩陣)進(jìn)行更新迭代,間接實(shí)現(xiàn)對(duì)潛在解集群的適應(yīng)性搜索。
算法不需要目標(biāo)函數(shù)的梯度信息,即無(wú)需計(jì)算導(dǎo)數(shù),因此具有很強(qiáng)的魯棒性,可應(yīng)用于非線性、非凸、多峰、多模態(tài)等各種復(fù)雜優(yōu)化問題。同時(shí),CMA-ES通過自適應(yīng)策略有效利用樣本信息,在保證全局性的同時(shí)提高了收斂速度。
CMA-ES已被廣泛應(yīng)用于機(jī)器學(xué)習(xí)、計(jì)算機(jī)視覺、計(jì)算生物學(xué)等諸多領(lǐng)域,并成為優(yōu)選的優(yōu)化算法之一。在Optuna、PyGMO等知名數(shù)值優(yōu)化庫(kù)中都有CMA-ES的實(shí)現(xiàn)版本。由于算法理論較為復(fù)雜,這里只是簡(jiǎn)要介紹,可參考文末的擴(kuò)展閱讀材料進(jìn)一步學(xué)習(xí)。
考慮二維 Rastrigin 函數(shù):
Rastrigin 函數(shù)
下面的熱圖顯示了該函數(shù)的值--顏色越亮表示值越高。該函數(shù)的全局最小值位于原點(diǎn)(0,0),但其中有許多局部極值。我們需要通過 CMA-ES 找出全局最小值。
Rastrigin 函數(shù)熱圖
CMA-ES利用多元正態(tài)分布作為基礎(chǔ)。它根據(jù)該分布在搜索空間中生成測(cè)試點(diǎn)。用戶需要猜測(cè)分布的原始均值和標(biāo)準(zhǔn)偏差,然后算法會(huì)反復(fù)調(diào)整這些參數(shù),并在搜索空間中掃描分布,以尋找最佳的目標(biāo)函數(shù)值。以下是測(cè)試點(diǎn)的初始分布:
CMA-ES 分布
xi 表示算法在搜索空間中產(chǎn)生的每一步的點(diǎn)集。lambda 是產(chǎn)生的點(diǎn)數(shù)。該分布的平均值在每一步中都會(huì)更新,最終希望收斂到真正的解決方案。sigma 是分布的標(biāo)準(zhǔn)偏差,即測(cè)試點(diǎn)的分布。C 是協(xié)方差矩陣,它定義了分布的形狀。根據(jù) C 的值,分布可能是“圓形”,也可能是拉長(zhǎng)的橢圓形。改變 C 可以讓CMA-ES進(jìn)入搜索空間的某些區(qū)域,或者避開其他區(qū)域。
第一代測(cè)試點(diǎn)
在圖中創(chuàng)建了一個(gè)包含6個(gè)點(diǎn)的群體,這是優(yōu)化器選擇的默認(rèn)群體大小。這是第一步。接下來(lái),算法需要:
僅僅更新分布的平均值是非常簡(jiǎn)單的。工作原理如下:計(jì)算每個(gè)測(cè)試點(diǎn)的目標(biāo)函數(shù)后,給這些點(diǎn)分配權(quán)重,目標(biāo)值較高的點(diǎn)權(quán)重較大,然后根據(jù)它們的位置計(jì)算出加權(quán)和,這就是新的平均值。實(shí)際上,CMA-ES(協(xié)方差矩陣自適應(yīng)演化策略)將分布均值向目標(biāo)值較好的點(diǎn)移動(dòng)。
更新 CMA-ES 分布均值
如果算法達(dá)到真實(shí)解決方案,分布的平均值將趨于該解決方案。協(xié)方差矩陣將導(dǎo)致分布的形狀發(fā)生變化(圓形或橢圓形),這取決于目標(biāo)函數(shù)的地理位置,會(huì)向有利的區(qū)域擴(kuò)展,而回避不利的區(qū)域。
二維Rastrigin函數(shù)相對(duì)簡(jiǎn)單,因?yàn)樗挥?個(gè)維度。然而,對(duì)于我們的特征選擇問題,我們面臨N=213個(gè)維度,并且空間并不是連續(xù)的。每個(gè)測(cè)試點(diǎn)都是一個(gè)N維向量,分量的值為"{0, 1}"。換句話說(shuō),每個(gè)測(cè)試點(diǎn)看起來(lái)像這樣:1,0,0,1,1,0,......一個(gè)二進(jìn)制向量。除此之外,問題是相同的:我們需要找到使目標(biāo)函數(shù)(即OLS模型的BIC參數(shù))最小化的點(diǎn)或向量。
對(duì)于具有 213 個(gè)維度和離散取值(0 或 1)的特征選擇問題,您可以考慮使用遺傳算法或者模擬退火算法來(lái)尋找最優(yōu)解。這些算法對(duì)于高維度和離散空間的優(yōu)化問題非常有效。
可以將特征選擇問題建模為一個(gè)多維的優(yōu)化問題,然后使用上述算法來(lái)尋找最小化目標(biāo)函數(shù)的解。這些算法可以幫助你在高維度和離散空間中尋找較優(yōu)的特征子集。
下面是使用 cmaes 庫(kù)[3] 進(jìn)行特征選擇的 CMA-ES 代碼的簡(jiǎn)單版本:
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 優(yōu)化器的定義涉及對(duì)均值和 sigma(標(biāo)準(zhǔn)偏差)進(jìn)行一些初始猜測(cè)。然后,優(yōu)化器會(huì)循環(huán)運(yùn)行多代,創(chuàng)建測(cè)試點(diǎn) x_for_eval ,并根據(jù)目標(biāo)評(píng)估其,然后修改分布(均值、sigma、協(xié)方差矩陣)等。每個(gè)x_for_eval點(diǎn)都是一個(gè)二進(jìn)制向量[1, 1, 1, 0, 0, 1, ...],用于從數(shù)據(jù)集中選擇特征。
請(qǐng)注意,這里使用的是 CMAwM() 優(yōu)化器(帶邊際的 CMA),而不是默認(rèn)的 CMA()。默認(rèn)優(yōu)化器在處理常規(guī)連續(xù)問題時(shí)效果很好,但這里的搜索空間是高維的,只允許兩個(gè)離散值(0 和 1)。默認(rèn)優(yōu)化器就會(huì)卡在這個(gè)空間里。CMAwM()擴(kuò)大了一點(diǎn)搜索空間(而它返回的解仍然是二進(jìn)制向量),這似乎足以解除它的阻礙。
這段簡(jiǎn)單的代碼確實(shí)有效,但還遠(yuǎn)未達(dá)到最佳狀態(tài)。下面是一個(gè)更復(fù)雜、更優(yōu)化的版本,它能更快地找到更好的解。
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']
下圖顯示了復(fù)雜、優(yōu)化的 CMA-ES 代碼搜索最佳解決方案的歷史。熱圖顯示了每一代中每種功能的流行程度(更亮 = 更流行)。可以看到,有些特征總是非常流行,而另一些則很快過時(shí),還有一些則是后來(lái)才 “發(fā)現(xiàn) ”的。根據(jù)這個(gè)問題的參數(shù),優(yōu)化器選擇的群體大小為 20 個(gè)點(diǎn)(個(gè)體),因此特征流行度是這 20 個(gè)點(diǎn)的平均值。
上下滑動(dòng)查看更多CMA-ES 優(yōu)化歷史
這些是經(jīng)過優(yōu)化的 CMA-ES 代碼的主要統(tǒng)計(jì)信息:
best objective: 33703.070530508514best generation: 921objective runs: 20000time to best: 48.326 sec
與 SFS 相比,它能找到更好(更小)的目標(biāo)值,調(diào)用目標(biāo)函數(shù)的次數(shù)更少(20k),所用時(shí)間也差不多。從所有指標(biāo)來(lái)看,這都是比 SFS 的凈收益。
同樣,特征選擇前的基準(zhǔn) BIC 為
baseline BIC: 34570.166173470934
順便提一句:在研究了傳統(tǒng)優(yōu)化算法(遺傳算法、模擬退火等)之后,CMA-ES 給我?guī)?lái)了驚喜。它幾乎沒有超參數(shù),計(jì)算量小,只需要少量個(gè)體(點(diǎn))來(lái)探索搜索空間,而且性能相當(dāng)不錯(cuò)。如果你需要解決優(yōu)化問題,值得把它放在你的工具箱里。
[1]Kaggle: https://www.kaggle.com/c/house-prices-advanced-regression-techniques/data
[2]mlxtend 庫(kù): https://rasbt.github.io/mlxtend/
[3]cmaes 庫(kù): https://github.com/CyberAgentAILab/cmaes
本文鏈接:http://www.www897cc.com/showinfo-26-101117-0.html協(xié)方差矩陣適應(yīng)進(jìn)化算法實(shí)現(xiàn)高效特征選擇
聲明:本網(wǎng)頁(yè)內(nèi)容旨在傳播知識(shí),若有侵權(quán)等問題請(qǐng)及時(shí)與本網(wǎng)聯(lián)系,我們將在第一時(shí)間刪除處理。郵件:2376512515@qq.com