盧偉冰長文解析K60至尊版 對Redmi有著里程碑式的意義
在今天的Redmi后性能時代戰略發布會結束之后，Redmi總經理盧偉冰又帶來了一篇長文，詳解了為什么 Redmi 要開啟后性能時代？為什么選擇和 MediaTek、Pixelworks 深度合作？以及后性
在建立模型時，特征選擇是一個重要環節，它指通過保留一部分特征子集來擬合模型，而舍棄其余特征。進行特征選擇有多重原因:
如果特征數量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.tebozhan.com/showinfo-26-101117-0.html協方差矩陣適應進化算法實現高效特征選擇
聲明：本網頁內容旨在傳播知識，若有侵權等問題請及時與本網聯系，我們將在第一時間刪除處理。郵件：2376512515@qq.com
Copyright ? 2016-2023 天津谷騏科技有限公司 版權所有 sitemap.xml
違法及侵權請聯系：2376512515@qq.com 津ICP備18001702號
津公網安備 12010102000574號