從數學到實現，全面回顧高斯過程當中的函數最優化

時間 2019-12-04

原文原文鏈接

高斯過程能夠被認爲是一種機器學習算法，它利用點與點之間同質性的度量做爲核函數，以從輸入的訓練數據預測未知點的值。本文從理論推導和實現詳細地介紹了高斯過程，並在後面提供了用它來近似求未知函數最優解的方法。文章選自efavdb，做者： Jonathan Landy，機器之心編譯。html

咱們回顧了高斯過程（GP）擬合數據所需的數學和代碼，最後得出一個經常使用應用的 demo——經過高斯過程搜索法快速實現函數最小化。下面的動圖演示了這種方法的動態過程，其中紅色的點是從紅色曲線採樣的樣本。使用這些樣本，咱們試圖利用 GP 儘快找到曲線的最小值。python

附錄包含（i）高斯迴歸後驗推導； (ii)SKLearn 的 GP 實現；(iii) GP 分類器的快速回顧。
git

前言

高斯過程（GP）是處理如下通常問題的一個工具：函數 f(x) 在 n 個點採樣，並獲得一組有噪聲 [1] 的函數度量值 {f(xi)=y_i ± σ_i,i=1,…,n}。那麼若給定這些可用的樣本，且 f hat 爲一些候選函數，咱們是否就能估計出 f =f hat 的機率?github

爲了逐步明確上述問題，咱們首先應用貝葉斯法則，算法

上式左邊的數值是所求機率的簡寫，即給定樣本函數值 {y} 的條件下 f=f hat 的機率。在上式的右邊，分子中的第一項須要咱們對測量過程當中的偏差來源作一些假設，分子中的第二項是先驗機率，在這裏咱們必須採用最合理的假設。例如，咱們將在下面看到先驗機率有效地決定了 f 函數在給定平滑度的機率。bash

在 GP 方法中，右式中兩個分子都服從多元正態/高斯分佈。模型能夠選擇高斯的具體參數來達到良好的擬合度，但特徵的正態族假設對於解決數學問題是必不可少的。採用這種方法，咱們能夠經過分析寫出後驗機率，而後用在一些應用和計算中。例如，咱們使用這種方法來得到文中最開始處圖片上的曲線，即經過擬合 GP 後驗機率中隨機抽樣而獲得曲線估計，在兩個收縮點處被固定爲相等的測量值。後驗樣本對於可視化和求蒙特卡洛的平均值都頗有用。dom

在本文中，咱們作的工做有：機器學習

（i）回顧計算上述後驗機率所需的數學運算；
函數

（ii）討論數值評估，並使用 GP 來擬合一些實例數據；工具

（iii）回顧擬合的 GP 如何快速最小化成本函數，例如機器學習中的交叉驗證分。

附錄包括高斯過程迴歸推導，SKLearn 的 GP 實現和 GP 分類器的快速回顧。

咱們的 GitHub 提供了簡單的高斯過程示例：github.com/EFavDB/gaus…

注意：爲了理解這篇文章中的數學細節，應該熟悉多元正態分佈。但若是主要對應用感興趣，能夠忽略這些細節。

後驗分析評估

爲了計算 (1) 式左邊的值，咱們要先計算右邊的值。由於分母不依賴 f hat，咱們只須要考慮分子中的項。這意味着分母必須是全部候選函數共有的歸一化因子。在本節中，咱們先將分子兩項的估計公式寫出來，而後考慮後驗機率。

咱們要作的第一個假設是，假如實際函數是 f hat，那麼咱們的測量值 y 關於 f hat 是獨立的而且服從高斯分佈。這個假設意味着方程 (1) 右邊的第一項是：

上式中的 y_i 是咱們樣本點的實際測量值，σ_i 方是它們的方差不肯定度。

第二個假設是，假設先驗機率 p（f hat）的公式。咱們將注意力集中在一組數據點點 {x_i : i=1,…,N} 上，其中前 n 個點是已被抽樣的點，剩下的（N-n）個點是其餘位置的測試點，即咱們用來估計 f 函數聯合統計數據的點。接下來，簡單地假設這些點服從 f 的多元正態分佈，由協方差矩陣Σ來控制，由此獲得

