跳到主要内容

高斯过程

高斯过程（Gaussian Process, GP）是一种强大的非参数监督学习方法，既可以用于回归问题，也可以用于概率分类。与传统的参数化模型不同，高斯过程不需要假设数据服从特定的函数形式，而是通过核函数定义函数空间上的概率分布。

高斯过程的核心优势在于它能够提供预测的不确定性估计——不仅告诉你预测值是多少，还能告诉你对预测有多大的信心。这种特性使其在主动学习、贝叶斯优化、传感器融合等场景中极具价值。

核心概念​

什么是高斯过程？​

从数学角度看，高斯过程是随机变量的集合，其中任意有限个子集的联合分布都是多元高斯分布。我们可以把高斯过程理解为"函数上的分布"——它不是对参数建模，而是直接对函数建模。

形式化定义：一个高斯过程由均值函数 m(x)m(x) 和协方差函数（核函数）k(x,x′)k(x, x') 完全定义：

f(x)∼GP(m(x),k(x,x′))f(x) \sim \mathcal{GP}(m(x), k(x, x'))

对于任意输入点集合 {x1,...,xn}\{x_1, ..., x_n\}，函数值 {f(x1),...,f(xn)}\{f(x_1), ..., f(x_n)\} 服从多元高斯分布：

[f(x1)⋮f(xn)]∼N([m(x1)⋮m(xn)],[k(x1,x1)⋯k(x1,xn)⋮⋱⋮k(xn,x1)⋯k(xn,xn)])\begin{bmatrix} f(x_1) \\ \vdots \\ f(x_n) \end{bmatrix} \sim \mathcal{N}\left(\begin{bmatrix} m(x_1) \\ \vdots \\ m(x_n) \end{bmatrix}, \begin{bmatrix} k(x_1, x_1) & \cdots & k(x_1, x_n) \\ \vdots & \ddots & \vdots \\ k(x_n, x_1) & \cdots & k(x_n, x_n) \end{bmatrix}\right)

直观理解：想象我们有无限多个可能的函数，高斯过程为每个函数分配一个概率。当我们观测到数据点后，不符合数据的函数概率降低，符合数据的函数概率增加。这就是高斯过程的"学习"过程。

高斯过程的优缺点​

优点：

  • 预测插值观测值：对于常规核函数，预测会穿过训练数据点
  • 概率输出：提供预测的均值和方差，可以量化不确定性
  • 灵活性：可以通过选择不同的核函数来编码先验知识
  • 非参数：不需要指定函数的具体形式

缺点：

  • 计算复杂度高：训练复杂度为 O(n3)O(n^3)，预测复杂度为 O(n2)O(n^2)
  • 高维效果下降：当特征数量超过几十个时，效率明显降低
  • 核函数选择困难：需要根据问题选择合适的核函数

高斯过程回归​

GaussianProcessRegressor 实现了高斯过程回归。它会结合先验分布和训练数据的似然函数，给出预测的均值和标准差。

基本用法​

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel
import numpy as np
import matplotlib.pyplot as plt

# 生成示例数据
np.random.seed(42)
X = np.linspace(0, 10, 20).reshape(-1, 1)
y = np.sin(X).ravel() + np.random.normal(0, 0.1, X.shape[0])

# 定义核函数
kernel = ConstantKernel(1.0) * RBF(length_scale=1.0)

# 创建高斯过程回归器
gpr = GaussianProcessRegressor(
    kernel=kernel,
    alpha=0.1,           # 噪声水平，用于数值稳定性
    n_restarts_optimizer=10,  # 优化器重启次数
    random_state=42
)

# 训练
gpr.fit(X, y)

# 预测（包含不确定性）
X_test = np.linspace(0, 12, 100).reshape(-1, 1)
y_pred, y_std = gpr.predict(X_test, return_std=True)

# 可视化
plt.figure(figsize=(10, 6))
plt.scatter(X, y, c='k', label='训练数据')
plt.plot(X_test, y_pred, 'b-', label='预测均值')
plt.fill_between(
    X_test.ravel(), 
    y_pred - 1.96 * y_std, 
    y_pred + 1.96 * y_std,
    alpha=0.2, color='b', label='95% 置信区间'
)
plt.xlabel('x')
plt.ylabel('y')
plt.title('高斯过程回归')
plt.legend()
plt.show()

print(f"优化后的核函数: {gpr.kernel_}")
print(f"对数边际似然: {gpr.log_marginal_likelihood_value_:.3f}")

核函数的作用​

核函数（协方差函数）决定了高斯过程的性质。它定义了两个输入点之间的相似度，进而决定了函数的平滑程度、周期性等特性。

核函数选择的影响：

from sklearn.gaussian_process.kernels import RBF, Matern, RationalQuadratic, ExpSineSquared

# 定义不同的核函数
kernels = {
    "RBF (平滑)": RBF(length_scale=1.0),
    "Matérn (ν=1.5)": Matern(length_scale=1.0, nu=1.5),
    "RationalQuadratic": RationalQuadratic(length_scale=1.0, alpha=1.0),
    "ExpSineSquared (周期)": ExpSineSquared(length_scale=1.0, periodicity=3.0)
}

plt.figure(figsize=(14, 10))

for i, (name, kernel) in enumerate(kernels.items(), 1):
    gpr = GaussianProcessRegressor(kernel=kernel, alpha=0.1, random_state=42)
    gpr.fit(X, y)
    y_pred, y_std = gpr.predict(X_test, return_std=True)
    
    plt.subplot(2, 2, i)
    plt.scatter(X, y, c='k', s=30)
    plt.plot(X_test, y_pred, 'b-')
    plt.fill_between(X_test.ravel(), y_pred - y_std, y_pred + y_std, alpha=0.2)
    plt.title(f'{name}')
    plt.xlabel('x')
    plt.ylabel('y')

plt.tight_layout()
plt.show()

处理噪声数据​

实际数据通常包含噪声。高斯过程可以通过两种方式处理噪声：

方式一：通过 alpha 参数

# alpha 参数控制噪声水平
# 较大的 alpha 值允许更多的噪声
gpr_noisy = GaussianProcessRegressor(
    kernel=RBF(),
    alpha=0.5,  # 假设噪声方差为 0.5
    random_state=42
)
gpr_noisy.fit(X, y)

方式二：通过 WhiteKernel 核

from sklearn.gaussian_process.kernels import WhiteKernel

# WhiteKernel 可以从数据中学习噪声水平
kernel_with_noise = RBF() + WhiteKernel(noise_level=0.1)

gpr_auto_noise = GaussianProcessRegressor(
    kernel=kernel_with_noise,
    n_restarts_optimizer=10,
    random_state=42
)
gpr_auto_noise.fit(X, y)

print(f"学习到的噪声水平: {gpr_auto_noise.kernel_.k2.noise_level:.4f}")

为什么需要处理噪声？

在数值计算中，核矩阵的条件数可能很大，导致计算不稳定。添加噪声相当于在核矩阵对角线上加一个小值（Tikhonov 正则化），可以提高数值稳定性。

超参数优化​

高斯过程通过最大化对数边际似然（Log Marginal Likelihood, LML）来优化核函数的超参数。

# 查看优化过程
gpr = GaussianProcessRegressor(
    kernel=RBF(length_scale=1.0) + WhiteKernel(noise_level=0.1),
    n_restarts_optimizer=10,
    random_state=42
)
gpr.fit(X, y)

print("优化前核函数:", RBF(length_scale=1.0) + WhiteKernel(noise_level=0.1))
print("优化后核函数:", gpr.kernel_)
print(f"对数边际似然: {gpr.log_marginal_likelihood_value_:.3f}")

# 手动计算不同超参数的对数边际似然
length_scales = np.logspace(-1, 1, 50)
lml_values = []

for ls in length_scales:
    kernel = RBF(length_scale=ls)
    gpr_temp = GaussianProcessRegressor(kernel=kernel, optimizer=None)
    gpr_temp.fit(X, y)
    lml_values.append(gpr_temp.log_marginal_likelihood_value_)

plt.figure(figsize=(8, 5))
plt.semilogx(length_scales, lml_values, 'b-')
plt.xlabel('length_scale')
plt.ylabel('对数边际似然')
plt.title('超参数优化：对数边际似然曲线')
plt.axvline(x=gpr.kernel_.k1.length_scale, color='r', linestyle='--', label='最优值')
plt.legend()
plt.grid(True)
plt.show()

从先验采样​

高斯过程的一个有趣特性是可以在没有训练数据的情况下从先验分布中采样函数：

# 定义核函数
kernel = RBF(length_scale=1.0)
gpr = GaussianProcessRegressor(kernel=kernel)

# 从先验采样（无需训练）
X_samples = np.linspace(0, 10, 100).reshape(-1, 1)
samples = gpr.sample_y(X_samples, n_samples=5, random_state=42)

plt.figure(figsize=(10, 6))
for i in range(5):
    plt.plot(X_samples, samples[:, i], alpha=0.7, label=f'样本 {i+1}')
plt.xlabel('x')
plt.ylabel('y')
plt.title('从高斯过程先验采样的函数')
plt.legend()
plt.show()

高斯过程分类​

GaussianProcessClassifier 实现了高斯过程分类，用于概率分类任务。与回归不同，分类问题的后验分布不是高斯分布，需要使用近似方法（如 Laplace 近似）。

原理简介​

对于二分类问题，高斯过程分类器：

  1. 在潜在函数 ff 上放置高斯过程先验
  2. 通过链接函数（如 logistic 函数）将潜在函数映射到概率
  3. 使用 Laplace 近似计算后验分布

潜在函数 ff 与类别概率的关系：

p(y=1∣x)=σ(f(x))=11+e−f(x)p(y=1|x) = \sigma(f(x)) = \frac{1}{1 + e^{-f(x)}}

基本用法​

from sklearn.gaussian_process import GaussianProcessClassifier
from sklearn.gaussian_process.kernels import RBF
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt

# 生成分类数据
X, y = make_classification(
    n_samples=100, n_features=2, n_redundant=0, 
    n_informative=2, random_state=42, n_clusters_per_class=1
)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 高斯过程分类
gpc = GaussianProcessClassifier(
    kernel=RBF(length_scale=1.0),
    random_state=42
)
gpc.fit(X_train, y_train)

# 预测概率
y_prob = gpc.predict_proba(X_test)

print(f"测试集准确率: {gpc.score(X_test, y_test):.2%}")

# 可视化决策边界
h = 0.02
x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))

Z = gpc.predict_proba(np.c_[xx.ravel(), yy.ravel()])[:, 1]
Z = Z.reshape(xx.shape)

plt.figure(figsize=(10, 6))
plt.contourf(xx, yy, Z, cmap='RdBu', alpha=0.8)
plt.scatter(X_train[:, 0], X_train[:, 1], c=y_train, cmap='RdBu', edgecolors='k')
plt.xlabel('特征 1')
plt.ylabel('特征 2')
plt.title('高斯过程分类 - 概率预测')
plt.colorbar(label='P(y=1)')
plt.show()

多分类策略​

对于多分类问题，GaussianProcessClassifier 支持两种策略：

一对多（One-vs-Rest）：为每个类别训练一个二分类器，将该类别与其他所有类别区分开。

一对一（One-vs-One）：为每对类别训练一个二分类器，通过投票决定最终类别。

from sklearn.datasets import load_iris

# 加载鸢尾花数据集（三分类）
iris = load_iris()
X, y = iris.data[:, :2], iris.target  # 只使用前两个特征便于可视化

# 一对多策略
gpc_ovr = GaussianProcessClassifier(
    kernel=RBF(),
    multi_class='one_vs_rest',
    random_state=42
)
gpc_ovr.fit(X, y)
print(f"一对多策略准确率: {gpc_ovr.score(X, y):.2%}")

# 一对一策略
gpc_ovo = GaussianProcessClassifier(
    kernel=RBF(),
    multi_class='one_vs_one',
    random_state=42
)
gpc_ovo.fit(X, y)
print(f"一对一策略准确率: {gpc_ovo.score(X, y):.2%}")

策略选择建议：

  • 一对多：类别数少时计算更快，支持概率输出
  • 一对一：类别数多时可能更快（每个子问题只涉及部分数据），但不支持概率输出

常用核函数​

核函数是高斯过程的核心组件，决定了模型的行为特性。sklearn 提供了丰富的核函数。

RBF 核（径向基函数核）​

RBF 核是最常用的核函数，也称为平方指数核。它产生的函数无限可微，非常平滑。

k(x,x′)=exp⁡(−d(x,x′)22l2)k(x, x') = \exp\left(-\frac{d(x, x')^2}{2l^2}\right)

