來源：DeepHub IMBA

本文約5000字，建議閱讀5分鐘

特征選擇旨在找到最具信息的特征，以防止過擬合，提高泛化能力，降低計算成本。

特征選擇是指從原始特征集中選擇一部分特征，以提高模型性能、減少計算開銷或改善模型的解釋性。特征選擇的目標是找到對目標變量預測最具信息量的特征，同時減少不必要的特征。這有助於防止過擬合、提高模型的泛化能力，並且可以減少訓練和推理的計算成本。

如果特征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:   44.850 sec
 GA:     157.604 sec
 CMA-ES: 46.912 sec

目標函數的最佳值-越少越好：

 SFS:   44.850 sec
 GA:     157.604 sec
 CMA-ES: 46.912 sec

CMA-ES找到了最佳目標函數。它的運行時間與SFS相當。它隻調用目標函數20k次，是所有方法中調用次數最少的。別忘了，他還是單線程的。

GA能夠在目標函數上超過SFS，但它是最慢的。它調用目標函數的次數比其他方法多一個數量級。這是因為我們可以認為他是一個半隨機的過程，因為遺傳突變這個東西沒有算法可解釋。

SFS速度很快(可以在所有CPU內核上運行)，但性能一般。但它是目前最簡單的算法。

如果你隻是想用一個簡單的算法快速估計出最佳的特征集，那麼SFS還不錯。如果你想要絕對最好的客觀價值，CMA-ES似乎是首選，並且它也不慢。

最後本文的代碼：

https://github.com/FlorinAndrei/fast_feature_selection

作者：Florin Andrei