這裏，咱們簡寫 f_i≡f（x_i）。請注意，咱們已經默認上面正態分佈的均值是零。這是爲了簡便起見：若是非零均值合適，那麼能夠與均值相加分析，或者從原始 f 中減去非零均值使新的分佈均值爲零。

Σ的特殊形式是在 GP 建模時最須要建模者觀察力和首創性的地方。對研究主題很是熟悉的研究者能夠構建很是好且很是複雜的先驗機率，而這種先驗一般是各項求和的形式，其中每一項都在所討論問題的數據基礎上加入了一些物理相關的貢獻。在這篇文章中，咱們假設一個簡單的公式，

注意，這個假設下，若是 x_i 和 x_j 很接近，那麼指數近似等於 1。這確保了附近的點高度相關，從而使全部高几率函數變得平滑。當兩測試點遠離時，式 (4) 中的衰減速率由長度參數 l 控制。若是 l 很大（小），曲線將在一個很長（短）的距離上平滑。咱們將在下一節中說明這些問題，並在下下節中解釋如何從已有的樣本數據中推斷合適的長度參數。

如今，若是咱們把式 (2) 和式 (3) 代入式 (1)，將獲得後驗機率 p(f1|{y}) 的表達式。這是一個指數函數，它的獨立變量是 f_i 函數中的二次項。也能夠說，與前驗機率同樣，後驗機率服從多變量正態分佈。經過簡單計算，就能夠獲得這個分佈均值和協方差的明確表達式：使用塊（符號），其中 0 對應於採樣點，1 對應於測試點，測試點的邊緣分佈是

y 是測量值的向量，長度爲 n,

方程 (5) 是高斯過程迴歸的一個主要結果——有了這個結果，咱們就能夠評估後驗機率了。注意，在採樣值 y 中全部點的均值是線性的，而且在測量值附近每一個點處的方差減少。若是你有興趣仔細推導這個結果，能夠參考咱們的附錄，在那裏有兩個推導。可是，在接下來的正文中，咱們僅簡單地探討這個公式的應用。

後驗機率的數值計算

在本節中，咱們將介紹式 (5) 的兩個典型應用：（i）在測試點 x 處評估後驗分佈的均值和標準差，（ii）從後驗機率中直接採樣函數 f_hat。前者能夠得到 f 函數在全部位置的置信區間，然後者能夠用來實現可視化和從後驗機率中得到通常的蒙特卡洛平均值。這兩個概念都在這篇文章的標題圖片中進行了說明：圖片中，咱們將 GP 擬合一個已有兩個測量點的一維函數。藍色陰影區域表示每一個位置函數值的一個σ置信區間，彩色曲線是後驗樣本。

咱們的 SimpleGP fitter 類的代碼能夠在 GitHub 上找到。咱們將在下文中解釋他是如何工做的，但若是對細節感興趣應該仔細查閱代碼。

區間

下面的代碼對咱們的 SimpleGP 類進行了初始化，定義了一些樣本位置、樣本值和不肯定性，而後評估了一組測試點後驗機率的均值和標準差。簡而言之，這個過程以下：經過擬合評估了出如今式（5）中的逆矩陣

，並保留結果供之後使用，這能夠避免在每一個測試點中從新評估這個逆矩陣。接下來，經過調用區間，針對每一個測試點再次評估式 (5)。

# Initialize fitter -- set covariance parameters
WIDTH_SCALE = 1.0
LENGTH_SCALE = 1.0
model = SimpleGP(WIDTH_SCALE, LENGTH_SCALE, noise=0) 

# Insert observed sample data here, fit 
sample_x = [-0.5, 2.5]
sample_y = [.5, 0]
sample_s = [0.01, 0.25]
model.fit(sample_x, sample_y, sample_s) 

# Get the mean and std at each point in x_test
test_x = np.arange(-5, 5, .05)
means, stds = model.interval(test_x)複製代碼

