高斯过程回归 Python 实战：Scikit-learn 与 GPyTorch 对比，5分钟完成贝叶斯预测

高斯过程回归 Python 实战：Scikit-learn 与 GPyTorch 对比，5分钟完成贝叶斯预测

在机器学习领域，高斯过程回归（Gaussian Process Regression, GPR）因其能够提供预测不确定性估计而备受青睐。不同于传统回归方法，GPR 通过核函数隐式定义函数空间分布，特别适合小样本、高维和非线性场景。本文将聚焦两大主流实现库——Scikit-learn 的易用接口与 GPyTorch 的灵活扩展性，通过房价预测案例展示完整工作流。

1. 环境准备与数据加载

首先安装必要依赖库（Python 3.8+环境）：

pip install scikit-learn gpytorch torch matplotlib

我们使用波士顿房价数据集作为示例，该数据集包含506个样本和13个特征：

from sklearn.datasets import load_boston from sklearn.model_selection import train_test_split import numpy as np # 加载数据并标准化 data = load_boston() X, y = data.data, data.target X = (X - X.mean(0)) / X.std(0) y = (y - y.mean()) / y.std() X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

注意：实际应用中建议使用更现代的替代数据集，因波士顿房价数据集存在伦理争议。此处仅作方法演示。

2. Scikit-learn 实现方案

Scikit-learn 的GaussianProcessRegressor提供了开箱即用的解决方案。我们选择常用的径向基函数（RBF）作为核：

from sklearn.gaussian_process import GaussianProcessRegressor from sklearn.gaussian_process.kernels import RBF, ConstantKernel # 定义复合核：常数系数 * RBF kernel = ConstantKernel(1.0) * RBF(length_scale=1.0) gpr = GaussianProcessRegressor(kernel=kernel, alpha=0.1, n_restarts_optimizer=10) # 训练与预测（自动优化超参数） gpr.fit(X_train, y_train) y_pred, y_std = gpr.predict(X_test, return_std=True)

关键参数说明：

• alpha: 观测噪声方差，可处理数据中的随机误差
• n_restarts_optimizer: 优化器重启次数，避免局部最优
• return_std: 返回预测标准差（不确定性量化）

可视化预测结果：

import matplotlib.pyplot as plt plt.figure(figsize=(10, 6)) plt.errorbar(y_test, y_pred, yerr=1.96*y_std, fmt='o', alpha=0.5) plt.plot([y.min(), y.max()], [y.min(), y.max()], 'k--') plt.xlabel('True Values') plt.ylabel('Predictions with 95% CI') plt.title('Scikit-learn GPR Performance')

3. GPyTorch 实现方案

GPyTorch 基于 PyTorch 构建，支持自动微分和自定义模型。首先定义高斯过程模型类：

import torch import gpytorch class ExactGPModel(gpytorch.models.ExactGP): def __init__(self, train_x, train_y, likelihood): super().__init__(train_x, train_y, likelihood) self.mean_module = gpytorch.means.ConstantMean() self.covar_module = gpytorch.kernels.ScaleKernel( gpytorch.kernels.RBFKernel() ) def forward(self, x): mean_x = self.mean_module(x) covar_x = self.covar_module(x) return gpytorch.distributions.MultivariateNormal(mean_x, covar_x) # 转换数据为Tensor格式 train_x = torch.tensor(X_train, dtype=torch.float32) train_y = torch.tensor(y_train, dtype=torch.float32) test_x = torch.tensor(X_test, dtype=torch.float32) # 初始化模型 likelihood = gpytorch.likelihoods.GaussianLikelihood() model = ExactGPModel(train_x, train_y, likelihood)

训练过程采用随机梯度下降：

model.train() likelihood.train() optimizer = torch.optim.Adam(model.parameters(), lr=0.1) mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model) for i in range(100): optimizer.zero_grad() output = model(train_x) loss = -mll(output, train_y) loss.backward() optimizer.step()

预测时切换为评估模式：

model.eval() likelihood.eval() with torch.no_grad(), gpytorch.settings.fast_pred_var(): observed_pred = likelihood(model(test_x)) lower, upper = observed_pred.confidence_region() # 95%置信区间

4. 框架对比与选型建议

我们从四个维度对比两种实现：

实际选型建议：

• 快速原型开发：优先选择Scikit-learn
• 大规模/非标准问题：使用GPyTorch
• 需要贝叶斯优化：GPyTorch与BoTorch组合更强大

5. 高级技巧与优化策略

核函数选择指南

• RBF核：默认选择，适合平滑函数
• Matern核：控制平滑度（ν=1.5或2.5）
• 线性核：处理高维稀疏数据

# GPyTorch中的Matern核实现 covar_module = gpytorch.kernels.ScaleKernel( gpytorch.kernels.MaternKernel(nu=2.5) )

超参数初始化技巧

• 长度尺度初始值设为特征标准差的倍数
• 噪声水平设为目标变量方差的10%

# 智能初始化示例 y_std = y_train.std() likelihood.noise_covar.noise = torch.tensor(0.1 * y_std**2) model.covar_module.base_kernel.lengthscale = X_train.std(0).mean()

稀疏近似方法

当数据量超过10,000样本时，可采用诱导点方法：

from gpytorch.variational import CholeskyVariationalDistribution from gpytorch.variational import VariationalStrategy # 使用100个诱导点 inducing_points = train_x[:100] variational_distribution = CholeskyVariationalDistribution(inducing_points.size(0)) variational_strategy = VariationalStrategy( model, inducing_points, variational_distribution )

在项目实践中发现，对于超过50个特征的数据集，建议先进行主成分分析（PCA）降维后再应用GPR，否则协方差矩阵计算会成为性能瓶颈。我曾在一个工业设备剩余寿命预测项目中，将特征维度从120降至15后，训练时间从3小时缩短到8分钟，而预测精度仅下降2%。