原文地址：https://borgwang.github.io/ml/2019/07/28/gaussian-processes.html
一元高斯分布
概率密度函數：$$p(x) = \frac{1}{\sigma\sqrt{2\pi}}\mathrm{exp}(-\frac{(x-\mu)^2}{2\sigma^2}) \tag{1}$$
其中\(\mu\)和\(\sigma\)分別表示均值和方差，這個概率密度函數曲線畫出來就是我們熟悉的鍾形曲線，均值和方差唯一地決定了曲線的形狀。
多元高斯分布
從一元高斯分布推廣到多元高斯分布，假設各維度之間相互獨立$$p(x_1, x_2, ..., x_n) = \prod_{i=1}^{n}p(x_i) \ =\frac{1}{(2\pi)^{\frac{n}{2}}\sigma_1\sigma_2...\sigma_n}\mathrm{exp}\left(-\frac{1}{2}\left [\frac{(x_1-\mu_1)^2}{\sigma_1^2} + \frac{(x_2-\mu_2)^2}{\sigma_2^2} + ... + \frac{(x_n-\mu_n)^2}{\sigma_n^2}\right]\right) \tag{2}$$其中\(\mu_1,\mu_2,\cdots\)和\(\sigma_1,\sigma_2,\cdots\)分別是第 1 維、第二維……的均值和方差。對上式向量和矩陣表示上式，令$$\boldsymbol{x - \mu}=[x_1-\mu_1, \ x_2-\mu_2,\ … \ x_n-\mu_n]^T \ K = \begin{bmatrix}
\sigma_1^2 & 0 & \cdots & 0\
0 & \sigma_2^2 & \cdots & 0\
\vdots & \vdots & \ddots & 0\
0 & 0 & 0 & \sigma_n^2
\end{bmatrix}$$
則：$$\sigma_1\sigma_2...\sigma_n = |K|^{\frac{1}{2}} \ \frac{(x_1-\mu_1)^2}{\sigma_1^2} + \frac{(x_2-\mu_2)^2}{\sigma_2^2} + ... + \frac{(x_n-\mu_n)^2}{\sigma_n^2}=(\boldsymbol{x-\mu})^TK^{-1}(\boldsymbol{x-\mu})$$
代入上式得到 $$p(\boldsymbol{x}) = (2\pi)^{-\frac{n}{2}}|K|^{-\frac{1}{2}}\mathrm{exp}\left( -\frac{1}{2}(\boldsymbol{x-\mu})^TK^{-1}(\boldsymbol{x-\mu}) \right) \tag{3}$$
其中\(\boldsymbol{\mu} \in \mathbb{R}^n\)是均值向量，\(K \in \mathbb{R}^{n \times n}\)為協方差矩陣，由於我們假設了各維度直接相互獨立，因此\(K\)是一個對角矩陣。在各維度變量相關時，上式的形式仍然一致，但此時協方差矩陣\(K\)不再是對角矩陣，只具備半正定和對稱的性質。上式通常也簡寫為$$x \sim \mathcal{N}(\boldsymbol{\mu}, K)$$
無限元高斯分布
用一個例子來展示這個擴展的過程:MLSS 2012: J. Cunningham - Gaussian Processes for Machine Learning
假設我們在周一到周四每天的 7:00 測試了 4 次心率，如圖中 4 個點，可能的高斯分布如圖所示（高瘦的曲線）。這是一個一元高斯分布，只有每天 7: 00 心率這個維度。
現在考慮不僅在每天的 7: 00 測心率，在 8:00 時也進行測量，這個時候變成兩個維度（二元高斯分布），如圖所示：
如果我們在每天的無限個時間點都進行測量，則變成了下圖的情況。注意下圖中把測量時間作為橫軸，則每個顏色的一條線代表一個（無限個時間點的測量）無限維的采樣。當對無限維進行采樣得到無限多個點時，其實可以理解為對函數進行采樣。
當從函數的視角去看待采樣，理解了每次采樣無限維相當於采樣一個函數之后，原本的概率密度函數不再是點的分布 ，而變成了函數的分布。這個無限元高斯分布即稱為高斯過程。
高斯過程正式地定義為：對於所有\(\boldsymbol{x} = [x_1, x_2, \cdots, x_n],f(\boldsymbol{x})=[f(x_1), f(x_2), \cdots, f(x_n)]\)都服從多元高斯分布，則稱\(f\)是一個高斯過程，表示為$$f(\boldsymbol{x}) \sim \mathcal{N}(\boldsymbol{\mu}(\boldsymbol{x}), \kappa(\boldsymbol{x},\boldsymbol{x})) \tag{4}$$這里\(\boldsymbol{\mu}(\boldsymbol{x}): \mathbb{R^{n}} \rightarrow \mathbb{R^{n}}\)表示均值函數（Mean function），返回各個維度的均值；\(\kappa(\boldsymbol{x},\boldsymbol{x}) : \mathbb{R^{n}} \times \mathbb{R^{n}} \rightarrow \mathbb{R^{n\times n}}\)為協方差函數 Covariance Function（也叫核函數 Kernel Function）返回各個維度之間的協方差矩陣。一個高斯過程為一個均值函數和協方差函數唯一地定義，並且一個高斯過程的有限維度的子集都服從一個多元高斯分布（為了方便理解，可以想象二元高斯分布兩個維度各自都服從一個高斯分布）。
核函數（協方差函數）
核函數是一個高斯過程的核心，核函數決定了一個高斯過程的性質。核函數在高斯過程中起的作用是生成一個協方差矩陣（相關系數矩陣），衡量任意兩個點之間的“距離”。最常用的一個核函數為高斯核函數，也成為徑向基函數 RBF。其基本形式如下。其中\(\sigma\)和\(l\)是高斯核的超參數。$$K(x_i,x_j)=\sigma^2\mathrm{exp}\left( -\frac{\left |x_i-x_j\right |_2^2}{2l^2}\right)$$
高斯核函數的 python 實現如下：
import numpy as np
def gaussian_kernel(x1, x2, l=1.0, sigma_f=1.0):
"""Easy to understand but inefficient."""
m, n = x1.shape[0], x2.shape[0]
dist_matrix = np.zeros((m, n), dtype=float)
for i in range(m):
for j in range(n):
dist_matrix[i][j] = np.sum((x1[i] - x2[j]) ** 2)
return sigma_f ** 2 * np.exp(- 0.5 / l ** 2 * dist_matrix)
def gaussian_kernel_vectorization(x1, x2, l=1.0, sigma_f=1.0):
"""More efficient approach."""
dist_matrix = np.sum(x1**2, 1).reshape(-1, 1) + np.sum(x2**2, 1) - 2 * np.dot(x1, x2.T)
return sigma_f ** 2 * np.exp(-0.5 / l ** 2 * dist_matrix)
x = np.array([700, 800, 1029]).reshape(-1, 1)
print(gaussian_kernel_vectorization(x, x, l=500, sigma=10))
輸出的協方差矩陣為
[[100. 98.02 80.53]
[ 98.02 100. 90.04]
[ 80.53 90.04 100. ]]
高斯過程更新
下圖是高斯過程的可視化，其中藍線是高斯過程的均值，淺藍色區域 95% 置信區間（由協方差矩陣的對角線得到），每條虛線代表一個函數采樣（這里用了 100 維模擬連續無限維）。左上角第一幅圖是高斯過程的先驗（這里用了零均值作為先驗），后面幾幅圖展示了當觀測到新的數據點的時候，高斯過程如何更新自身的均值函數和協方差函數。
將高斯過程的先驗表示為\(f(\boldsymbol{x}) \sim \mathcal{N}(\boldsymbol{\mu}_{f}, K_{ff})\)，如果現在我們觀測到一些數據\((\boldsymbol{x'}, \boldsymbol{y})\)，並且假設\(\boldsymbol{y}\)與\(f(\boldsymbol{x})\)服從聯合高斯分布$$\begin{bmatrix}
f(\boldsymbol{x})\
\boldsymbol{y}\
\end{bmatrix} \sim \mathcal{N} \left(
\begin{bmatrix}
\boldsymbol{\mu_f}\
\boldsymbol{\mu_y}\
\end{bmatrix},
\begin{bmatrix}
K_{ff} & K_{fy}\
K_{fy}^T & K_{yy}\
\end{bmatrix}
\right)$$
其中\(K_{ff} = \kappa(\boldsymbol{x}, \boldsymbol{x})\),\(K_{fy}=\kappa(\boldsymbol{x}, \boldsymbol{x'})\),\(K_{yy} = \kappa(\boldsymbol{x'}, \boldsymbol{x'})\)，則有$$f|\boldsymbol{y} \sim \mathcal{N}(K_{fy}K_{yy}^{-1}\boldsymbol{y}+\boldsymbol{\mu_f},K_{ff}-K_{fy}K_{yy}^{-1}K_{fy}^T) \tag{5}$$
公式（5）表明了給定數據\((\boldsymbol{x'}, \boldsymbol{y})\)之后函數的分布\(f\)仍然是一個高斯過程。
上式其實就是高斯過程回歸的基本公式，首先有一個高斯過程先驗分布，觀測到一些數據（機器學習中的訓練數據），基於先驗和一定的假設（聯合高斯分布）計算得到高斯過程后驗分布的均值和協方差。
簡單高斯過程回歸實現
考慮代碼實現一個高斯過程回歸，API 接口風格采用 sciki-learn fit-predict 風格。由於高斯過程回歸是一種非參數化的方法，每次的 inference 都需要利用所有的訓練數據進行計算得到結果，因此並沒有一個顯式的訓練模型參數的過程，所以 fit 方法只需要將訓練數據記錄下來，實際的 inference 在 predict 方法中進行。Python 代碼如下
from scipy.optimize import minimize
class GPR:
def __init__(self, optimize=True):
self.is_fit = False
self.train_X, self.train_y = None, None
self.params = {"l": 0.5, "sigma_f": 0.2}
self.optimize = optimize
def fit(self, X, y):
# store train data
self.train_X = np.asarray(X)
self.train_y = np.asarray(y)
self.is_fit = True
def predict(self, X):
if not self.is_fit:
print("GPR Model not fit yet.")
return
X = np.asarray(X)
Kff = self.kernel(X, X) # (N,N)
Kyy = self.kernel(self.train_X, self.train_X) # (k,k)
Kfy = self.kernel(X, self.train_X) # (N,k)
Kyy_inv = np.linalg.inv(Kyy + 1e-8 * np.eye(len(self.train_X))) # (k,k)
mu = Kfy.dot(Kyy_inv).dot(self.train_y)
cov = self.kernel(X, X) - Kfy.dot(Kyy_inv).dot(Kfy.T)
return mu, cov
def kernel(self, x1, x2):
dist_matrix = np.sum(x1**2, 1).reshape(-1, 1) + np.sum(x2**2, 1) - 2 * np.dot(x1, x2.T)
return self.params["sigma_f"] ** 2 * np.exp(-0.5 / self.params["l"] ** 2 * dist_matrix)
def y(x, noise_sigma=0.0):
x = np.asarray(x)
y = np.cos(x) + np.random.normal(0, noise_sigma, size=x.shape)
return y.tolist()
train_X = np.array([3, 1, 4, 5, 9]).reshape(-1, 1)
train_y = y(train_X, noise_sigma=1e-4)
test_X = np.arange(0, 10, 0.1).reshape(-1, 1)
gpr = GPR()
gpr.fit(train_X, train_y)
mu, cov = gpr.predict(test_X)
test_y = mu.ravel()
uncertainty = 1.96 * np.sqrt(np.diag(cov))
plt.figure()
plt.title("l=%.2f sigma_f=%.2f" % (gpr.params["l"], gpr.params["sigma_f"]))
plt.fill_between(test_X.ravel(), test_y + uncertainty, test_y - uncertainty, alpha=0.1)
plt.plot(test_X, test_y, label="predict")
plt.scatter(train_X, train_y, label="train", c="red", marker="x")
plt.legend()
結果如下圖，紅點是訓練數據，藍線是預測值，淺藍色區域是 95% 置信區間。真實的函數是一個 cosine 函數，可以看到在訓練數據點較為密集的地方，模型預測的不確定性較低，而在訓練數據點比較稀疏的區域，模型預測不確定性較高。
超參數優化
上文提到高斯過程是一種非參數模型，沒有訓練模型參數的過程，一旦核函數、訓練數據給定，則模型就被唯一地確定下來。但是核函數本身是有參數的，比如高斯核的參數\(\sigma\)和\(l\)，我們稱為這種參數為模型的超參數（類似於 k-NN 模型中 k 的取值）。
核函數本質上決定了樣本點相似性的度量方法，進行影響到了整個函數的概率分布的形狀。上面的高斯過程回歸的例子中使用了\(\sigma=0.2\)和\(l=0.5\)的超參數，我們可以選取不同的超參數看看回歸出來的效果。
從上圖可以看出，\(l\)越大函數更加平滑，同時訓練數據點之間的預測方差更小，反之\(l\)越小則函數傾向於更加“曲折”，訓練數據點之間的預測方差更大；\(\sigma\)則直接控制方差大小，\(\sigma\)越大方差越大，反之亦然。
如何選擇最優的核函數參數\(\sigma\)和\(l\)呢？答案最大化在這兩個超參數下\(y\)出現的概率，通過最大化邊緣對數似然（Marginal Log-likelihood）來找到最優的參數，邊緣對數似然表示為$$\mathrm{log}\ p(\boldsymbol{y}|\sigma, l) = \mathrm{log} \ \mathcal{N}(\boldsymbol{0}, K_{yy}(\sigma, l)) = -\frac{1}{2}\boldsymbol{y}^T K_{yy}^{-1}\boldsymbol{y} - \frac{1}{2}\mathrm{log}\ |K_{yy}| - \frac{N}{2}\mathrm{log} \ (2\pi) \tag{6}$$
具體的實現中，我們在 fit 方法中增加超參數優化這部分的代碼，最小化負邊緣對數似然。
from scipy.optimize import minimize
class GPR:
def __init__(self, optimize=True):
self.is_fit = False
self.train_X, self.train_y = None, None
self.params = {"l": 0.5, "sigma_f": 0.2}
self.optimize = optimize
def fit(self, X, y):
# store train data
self.train_X = np.asarray(X)
self.train_y = np.asarray(y)
# hyper parameters optimization
def negative_log_likelihood_loss(params):
self.params["l"], self.params["sigma_f"] = params[0], params[1]
Kyy = self.kernel(self.train_X, self.train_X) + 1e-8 * np.eye(len(self.train_X))
return 0.5 * self.train_y.T.dot(np.linalg.inv(Kyy)).dot(self.train_y) + 0.5 * np.linalg.slogdet(Kyy)[1] + 0.5 * len(self.train_X) * np.log(2 * np.pi)
if self.optimize:
res = minimize(negative_log_likelihood_loss, [self.params["l"], self.params["sigma_f"]],
bounds=((1e-4, 1e4), (1e-4, 1e4)),
method='L-BFGS-B')
self.params["l"], self.params["sigma_f"] = res.x[0], res.x[1]
self.is_fit = True
將訓練、優化得到的超參數、預測結果可視化如下圖，可以看到最優的\(l=1.2\)，\(\sigma_f=0.8\)
這里用 scikit-learn 的 GaussianProcessRegressor 接口進行對比
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, RBF
# fit GPR
kernel = ConstantKernel(constant_value=0.2, constant_value_bounds=(1e-4, 1e4)) * RBF(length_scale=0.5, length_scale_bounds=(1e-4, 1e4))
gpr = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=2)
gpr.fit(train_X, train_y)
mu, cov = gpr.predict(test_X, return_cov=True)
test_y = mu.ravel()
uncertainty = 1.96 * np.sqrt(np.diag(cov))
# plotting
plt.figure()
plt.title("l=%.1f sigma_f=%.1f" % (gpr.kernel_.k2.length_scale, gpr.kernel_.k1.constant_value))
plt.fill_between(test_X.ravel(), test_y + uncertainty, test_y - uncertainty, alpha=0.1)
plt.plot(test_X, test_y, label="predict")
plt.scatter(train_X, train_y, label="train", c="red", marker="x")
plt.legend()
得到結果為\(l=1.2\)，\(\sigma_f=0.6\)，這個與我們實現的優化得到的超參數有一點點不同，可能是實現的細節有所不同導致。
多維輸入
我們上面討論的訓練數據都是一維的，高斯過程直接可以擴展於多維輸入的情況，直接將輸入維度增加即可。
def y_2d(x, noise_sigma=0.0):
x = np.asarray(x)
y = np.sin(0.5 * np.linalg.norm(x, axis=1))
y += np.random.normal(0, noise_sigma, size=y.shape)
return y
train_X = np.random.uniform(-4, 4, (100, 2)).tolist()
train_y = y_2d(train_X, noise_sigma=1e-4)
test_d1 = np.arange(-5, 5, 0.2)
test_d2 = np.arange(-5, 5, 0.2)
test_d1, test_d2 = np.meshgrid(test_d1, test_d2)
test_X = [[d1, d2] for d1, d2 in zip(test_d1.ravel(), test_d2.ravel())]
gpr = GPR(optimize=True)
gpr.fit(train_X, train_y)
mu, cov = gpr.predict(test_X)
z = mu.reshape(test_d1.shape)
fig = plt.figure(figsize=(7, 5))
ax = Axes3D(fig)
ax.plot_surface(test_d1, test_d2, z, cmap=cm.coolwarm, linewidth=0, alpha=0.2, antialiased=False)
ax.scatter(np.asarray(train_X)[:,0], np.asarray(train_X)[:,1], train_y, c=train_y, cmap=cm.coolwarm)
ax.contourf(test_d1, test_d2, z, zdir='z', offset=0, cmap=cm.coolwarm, alpha=0.6)
ax.set_title("l=%.2f sigma_f=%.2f" % (gpr.params["l"], gpr.params["sigma_f"]))
下面是一個二維輸入數據的告你過程回歸，左圖是沒有經過超參優化的擬合效果，右圖是經過超參優化的擬合效果。