其中 ll 是长度尺度参数。

from sklearn.gaussian_process.kernels import RBF

# 各向同性 RBF（所有维度使用相同长度尺度）
rbf_iso = RBF(length_scale=1.0)

# 各向异性 RBF（每个维度使用不同长度尺度）
rbf_aniso = RBF(length_scale=[1.0, 2.0, 0.5])  # 3维特征

# 使用示例
gpr = GaussianProcessRegressor(kernel=rbf_iso)

长度尺度的含义：

  • 较小的长度尺度：函数变化快，可以拟合高频模式
  • 较大的长度尺度：函数平滑，倾向于低频模式

Matérn 核​

Matérn 核是 RBF 核的推广，通过参数 ν\nu 控制函数的平滑程度。

k(x,x′)=21−νΓ(ν)(2νdl)νKν(2νdl)k(x, x') = \frac{2^{1-\nu}}{\Gamma(\nu)}\left(\frac{\sqrt{2\nu}d}{l}\right)^\nu K_\nu\left(\frac{\sqrt{2\nu}d}{l}\right)

其中 KνK_\nu 是修正贝塞尔函数。

from sklearn.gaussian_process.kernels import Matern

# 常用的 ν 值
matern_12 = Matern(nu=0.5)   # 绝对指数核，不可微
matern_32 = Matern(nu=1.5)   # 一阶可微
matern_52 = Matern(nu=2.5)   # 二阶可微
matern_inf = RBF()           # ν→∞ 时等价于 RBF

选择建议：

  • 如果目标函数足够平滑，使用 RBF 或 Matérn (ν=2.5)
  • 如果目标函数有突变或不连续，使用 Matérn (ν=0.5 或 1.5)

有理二次核​

有理二次核可以看作是不同长度尺度的 RBF 核的无限加权和，适合建模多尺度变化。

k(x,x′)=(1+d(x,x′)22αl2)−αk(x, x') = \left(1 + \frac{d(x, x')^2}{2\alpha l^2}\right)^{-\alpha}

