強大而精致的機器學習調參方法：貝葉斯優化

一、簡介

貝葉斯優化用於機器學習調參由J. Snoek(2012)提出，主要思想是，給定優化的目標函數(廣義的函數，只需指定輸入和輸出即可，無需知道內部結構以及數學性質)，通過不斷地添加樣本點來更新目標函數的后驗分布(高斯過程,直到后驗分布基本貼合於真實分布。簡單的說，就是考慮了上一次參數的信息**，從而更好的調整當前的參數。

他與常規的網格搜索或者隨機搜索的區別是：

• 貝葉斯調參采用高斯過程，考慮之前的參數信息，不斷地更新先驗；網格搜索未考慮之前的參數信息
• 貝葉斯調參迭代次數少，速度快；網格搜索速度慢,參數多時易導致維度爆炸
• 貝葉斯調參針對非凸問題依然穩健；網格搜索針對非凸問題易得到局部優最

二、理論

介紹貝葉斯優化調參，必須要從兩個部分講起：

• 高斯過程，用以擬合優化目標函數
• 貝葉斯優化，包括了“開采”和“勘探”，用以花最少的代價找到最優值

2.1 高斯過程

高斯過程可以用於非線性回歸、非線性分類、參數尋優等等。以往的建模需要對 \(p(y|X)\) 建模，當用於預測時，則是

\[p(y_{N+1} | X_{N+1}) \]

而高斯過程則， 還考慮了 \(y_N\) 和 \(y_{N+1}\) 之間的關系，即：

\[p(y_{N+1} | X_{N+1}, y_{N}) \]

高斯過程通過假設 \(Y\) 值服從聯合正態分布，來考慮 \(y_N\) 和 \(y_{N+1}\) 之間的關系，因此需要給定參數包括：均值向量和協方差矩陣，即：

\[\begin{bmatrix} y_1 \\ y_2 \\ ... \\ y_n \\ \end{bmatrix} \sim N( \mathbf{0}, \begin{bmatrix} k(x_1, x_1) , k(x_1, x_2), ..., k(x_1, x_n) \\ k(x_2, x_1) , k(x_2, x_2), ..., k(x_2, x_n) \\ ... \\ k(x_n, x_1) , k(x_n, x_2), ..., k(x_n, x_n) \end{bmatrix} ) \]

其中協方差矩陣又叫做 核矩陣, 記為 \(\mathbf{K}\) ，僅和特征 \(x\) 有關，和 \(y\) 無關。

高斯過程的思想是： 假設 \(Y\) 服從高維正態分布（先驗），而根據訓練集可以得到最優的核矩陣 ，從而得到后驗以估計測試集 \(Y*\)

我們有后驗：

\[p(y_*| \mathbf{y} \sim N(K_* K^{-1} \mathbf{y}, ~ K_{**} - K_* K^{-1} K_*^T) \]

其中，\(K_*\)為訓練集的核向量，有如下關系：

\[\begin{bmatrix} \mathbf{y} \\ y_* \end{bmatrix} \sim N(\mathbf{0}, \begin{bmatrix} K, K_*^T \\ K_*, K_{**} \\ \end{bmatrix}) \]

可以發現，在后驗公式中，只有均值和訓練集 \(Y\) 有關，方差則僅僅和核矩陣，也就是訓練集和測試集的 \(X\) 有關，與訓練集 \(Y\) 無關

高斯過程的估計（訓練）方法

假設使用平方指數核(Squared Exponential Kernel)，那么有：

\[k(x_1, x_2) = \sigma^2_f exp(\frac{-(x_1 - x_2)^2}{2 l^2}) \]

那么所需要的確定的超參數 \(\theta = [\sigma^2_f, l]\) ，由於 \(Y\) 服從多維正態分布，因此似然函數為：

\[L = log p(y| x, \theta) = - \frac{1}{2} log|\mathbf{K}| - \frac{1}{2} (y - \mu)^T \mathbf{K}^{-1} (y - \mu) - n*log(2\pi)/2 \]

由於 \(K\) 是由 \(\theta\) 決定的，所以通過梯度下降即可求出超參數 \(\theta\)，而根據核矩陣的計算方式也可以進行預測。

上圖是一張高斯分布擬合函數的示意圖，可以看到，它只需要九個點，就可以大致擬合出整個函數形狀（圖片來自：https://github.com/fmfn/BayesianOptimization）

2.2 貝葉斯優化理論

貝葉斯優化是一種逼近思想，當計算非常復雜、迭代次數較高時能起到很好的效果，多用於超參數確定

基本思想

是基於數據使用貝葉斯定理估計目標函數的后驗分布，然后再根據分布選擇下一個采樣的超參數組合。它充分利用了前一個采樣點的信息，其優化的工作方式是通過對目標函數形狀的學習，並找到使結果向全局最大提升的參數

高斯過程 用於在貝葉斯優化中對目標函數建模，得到其后驗分布

通過高斯過程建模之后，我們嘗試抽樣進行樣本計算，而貝葉斯優化很容易在局部最優解上不斷采樣，這就涉及到了開發和探索之間的權衡。

• 開發 (exploitation)： 根據后驗分布，在最可能出現全局最優解的區域進行采樣, 開發高意味着均值高
• 探索 (exploration): 在還未取樣的區域獲取采樣點， 探索高意味着方差高

而如何高效的采樣，即開發和探索，我們需要用到 Acquisition Function, 它是用來尋找下一個 x 的函數。

Acquistion Function

一般形式的Acquisition Funtion是關於x的函數，映射到實數空間R，表示改點的目標函數值能夠比當前最優值大多少的概率，目前主要有以下幾種主流的Acquisition Function

POI(probability of improvement)

\[POI(X) = P(f(X) \ge f(X^+) + \xi) = \Phi(\frac{\mu(x) - f(X^+) - \xi}{\sigma(x)}) \]

其中， \(f(X)\) 為X的目標函數值， \(f(X^+)\) 為 到目前為止 最優的X的目標函數值， \(\mu(x), \sigma(x)\) 分別是高斯過程所得到的目標函數的均值和方差，即 \(f(X)\) 的后驗分布。 \(\xi\) 為trade-off系數,如果沒有該系數，POI函數會傾向於取在 \(X^+\) 周圍的點，即傾向於exploit而不是explore，因此加入該項進行權衡。

而我們要做的，就是嘗試新的X，使得 \(POI(X)\) 最大，則采取該\(X\) （因為\(f(X)\)的計算代價非常大），通常我們使用 蒙特卡洛模擬 的方法進行。

詳細情況見下圖（圖片來自 Ref[5])

Expected Improvement

POI是一個概率函數，因此只考慮了f(x) 比 \(f(x^+)\) 大的概率，而EI則是一個期望函數，因此考慮了 f(x) 比 \(f(x^+)\) 大多少。我們通過下式獲取x

\[x = argmax_x \ \ E(\max\{0, f_{t+1}(x) - f(X^+)\}| D_t) \]

其中 \(D_t\) 為前t個樣本，在正態分布的假定下，最終得到:

\[EI(x) = \begin{cases} (\mu(x) - f(x^+)) \Phi(Z) + \sigma(x) \phi(Z), if \ \sigma(x) > 0 \\ 0, if \ \sigma(x) = 0 \end{cases} \]

其中 \(Z= \frac{\mu(x) - f(x^+)}{\sigma(x)}\)

Confidence bound criteria

\[LCB(x) = \mu(x) - \kappa \sigma(x) \]

\[UCB(x) = \mu(x) + \kappa \sigma(x) \]

2.3 缺點和不足

• 高斯過程核矩陣不好選

三、例子

目前可以做貝葉斯優化的包非常多,光是python就有:

• BayesianOptimization
• bayesopt
• skopt
• ...

本文使用BayesianOptimization為例，利用sklearn的隨機森林模型進行分類

安裝

pip install bayesian-optimization

前期准備

from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.cross_validation import cross_val_score
from bayes_opt import BayesianOptimization

# 產生隨機分類數據集，10個特征， 2個類別
x, y = make_classification(n_samples=1000,n_features=10,n_classes=2)

我們先看看不調參的結果：

rf = RandomForestClassifier()
print(np.mean(cross_val_score(rf, x, y, cv=20, scoring='roc_auc')))

>>> 0.965162

可以看到，不調參的話模型20此交叉驗證AUC均值是0.965162，算是一個不錯的模型，那么如果用bayes調參結果會怎么樣呢

bayes調參初探

我們先定義一個目標函數，里面放入我們希望優化的函數。比如此時，函數輸入為隨機森林的所有參數，輸出為模型交叉驗證5次的AUC均值，作為我們的目標函數。因為bayes_opt庫只支持最大值，所以最后的輸出如果是越小越好，那么需要在前面加上負號，以轉為最大值。由於bayes優化只能優化連續超參數，因此要加上int()轉為離散超參數。

def rf_cv(n_estimators, min_samples_split, max_features, max_depth):
    val = cross_val_score(
        RandomForestClassifier(n_estimators=int(n_estimators),
            min_samples_split=int(min_samples_split),
            max_features=min(max_features, 0.999), # float
            max_depth=int(max_depth),
            random_state=2
        ),
        x, y, scoring='roc_auc', cv=5
    ).mean()
    return val

然后我們就可以實例化一個bayes優化對象了：

 rf_bo = BayesianOptimization(
        rf_cv,
        {'n_estimators': (10, 250),
        'min_samples_split': (2, 25),
        'max_features': (0.1, 0.999),
        'max_depth': (5, 15)}
    )

里面的第一個參數是我們的優化目標函數，第二個參數是我們所需要輸入的超參數名稱，以及其范圍。超參數名稱必須和目標函數的輸入名稱一一對應。

完成上面兩步之后，我們就可以運行bayes優化了！

rf_bo.maximize()

完成的時候會不斷地輸出結果，如下圖所示：

等到程序結束，我們可以查看當前最優的參數和結果：

rf_bo.res['max']

>>> {'max_params': {'max_depth': 5.819908283575526,
  'max_features': 0.4951745603509127,
  'min_samples_split': 2.3110014720414958,
  'n_estimators': 249.73529231990733},
 'max_val': 0.9774079407940794}

bayes調參進階

上面bayes算法得到的參數並不一定最優，當然我們會遇到一種情況，就是我們已經知道有一組或是幾組參數是非常好的了，我們想知道其附近有沒有更好的。這個操作相當於上文bayes優化中的Explore操作，而bayes_opt庫給了我們實現此方法的函數：

rf_bo.explore(
    {'n_estimators': [10, 100, 200],
        'min_samples_split': [2, 10, 20],
        'max_features': [0.1, 0.5, 0.9],
        'max_depth': [5, 10, 15]
    }
)

這里我們添加了三組較優的超參數，讓其在該參數基礎上進行explore，可能會得到更好的結果。

同時，我們還可以修改高斯過程的參數，高斯過程主要參數是核函數(kernel)，還有其他參數可以參考sklearn.gaussianprocess

gp_param={'kernel':None}
rf_bo.maximize(**gp_param)

最終我們的到參數如下：

{'max_params': {'max_depth': 5.819908283575526,
  'max_features': 0.4951745603509127,
  'min_samples_split': 2.3110014720414958,
  'n_estimators': 249.73529231990733},
 'max_val': 0.9774079407940794}

運行交叉驗證測試一下：

rf = RandomForestClassifier(max_depth=6, max_features=0.39517, min_samples_split=2, n_estimators=250)
np.mean(cross_val_score(rf, x, y, cv=20, scoring='roc_auc'))
>>> 0.9754953

得到最終結果是0.9755，比之前的0.9652提高了約0.01，做過kaggle的朋友都懂，這在后期已經是非常大的提高了！到后面想提高0.001都極其困難，因此bayes優化真的非常強大！

結束！

Reference

• [1] J. Snoek, H. Larochelle, and R. P. Adams, “Practical bayesianoptimization of machine learning algorithms,” in Advances in neural information processing systems, 2012, pp. 2951–2959.
• [2] 高斯過程：http://www.gaussianprocess.org/gpml/
• [3] 高斯過程：https://www.zhihu.com/question/46631426?sort=created
• [4] 高斯過程：http://www.360doc.com/content/17/0810/05/43535834_678049865.shtml
• [5] Brochu E, Cora V M, De Freitas N. A tutorial on Bayesian optimization of expensive cost functions, with application to active user modeling and hierarchical reinforcement learning[J]. arXiv preprint arXiv:1012.2599, 2010.