在以上代碼塊中，WIDTH_SCALE 和 LENGTH_SCALE 用來指定協方差矩陣式（4）。前者對應於等式中的σ，後者對應於 l。增長 WIDTH_SCALE 對應於未知函數大小不肯定度的增長，增長 LENGTH_SCALE 對應於增長咱們指望的函數平滑程度。下圖說明這些問題：這裏，經過設置 WIDTH_SCALE = LENGTH_SCALE = 1 得到藍色區間，經過設置 WIDTH_SCALE = 0.5 和 LENGTH_SCALE = 2 得到橙色區間。結果是橙色相對藍色後驗估計更加緊密平滑。在這兩幅圖中，實曲線表示後驗分佈均值，豎線表示一個σ置信區間。

後驗採樣

爲了從後驗機率中採樣實際函數，咱們將再次簡單地評估式 (5) 中的均值和協方差矩陣，此次是對咱們所求採樣函數的多個測試點進行。一旦咱們有了這些測試點後驗機率的均值和協方差矩陣，咱們可使用多元正態採樣的外部庫從 (5) 中抽取樣本——爲此，咱們使用了 python 中 numpy。下面代碼的最後一步執行這些步驟。

# Insert observed sample data here. 
sample_x = [-1.5, -0.5, 0.7, 1.4, 2.5, 3.0]
sample_y = [1, 2, 2, .5, 0, 0.5]
sample_s = [0.01, 0.25, 0.5, 0.01, 0.3, 0.01]

# Initialize fitter -- set covariance parameters
WIDTH_SCALE = 1.0
LENGTH_SCALE = 1.0
model = SimpleGP(WIDTH_SCALE, LENGTH_SCALE, noise=0)
model.fit(sample_x, sample_y, sample_s)

# Get the mean and std at each point in test_x
test_x = np.arange(-5, 5, .05)
means, stds = model.interval(test_x)

# Sample here
SAMPLES = 10
samples = model.sample(test_x, SAMPLES)
複製代碼

注意在第 2-4 行中，咱們添加了一些附加的函數樣本點（爲了好玩）。最後的區間和後驗樣本以下圖所示。注意到在採樣點附近，後驗結果很是理想。然而，在圖的左側，一旦咱們移動的距離≥1（即協方差矩陣 (4) 中的長度參數），後驗就趨近於先驗。

選擇協方差超參數

以前，咱們證實了咱們的協方差的長度參數顯着地影響後驗機率的區間形狀以及其中的樣本。適當設置這些參數是使用 GP 的一個廣泛難點。在這裏，咱們描述兩種方法，能夠巧妙地設置超參數，並給出一些採樣數據。

交叉驗證

交叉驗證是設置超參數的標準方法。這須要將可用的樣本數據分爲訓練集和驗證集。訓練集經過一系列超參數進行 GP 擬合，而後在已有的驗證集上評估模型的準確性。而後，經過選擇不一樣的超參數重複這個過程，選擇可使驗證集表現最優的一組。

邊緣似然最大化

一般狀況下，人們傾向於將 GP 應用於評估樣本有較高成本的狀況。這意味着人們一般在只有少數樣本可用的狀況下使用 GP。這種狀況下，隨着訓練點數量的增長，最優超參數能夠快速變化，意味着經過交叉驗證獲得的最優選擇可能遠不如訓練一個完整樣本集獲得的最優集合 [3]。

設置超參數的另外一種常見方法是使邊緣似然最大化。這就是說，咱們試圖最大化已關察樣本的可能性，於是在可用超參數的基礎上進行優化。具體來講，就是經過對未知的 f hat 進行積分來評估邊緣似然 [4]。

咱們能夠像附錄中評估後驗分佈那樣直接進行積分。但更快的方法是注意到 f 積分後，y 值服從以下的正態分佈

其中σ^2 * I_00 在式（6）中定義，由此得出，

上面兩項是矛盾的：第二項經過找出使指數最大的協方差矩陣來減少，最大化指數使數據過分擬合。然而，第一項與第二項相反，第一項是高斯積分的歸一化因子，它隨着衰減長度變短和對角線誤差下降而變大，至關於抑制複雜度的正則項。