from sklearn.gaussian_process.kernels import RationalQuadratic

# α 控制长度尺度的混合程度
# α → ∞ 时退化为 RBF
rq = RationalQuadratic(length_scale=1.0, alpha=1.0)

周期核（Exp-Sine-Squared）​

周期核用于建模周期性函数。

k(x,x′)=exp⁡(−2sin⁡2(πd/p)l2)k(x, x') = \exp\left(-\frac{2\sin^2(\pi d/p)}{l^2}\right)

其中 pp 是周期参数。

from sklearn.gaussian_process.kernels import ExpSineSquared

# periodicity 控制周期
periodic = ExpSineSquared(length_scale=1.0, periodicity=2*np.pi)

点积核​

点积核是非平稳核，适合建模线性函数或多项式函数。

k(x,x′)=σ02+x⋅x′k(x, x') = \sigma_0^2 + x \cdot x'

from sklearn.gaussian_process.kernels import DotProduct, ConstantKernel

# 线性核
linear = DotProduct(sigma_0=1.0)

# 多项式核（通过组合实现）
polynomial = ConstantKernel(1.0) * DotProduct(sigma_0=1.0) ** 2  # 二次多项式

核函数组合​

核函数可以通过加法和乘法组合，构建复杂的先验：

from sklearn.gaussian_process.kernels import RBF, ConstantKernel, WhiteKernel, ExpSineSquared

# 加法组合：建模多个独立信号源
kernel_sum = RBF(length_scale=1.0) + ExpSineSquared(length_scale=1.0, periodicity=5.0)

# 乘法组合：建模信号的调制
kernel_prod = ConstantKernel(1.0) * RBF(length_scale=1.0)

# 复杂组合
kernel_complex = (
    ConstantKernel(1.0) * RBF(length_scale=10.0) +  # 长期趋势
    RBF(length_scale=1.0) * ExpSineSquared(length_scale=1.0, periodicity=1.0) +  # 周期分量
    WhiteKernel(noise_level=0.1)  # 噪声
)

print("复杂核函数:", kernel_complex)

组合原则：

  • 加法：叠加不同模式的信号（如趋势 + 周期）
  • 乘法：调制信号（如周期信号的振幅随位置变化）

完整示例：贝叶斯优化​

高斯过程的一个重要应用是贝叶斯优化，用于高效地寻找黑盒函数的最优解。

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern
import matplotlib.pyplot as plt

# 定义目标函数（假设是黑盒函数）
def target_function(x):
    return np.sin(x) + 0.1 * np.cos(3 * x) - 0.05 * x

# 初始观测点
X_observed = np.array([[1.0], [3.0], [5.0], [7.0]])
y_observed = target_function(X_observed.ravel())

# 高斯过程模型
kernel = Matern(nu=2.5)
gpr = GaussianProcessRegressor(
    kernel=kernel,
    alpha=0.01,
    n_restarts_optimizer=5,
    random_state=42
)

# 贝叶斯优化迭代
n_iterations = 5
X_all = X_observed.copy()
y_all = y_observed.copy()

plt.figure(figsize=(15, 10))

for i in range(n_iterations):
    # 训练模型
    gpr.fit(X_all, y_all)
    
    # 在整个范围内预测
    X_grid = np.linspace(0, 10, 200).reshape(-1, 1)
    y_pred, y_std = gpr.predict(X_grid, return_std=True)
    
    # 计算采集函数（Expected Improvement）
    y_best = y_all.max()
    with np.errstate(divide='warn'):
        imp = y_pred - y_best
        Z = imp / y_std
        ei = imp * norm.cdf(Z) + y_std * norm.pdf(Z)
        ei[y_std < 1e-8] = 0.0
    
    # 选择下一个采样点
    next_x = X_grid[np.argmax(ei)]
    next_y = target_function(next_x.item())
    
    # 可视化
    plt.subplot(2, 3, i + 1)
    plt.plot(X_grid, target_function(X_grid.ravel()), 'r:', label='真实函数')
    plt.plot(X_grid, y_pred, 'b-', label='GP 预测')
    plt.fill_between(X_grid.ravel(), y_pred - y_std, y_pred + y_std, alpha=0.2)
    plt.scatter(X_all, y_all, c='k', s=50, zorder=5, label='观测点')
    plt.scatter([next_x], [next_y], c='g', s=100, marker='*', zorder=6, label='新采样点')
    plt.title(f'迭代 {i+1}')
    
    # 更新数据
    X_all = np.vstack([X_all, next_x.reshape(1, -1)])
    y_all = np.append(y_all, next_y)

plt.subplot(2, 3, 6)
plt.plot(X_grid, target_function(X_grid.ravel()), 'r:', label='真实函数')
plt.plot(X_grid, y_pred, 'b-', label='最终 GP')
plt.fill_between(X_grid.ravel(), y_pred - y_std, y_pred + y_std, alpha=0.2)
plt.scatter(X_all[:-1], y_all[:-1], c='k', s=50, zorder=5)
plt.scatter([X_all[-1]], [y_all[-1]], c='g', s=100, marker='*', zorder=6)
plt.title('最终结果')

plt.tight_layout()
plt.show()

print(f"找到的最大值: {y_all.max():.4f}")
print(f"对应的位置: {X_all[np.argmax(y_all)][0]:.4f}")

from scipy.stats import norm  # 确保导入

高斯过程 vs 其他方法​

特性高斯过程神经网络随机森林
不确定性估计原生支持需要额外方法需要额外方法
数据量要求小样本表现好需要大量数据中等数据量
计算复杂度O(n3)O(n^3)训练慢，预测快训练中等，预测快
可解释性核函数可解释黑盒特征重要性
外推能力差好差
适用场景小样本、需要不确定性大数据、复杂模式通用场景

最佳实践​

数据预处理​

from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

# 高斯过程对特征尺度敏感，建议标准化
pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('gpr', GaussianProcessRegressor(kernel=RBF()))
])

