高斯過程可以被認為是一種機器學習算法，它利用點與點之間同質性的度量作為核函數，以從輸入的訓練數據預測未知點的值。本文從理論推導和實現詳細地介紹了高斯過程，並在後面提供了用它來近似求未知函數最優解的方法。
我們回顧了高斯過程（GP）擬合數據所需的數學和代碼，最後得出一個常用應用的 demo——通過高斯過程搜索法快速實現函數最小化。下面的動圖演示了這種方法的動態過程，其中紅色的點是從紅色曲線採樣的樣本。使用這些樣本，我們試圖利用 GP 儘快找到曲線的最小值。
附錄包含（i）高斯回歸後驗推導； (ii)SKLearn 的 GP 實現；(iii) GP 分類器的快速回顧。
前言
高斯過程（GP）是處理以下一般問題的一個工具：函數 f(x) 在 n 個點採樣，並得到一組有噪聲 [1] 的函數度量值 {f(xi)=y_i ± σ_i,i=1,…,n}。那麼若給定這些可用的樣本，且 f hat 為一些候選函數，我們是否就能估計出 f =f hat 的概率?
為了逐步明確上述問題，我們首先應用貝葉斯法則，
上式左邊的數值是所求概率的簡寫，即給定樣本函數值 {y} 的條件下 f=f hat 的概率。在上式的右邊，分子中的第一項需要我們對測量過程中的誤差來源做一些假設，分子中的第二項是先驗概率，在這裡我們必須採用最合理的假設。例如，我們將在下面看到先驗概率有效地決定了 f 函數在給定平滑度的概率。
在 GP 方法中，右式中兩個分子都服從多元正態/高斯分布。模型可以選擇高斯的具體參數來達到良好的擬合度，但特徵的正態族假設對於解決數學問題是必不可少的。採用這種方法，我們可以通過分析寫出後驗概率，然後用在一些應用和計算中。例如，我們使用這種方法來獲得文中最開始處圖片上的曲線，即通過擬合 GP 後驗概率中隨機抽樣而得到曲線估計，在兩個收縮點處被固定為相等的測量值。後驗樣本對於可視化和求蒙特卡洛的平均值都很有用。
在本文中，我們做的工作有：
（i）回顧計算上述後驗概率所需的數學運算；
（ii）討論數值評估，並使用 GP 來擬合一些實例數據；
（iii）回顧擬合的 GP 如何快速最小化成本函數，例如機器學習中的交叉驗證分。
附錄包括高斯過程回歸推導，SKLearn 的 GP 實現和 GP 分類器的快速回顧。
我們的 GitHub 提供了簡單的高斯過程示例：https://github.com/EFavDB/gaussian_processes
注意：為了理解這篇文章中的數學細節，應該熟悉多元正態分布。但如果主要對應用感興趣，可以忽略這些細節。
後驗分析評估
為了計算 (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 +xixj。它不是穩定的，就是說如果加入一個常量的平移，結果就會改變。如果把 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)
關於 sklearn 實現的更多細節可以在這裡找到：http://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.GaussianProcessRegressor.html
附錄 C：GP 分類器
這裡，我們說明通常 GP 如何被用來擬合二進位類型數據，響應變量 y 可以取 0 或 1 的數據。GP 分類器的數學運算不像 GP 回歸那樣清楚。因為 0/1 響應不是高斯分布的，意味著後驗概率也不是。為了利用該程序，可以通過拉普拉斯（Laplace）近似正常地對後驗概率近似。
首先寫這樣一個公式，表示在 x 處看到一個給定的 y 值的概率。具體如下，
這個公式是 logistic 回歸的一個正常非線性泛化。此外，f 的先驗概率再次得到等式（3）。使用此式和（A8），我們可以得到 f 的後驗概率