實際中，爲了使式 (10) 最大，一般利用梯度的解析表達式和梯度降低法，這是 SKLearn 採起的方法。模型的一個優勢是可以優化 GP 的超參數。然而，式（10）不必定是凸的，一般存在多個局部最小值。爲了得到一個好的最小值，能夠嘗試在一些優秀的初始點上進行初始化。或者能夠在隨機點重複初始化梯度降低，最後選擇最優解。

函數最小搜索和機器學習

如今咱們將介紹 GP 的一個經常使用的應用：快速地搜索函數最小值。在這個問題中，咱們能夠迭代得到函數的噪聲樣本，從而儘快識別函數的全局最小值。梯度降低能夠應用於這種狀況，可是若是函數不具有凸性，一般須要重複採樣。爲了減小所需的步驟/樣本的數量，能夠嘗試應用更通常的探索式策略，即平衡「優化當前已知最小值的目標」與「尋找可能更小的新局部最小值的目標」。GP 後驗爲開發這樣的策略的提供了一個自然的起點。

GP 搜索法的想法是在 GP 後驗的基礎上得到一個得分函數。這個得分函數用來對搜索給定點的信息進行編碼，它能夠對探索（explore）和利用（exploit）造成一種權衡。一旦每一個點都進行評分，那麼具備最大（或最小，最合適的）分數的點將會被採樣。而後迭代重複該過程直到找到一個符合要求的解爲止。咱們將在下面討論四種可能的選擇，並給出一個例子。

高斯置信下界（GLCB）

GLCB 在每點的評分方式爲

這裏，μ和σ是函數在 x 處的均值和標準差的 GP 後驗估計值，κ是控制參數。請注意，第一項μ（x）鼓勵利用最可靠的局部最小值，並在它的周圍執行搜索。相似地，第二項κσ鼓勵在當前 GP 最不肯定真實函數值的點上進行探索。

改進的高斯機率（GPI）

若是目前爲止所看到的最小值是 y，則能夠利用該點處的真實函數值小於 y 的機率來給每一個點評分。也就是說，咱們能夠寫爲

高斯預期改進（EI）

上式常見的變形叫作預期改進，定義爲

這個得分函數傾向於鼓勵更多地去探索而不是改善機率，由於它更重視不肯定性。

機率最小值

要獲得的最終得分函數是問題中最小值的機率。得到這個分數的一個方法是進行屢次後驗採樣。對於每一個樣本，首先標記它的全局最小值，而後採起多數投票方法來決定接下來的樣本。

本文最開始處的動圖展現了一個實際的 GP 搜索，使用 skopt[5] 在 python 中執行。左邊的紅色曲線是正在尋找全局最小值的（隱藏）曲線 f。紅點是目前已經得到的樣本，綠色陰影曲線是每一個點的 GP 後驗置信區間，該置信區間會隨着更多樣本的得到而逐漸改進。右邊是經過在 GP 後驗基礎上分析獲得的每點的預期改進（EI）得分函數——該例中用於指導搜索的得分函數。該過程用五個隨機樣本進行初始化，而後進行引導搜索。請注意，隨着過程的演變，前幾個樣本集中利用已知的局部最小值。但通過幾回迭代後，繼續對這些位置採樣的收益減少，探索中間點的需求佔了上風——中間點是發現的實際全局最小值的點。

討論

在這篇文章中，咱們歸納了大部分 GP 的數學運算：獲得後驗機率所需的數學，如何進行後驗採樣，最後討論瞭如何實際應用後驗機率。

總的來講，GP 表明了一個能夠擬合任何函數的強大工具。實際中，使用這個工具的挑戰主要在於合適超參數的選擇，尋找合適的參數常常被困在局部最小，使擬合失效。不過，若是選擇得當，GP 的應用能夠提供一些有價值的性能提高。

附錄中討論了關於 GP 的其餘話題。若是對更多的細節感興趣，咱們能夠推薦 Rasmussen 和 Williams 的免費在線文本 [6]。

附錄 A：後驗推導

本附錄中，咱們提出後驗推導（5）的兩種方法。

方法 1

先平方，結合式（2）和式（3），簡單地計算得出

這裏，(1/σ^2) * I 在式（6）中被定義，但在樣本集外的全部行中都爲零。爲了獲得式（5），咱們必須統一爲正文中逆矩陣的分塊結構。