核函数选择​

数据特点推荐核函数
平滑函数RBF
有突变或不连续Matérn (低 ν)
周期性数据ExpSineSquared
多尺度变化RationalQuadratic
线性趋势DotProduct 或 RBF + DotProduct
混合模式核函数组合

计算优化​

对于大数据集，高斯过程的计算成本很高。可以考虑：

# 1. 使用稀疏近似（sklearn 不直接支持，可以使用 GPy 或 GPyTorch）

# 2. 限制训练数据
from sklearn.model_selection import train_test_split

X_subset, _, y_subset, _ = train_test_split(
    X_train, y_train, train_size=1000, random_state=42
)

# 3. 调整核函数复杂度
simple_kernel = RBF(length_scale=1.0)  # 简单核计算更快

小结​

高斯过程是一种强大而优雅的机器学习方法：

  1. 概率输出：不仅给出预测值，还提供不确定性估计
  2. 核函数设计：通过核函数编码先验知识，灵活适应不同问题
  3. 小样本优势：在数据稀缺时表现优异
  4. 计算成本：O(n3)O(n^3) 的复杂度限制了其在大数据上的应用
  5. 应用场景：贝叶斯优化、主动学习、传感器融合、时间序列预测等

高斯过程最适合中小规模数据集且需要不确定性估计的场景。当数据量很大时，建议考虑其他方法或使用稀疏高斯过程近似。