首先，咱們能夠得出

這裏咱們使用了分塊表示法。爲了計算上述的逆矩陣，咱們將利用分塊矩陣求逆公式，

矩陣（A2）中塊 C = 0、D = 1，這大大簡化了上述過程。代入後獲得

利用這個結果和（A1），咱們能夠獲得測試集的平均值

其中第二行的分子分母同時乘了 (1/σ^2) * I_00 的逆。相似地，測試集的協方差由（A3）的右下塊給出。獲得，

由（A4）和（A5）得出式（5）。

方法 2

在第二種方法中，咱們考慮一組測試點 f_1 和一組觀測樣本 f_0 的聯合分佈，咱們再次假設密度函數的均值爲零。那麼二者的聯合機率密度是

如今，咱們利用結果

式（5）的得出須要利用上述兩個表達式。主要的困難在於平方，相似於以前的推導，這能夠經過分塊矩陣求逆公式來完成。

附錄 B：SKLearn 實現和其餘內核

SKLearn 提供了包含 GaussianProcessRegressor 的類。使得能夠在任何維度上進行擬合和採樣——即它比咱們的最小類更通常，由於它能夠在多維上擬合特徵向量。另外，SKLearn 類的擬合方法嘗試爲一組給定的數據找到一組最佳的超參數。如上所述，這是經過邊緣似然的最大化來完成的。這裏，咱們提供一些關於這個類的基本註釋和能夠用來定義（3）中協方差矩陣Σ的內核函數，以及一段說明調用的簡單代碼。

預約義的核函數

徑向基函數（RBF）：這是默認值，至關於式（4）。RBF 由一個尺度參數 l 表徵，多維狀況下，能夠是一個容許各向異性相關長度的向量。
White kernel：The White Kernel 用於噪聲估計——文檔建議用於估計全局噪聲水平，但不是逐點。
Matern：這是一個廣義的指數衰減，其中指數是分離距離的冪律。特殊限制包括 RBF 和絕對距離指數衰減。
有理二次方程：（1+（d / l）2）α。
Exp-Sine-Squared：它容許模擬週期性函數。相似於 RBF，但其距離是實際距離的正弦。存在週期性參數和「方差」——高斯抑制（Gaussian suppression）的尺度。
點積核函數：格式爲 1 +xi⋅xj。它不是穩定的，就是說若是加入一個常量的平移，結果就會改變。若是把 N（0,1）的先驗值放在係數上，將獲得線性迴歸分析的結果。
核函數做爲對象：能夠支持核函數之間的二進制操做以建立更復雜的核函數，例如加法、乘法和指數（後者只是將初始核函數提高爲冪）。你能夠經過一些輔助函數來訪問核函數中的全部參數，例如 kernel.get_params().kernel.hyperparameters 是全部超參數的列表。

參數

n_restarts_optimizer：從新擬合的次數，用於探索多個局部最小值，默認值是零。
alpha：這個可選參數容許每一個測量都傳遞不肯定性。
normalize_y：用來表示咱們正在尋找的 y 的平均值不必定是零。

調用示例

下面的代碼進行了一次簡單擬合，結果是本文最開始展現的圖片。

from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C
from sklearn.gaussian_process import GaussianProcessRegressor
import numpy as np

# Build a model
kernel = C(1.0, (1e-3, 1e3)) * RBF(10, (0.5, 2))
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9)

# Some data
xobs = np.array([[1], [1.5], [-3]])
yobs = np.array([3, 0, 1])

# Fit the model to the data (optimize hyper parameters)
gp.fit(xobs, yobs)

# Plot points and predictions
x_set = np.arange(-6, 6, 0.1)
x_set = np.array([[i] for i in x_set])
means, sigmas = gp.predict(x_set, return_std=True)

plt.figure(figsize=(8, 5))
plt.errorbar(x_set, means, yerr=sigmas, alpha=0.5)
plt.plot(x_set, means, 'g', linewidth=4)
colors = ['g', 'r', 'b', 'k']
for c in colors:
y_set = gp.sample_y(x_set, random_state=np.random.randint(1000))
 plt.plot(x_set, y_set, c + '--', alpha=0.5)
複製代碼