第7章：稀疏核机（Sparse Kernel Machines）
一、导读
7.1 本章学习目标
第7章聚焦于核方法的一个关键问题：稀疏性。第6章中介绍的核方法（如高斯过程回归、核岭回归）在测试阶段需要计算与所有训练数据点的核函数值，这在训练集很大时计算代价高昂。稀疏核机通过选择少量代表性数据点（相关向量）来构建模型，在保持预测精度的同时大幅降低计算成本。
本章的核心学习目标包括：
- 理解稀疏性的动机：为什么我们需要稀疏的核模型？全核模型的计算瓶颈在哪里？
- 掌握相关向量机（RVM）的原理：RVM如何通过自动相关性确定（ARD）实现稀疏性？
- 理解证据近似与模型选择：如何自动确定超参数？边缘似然在模型选择中的作用。
- 掌握稀疏性的度量：如何衡量模型的稀疏程度？稀疏性与泛化能力有什么关系？
- 理解RVM与SVM的区别：两种稀疏核机在设计哲学和实际表现上的异同。
- 了解贝叶斯模型比较：如何使用边缘似然比较不同核函数和超参数？
7.2 为什么本章如此重要
稀疏性是机器学习中最重要的概念之一。一个稀疏模型只使用少量参数或数据点，具有计算效率高、可解释性强、不容易过拟合等优点。相关向量机展示了如何通过贝叶斯方法自然地获得稀疏性——不需要显式的L1正则化，而是通过自动相关性确定（ARD）先验让不相关的基函数的权重自动趋近于零。
RVM的理论贡献在于：它证明了贝叶斯框架可以自然地产生稀疏解，而且通常比SVM更稀疏。这为后续的稀疏贝叶斯学习方法奠定了基础。
7.3 核心问题引导
在阅读本章时，请带着以下问题思考：
为什么核岭回归的解不是稀疏的？什么条件下解会变得稀疏？
ARD先验如何导致稀疏性？它与L1正则化有什么本质区别？
相关向量机为什么通常比SVM产生更稀疏的模型？
证据近似中，超参数的优化为什么不会过拟合？
如何自动选择相关向量？选择标准是什么？
RVM的预测概率与SVM的决策函数有什么本质区别？
二、核心概念详解
7.1 稀疏性的动机
7.1.1 全核模型的计算瓶颈
在第6章中，核岭回归的预测为：
$$y(\mathbf{x}) = \sum_{n=1}^{N} a_n k(\mathbf{x}, \mathbf{x}_n)$$
测试阶段需要计算与所有 $N$ 个训练点的核函数值。当 $N$ 很大时（如 $N=100000$），每次预测需要计算 $100000$ 次核函数，这在实时应用中是不可接受的。
此外，训练阶段需要求解 $(N \times N)$ 的线性系统，复杂度为 $O(N^3)$。
7.1.2 稀疏性的目标
稀疏核机的目标是找到一个模型：
$$y(\mathbf{x}) = \sum_{n \in \mathcal{S}} a_n k(\mathbf{x}, \mathbf{x}_n)$$
其中 $\mathcal{S} \subset \{1, \ldots, N\}$ 是一个小集合（$|\mathcal{S}| \ll N$），对应的数据点称为相关向量（Relevance Vectors）。
稀疏性的好处：
- 测试速度快：只需要计算与 $|\mathcal{S}|$ 个点的核函数
- 存储效率高：只需要存储 $|\mathcal{S}|$ 个数据点
- 可解释性：相关向量可以被视为数据的代表性样本
7.2 相关向量机（RVM）
7.2.1 模型设定
相关向量机（Relevance Vector Machine, RVM）是Tipping（2001）提出的贝叶斯稀疏核方法。
对于回归问题，RVM的模型为：
$$y(\mathbf{x}) = \sum_{i=1}^{M} w_i \phi_i(\mathbf{x}) + \epsilon$$
其中 $\phi_i(\mathbf{x}) = k(\mathbf{x}, \mathbf{x}_i)$ 是以训练点为中心的核函数，$\epsilon \sim \mathcal{N}(0, \sigma^2)$ 是噪声。
7.2.2 ARD先验
RVM的关键创新是自动相关性确定（Automatic Relevance Determination, ARD）先验。每个权重 $w_i$ 有独立的先验：
$$p(w_i|\alpha_i) = \mathcal{N}(w_i|0, \alpha_i^{-1})$$
其中 $\alpha_i$ 是每个权重独立的超参数（精度）。
在训练过程中，大部分 $\alpha_i$ 会趋向无穷大，对应的权重 $w_i$ 精确为零。只有少量 $\alpha_i$ 保持有限值，对应的权重非零，这些训练点成为相关向量。
7.2.3 后验分布
给定数据，权重的后验分布为：
$$p(\mathbf{w}|\mathbf{t}, \boldsymbol{\alpha}, \sigma^2) = \mathcal{N}(\mathbf{w}|\mathbf{m}, \boldsymbol{\Sigma})$$
其中：
- $\boldsymbol{\Sigma} = (\sigma^{-2} \mathbf{\Phi}^T \mathbf{\Phi} + \mathbf{A})^{-1}$
- $\mathbf{m} = \sigma^{-2} \boldsymbol{\Sigma} \mathbf{\Phi}^T \mathbf{t}$
- $\mathbf{A} = \text{diag}(\alpha_1, \ldots, \alpha_M)$
7.2.4 超参数优化
通过最大化边缘似然 $p(\mathbf{t}|\boldsymbol{\alpha}, \sigma^2)$ 来确定超参数：
$$\ln p(\mathbf{t}|\boldsymbol{\alpha}, \sigma^2) = -\frac{1}{2} \left[N \ln(2\pi) + \ln|\mathbf{C}| + \mathbf{t}^T \mathbf{C}^{-1} \mathbf{t}\right]$$
其中 $\mathbf{C} = \sigma^2 \mathbf{I} + \mathbf{\Phi} \mathbf{A}^{-1} \mathbf{\Phi}^T$。
对 $\alpha_i$ 的更新公式：
$$\alpha_i^{new} = \frac{\gamma_i}{m_i^2}$$
其中 $\gamma_i = 1 - \alpha_i \Sigma_{ii}$，衡量了第 $i$ 个基函数对模型的"相关程度"。
import numpy as np
import matplotlib.pyplot as plt
class RelevanceVectorMachine:
"""相关向量机（回归）"""
def __init__(self, kernel_func, max_iter=100, tol=1e-5):
self.kernel_func = kernel_func
self.max_iter = max_iter
self.tol = tol
def fit(self, X, t):
N = len(X)
# 初始化
# 基函数：以每个训练点为中心的核函数
self.X_train = X.copy()
Phi = self.kernel_func(X, X) # N x N
# 初始化超参数
alpha = np.ones(N) * 0.01
sigma2 = np.var(t) * 0.1
for iteration in range(self.max_iter):
# 后验参数
A = np.diag(alpha)
Sigma_inv = Phi.T @ Phi / sigma2 + A
Sigma = np.linalg.inv(Sigma_inv)
m = Phi.T @ t / sigma2
m = Sigma @ m
# 更新alpha
gamma = 1 - alpha * np.diag(Sigma)
gamma = np.maximum(gamma, 1e-10)
alpha_new = gamma / (m**2 + 1e-10)
# 更新sigma2
sigma2_new = np.sum((t - Phi @ m)**2) / (N - np.sum(gamma) + 1e-10)
sigma2_new = max(sigma2_new, 1e-10)
# 检查收敛
alpha_change = np.max(np.abs(alpha_new - alpha) / (alpha + 1e-10))
sigma2_change = abs(sigma2_new - sigma2) / (sigma2 + 1e-10)
alpha = alpha_new
sigma2 = sigma2_new
if max(alpha_change, sigma2_change) < self.tol:
print(f"收敛于第 {iteration+1} 次迭代")
break
# 保存相关向量
self.relevance_indices = np.where(alpha < 1e6)[0]
self.alpha = alpha[self.relevance_indices]
self.m = m[self.relevance_indices]
self.sigma2 = sigma2
self.X_relevance = X[self.relevance_indices]
print(f"相关向量数量: {len(self.relevance_indices)} / {N}")
print(f"稀疏度: {1 - len(self.relevance_indices)/N:.2%}")
return self
def predict(self, X):
Phi = self.kernel_func(X, self.X_relevance)
return Phi @ self.m
def predict_with_uncertainty(self, X):
Phi = self.kernel_func(X, self.X_relevance)
mean = Phi @ self.m
# 预测方差
A_inv = np.diag(1.0 / self.alpha)
# 近似后验协方差
Phi_train = self.kernel_func(self.X_relevance, self.X_relevance)
Sigma = np.linalg.inv(Phi_train.T @ Phi_train / self.sigma2 + np.diag(self.alpha))
var = np.array([self.sigma2 + Phi[i] @ Sigma @ Phi[i] for i in range(len(X))])
return mean, var
def rbf_kernel_rvm(x1, x2, gamma=1.0):
"""RBF核函数"""
sq_dists = np.sum(x1**2, axis=1, keepdims=True) + np.sum(x2**2, axis=1) - 2 * x1 @ x2.T
return np.exp(-gamma * sq_dists)
# 示例
np.random.seed(42)
N = 50
X_rvm = np.sort(np.random.uniform(-3, 3, N)).reshape(-1, 1)
t_rvm = np.sin(X_rvm).ravel() + 0.1 * np.random.randn(N)
rvm = RelevanceVectorMachine(kernel_func=lambda x1, x2: rbf_kernel_rvm(x1, x2, gamma=2.0))
rvm.fit(X_rvm, t_rvm)
# 预测
X_test_rvm = np.linspace(-4, 4, 200).reshape(-1, 1)
mean_rvm, var_rvm = rvm.predict_with_uncertainty(X_test_rvm)
std_rvm = np.sqrt(var_rvm)
# 可视化
fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(X_rvm, t_rvm, c='red', s=30, zorder=5, label='Training data')
ax.scatter(rvm.X_relevance, rvm.predict(rvm.X_relevance), c='blue', s=80,
marker='x', zorder=6, label=f'Relevance vectors ({len(rvm.relevance_indices)})')
ax.plot(X_test_rvm, mean_rvm, 'b-', linewidth=2, label='RVM prediction')
ax.fill_between(X_test_rvm.ravel(), mean_rvm - 2*std_rvm, mean_rvm + 2*std_rvm,
alpha=0.2, color='blue', label='95% CI')
ax.plot(X_test_rvm, np.sin(X_test_rvm), 'g--', linewidth=2, label='True function')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Relevance Vector Machine')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('rvm.png', dpi=150)
plt.close()7.3 RVM vs SVM
7.3.1 设计哲学对比
| 特性 | SVM | RVM |
|---|---|---|
| 框架 | 频率学派（结构风险最小化） | 贝叶斯 |
| 稀疏性来源 | 间隔最大化（几何） | ARD先验（概率） |
| 超参数选择 | 交叉验证 | 边缘似然 |
| 输出 | 决策函数（无概率） | 预测概率 |
| 核函数约束 | Mercer条件 | 不需要 |
| 稀疏程度 | 通常较多支持向量 | 通常更稀疏 |
7.3.2 性能比较
def compare_rvm_svm():
"""比较RVM和SVM的性能"""
from sklearn.svm import SVR
from sklearn.model_selection import cross_val_score
np.random.seed(42)
# 生成数据
N = 100
X = np.sort(np.random.uniform(-3, 3, N)).reshape(-1, 1)
t = np.sin(X).ravel() + 0.2 * np.random.randn(N)
# RVM
rvm = RelevanceVectorMachine(kernel_func=lambda x1, x2: rbf_kernel_rvm(x1, x2, gamma=1.0))
rvm.fit(X, t)
n_rv = len(rvm.relevance_indices)
# SVM
svm = SVR(kernel='rbf', gamma=1.0, C=10)
svm.fit(X, t)
n_sv = np.sum(svm.support_)
# 测试精度
X_test = np.linspace(-3, 3, 200).reshape(-1, 1)
t_true = np.sin(X_test).ravel()
rvm_pred = rvm.predict(X_test)
svm_pred = svm.predict(X_test)
rvm_mse = np.mean((rvm_pred - t_true) ** 2)
svm_mse = np.mean((svm_pred - t_true) ** 2)
print(f"RVM: 相关向量={n_rv}, MSE={rvm_mse:.4f}")
print(f"SVM: 支持向量={n_sv}, MSE={svm_mse:.4f}")
print(f"\nRVM稀疏度: {1 - n_rv/N:.2%}")
print(f"SVM稀疏度: {1 - n_sv/N:.2%}")
compare_rvm_svm()7.4 证据近似与模型选择
7.4.1 类型II最大似然
RVM使用类型II最大似然（Type-II Maximum Likelihood）或证据近似来确定超参数。核心思想是：
给定超参数，计算参数的后验分布
通过积分掉参数，得到边缘似然 $p(\mathbf{t}|\boldsymbol{\alpha}, \sigma^2)$
最大化边缘似然来确定超参数
边缘似然关于超参数的优化不会过拟合，因为：
- 边缘似然自动平衡了模型拟合度和复杂度
- 过于复杂的模型（超参数过多）会被边缘似然惩罚
7.4.2 贝叶斯模型比较
边缘似然可以用于比较不同的模型（不同的核函数、不同的超参数设置）：
$$p(\mathcal{M}_i|\mathbf{t}) \propto p(\mathbf{t}|\mathcal{M}_i) p(\mathcal{M}_i)$$
选择边缘似然最大的模型。
def bayesian_model_comparison_rvm():
"""使用边缘似然比较不同gamma值的RVM"""
np.random.seed(42)
N = 50
X = np.sort(np.random.uniform(-3, 3, N)).reshape(-1, 1)
t = np.sin(X).ravel() + 0.2 * np.random.randn(N)
gammas = [0.1, 0.5, 1.0, 2.0, 5.0, 10.0]
log_evidences = []
for gamma in gammas:
Phi = rbf_kernel_rvm(X, X, gamma)
# 简单的证据计算
alpha = np.ones(N) * 0.1
sigma2 = np.var(t) * 0.1
for _ in range(20):
A = np.diag(alpha)
Sigma = np.linalg.inv(Phi.T @ Phi / sigma2 + A)
m = Sigma @ Phi.T @ t / sigma2
gamma_i = 1 - alpha * np.diag(Sigma)
gamma_i = np.maximum(gamma_i, 1e-10)
alpha = gamma_i / (m**2 + 1e-10)
sigma2 = np.sum((t - Phi @ m)**2) / (N - np.sum(gamma_i) + 1e-10)
sigma2 = max(sigma2, 1e-10)
# 对数边缘似然
C = sigma2 * np.eye(N) + Phi @ np.diag(1.0/alpha) @ Phi.T
log_evidence = -0.5 * (N * np.log(2*np.pi) + np.log(np.linalg.det(C) + 1e-300)
+ t @ np.linalg.inv(C + 1e-6*np.eye(N)) @ t)
log_evidences.append(log_evidence)
n_rv = np.sum(alpha < 1e6)
print(f"γ={gamma:5.1f}: ln(evidence)={log_evidence:.2f}, 相关向量={n_rv}")
best_gamma = gammas[np.argmax(log_evidences)]
print(f"\n最优gamma: {best_gamma}")
bayesian_model_comparison_rvm()7.5 稀疏性的分析
7.5.1 为什么ARD产生稀疏性
ARD先验的关键在于每个权重有独立的超参数 $\alpha_i$。在边缘似然优化过程中：
如果第 $i$ 个基函数对解释数据的贡献很小，增大 $\alpha_i$（减小权重先验方差）会减小边缘似然中的复杂度惩罚
当 $\alpha_i \to \infty$ 时，$w_i$ 被精确地压缩到零
只有对数据有显著贡献的基函数，其 $\alpha_i$ 保持有限值
这种机制类似于"赢家通吃"：只有少数最相关的基函数被保留。
7.5.2 稀疏性与泛化
稀疏模型通常具有更好的泛化能力，原因有三：
参数少：减少了过拟合的风险
自动模型选择：ARD先验自动选择了合适的模型复杂度
代表性样本：相关向量捕获了数据的最重要特征
def analyze_sparsity():
"""分析稀疏性与泛化的关系"""
np.random.seed(42)
# 不同训练集大小下的稀疏度
N_values = [20, 50, 100, 200]
print(f"{'N':>5} {'相关向量':>8} {'稀疏度':>8}")
print("-" * 30)
for N in N_values:
X = np.sort(np.random.uniform(-3, 3, N)).reshape(-1, 1)
t = np.sin(X).ravel() + 0.2 * np.random.randn(N)
rvm = RelevanceVectorMachine(
kernel_func=lambda x1, x2: rbf_kernel_rvm(x1, x2, gamma=1.0),
max_iter=50
)
rvm.fit(X, t)
n_rv = len(rvm.relevance_indices)
sparsity = 1 - n_rv / N
print(f"{N:5d} {n_rv:8d} {sparsity:8.2%}")
analyze_sparsity()三、实战案例
7.1 回归任务中的RVM
def rvm_regression_case_study():
"""RVM在实际回归任务中的应用"""
np.random.seed(42)
# 生成带有趋势和季节性的数据
N = 100
X = np.linspace(0, 10, N).reshape(-1, 1)
t = 0.5 * X.ravel() + 2 * np.sin(2 * np.pi * X.ravel() / 3) + np.random.randn(N) * 0.5
# 训练RVM
rvm = RelevanceVectorMachine(
kernel_func=lambda x1, x2: rbf_kernel_rvm(x1, x2, gamma=0.5),
max_iter=100
)
rvm.fit(X, t)
# 预测
X_test = np.linspace(-1, 11, 300).reshape(-1, 1)
mean, var = rvm.predict_with_uncertainty(X_test)
std = np.sqrt(var)
# 可视化
fig, ax = plt.subplots(figsize=(12, 6))
ax.scatter(X, t, c='red', s=20, alpha=0.5, label='Training data')
ax.scatter(rvm.X_relevance, t[rvm.relevance_indices], c='blue', s=60,
marker='*', zorder=5, label=f'Relevance vectors ({len(rvm.relevance_indices)})')
ax.plot(X_test, mean, 'b-', linewidth=2, label='RVM prediction')
ax.fill_between(X_test.ravel(), mean - 2*std, mean + 2*std,
alpha=0.15, color='blue', label='95% CI')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('RVM Regression with Trend and Seasonality')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('rvm_regression.png', dpi=150)
plt.close()
rvm_regression_case_study()7.2 分类任务中的RVM
def rvm_classification_case_study():
"""RVM在分类任务中的应用"""
np.random.seed(42)
# 生成二分类数据
N = 100
X_cls = np.vstack([
np.random.randn(N//2, 2) * 0.8 + [1.5, 1.5],
np.random.randn(N//2, 2) * 0.8 + [-1.5, -1.5]
])
t_cls = np.array([1]*(N//2) + [0]*(N//2))
# 简化的RVM分类（使用逻辑sigmoid似然）
# 这里使用sklearn的近似实现
from sklearn.linear_model import LogisticRegression
# 使用核特征近似
gamma = 0.5
K_train = rbf_kernel_rvm(X_cls, X_cls, gamma)
# 使用核PCA进行降维
from sklearn.decomposition import KernelPCA
kpca = KernelPCA(n_components=20, kernel='rbf', gamma=gamma)
X_kpca = kpca.fit_transform(X_cls)
# 带L1正则化的逻辑回归（产生稀疏解）
lr = LogisticRegression(penalty='l1', C=1.0, solver='saga', max_iter=1000)
lr.fit(X_kpca, t_cls)
# 找到非零系数对应的原始数据点
non_zero = np.any(lr.coef_ != 0, axis=0)
print(f"非零特征数: {np.sum(non_zero)} / {X_kpca.shape[1]}")
accuracy = lr.score(X_kpca, t_cls)
print(f"训练准确率: {accuracy:.4f}")
rvm_classification_case_study()7.3 模型复杂度分析
def model_complexity_analysis():
"""分析模型复杂度与泛化的关系"""
np.random.seed(42)
# 固定测试集
N_test = 200
X_test = np.linspace(-3, 3, N_test).reshape(-1, 1)
t_test = np.sin(X_test).ravel()
# 不同训练集大小
N_values = [10, 20, 50, 100, 200]
train_errors = []
test_errors = []
n_rvs = []
for N in N_values:
X = np.sort(np.random.uniform(-3, 3, N)).reshape(-1, 1)
t = np.sin(X).ravel() + 0.2 * np.random.randn(N)
rvm = RelevanceVectorMachine(
kernel_func=lambda x1, x2: rbf_kernel_rvm(x1, x2, gamma=1.0),
max_iter=50
)
rvm.fit(X, t)
train_pred = rvm.predict(X)
test_pred = rvm.predict(X_test)
train_err = np.sqrt(np.mean((train_pred - t)**2))
test_err = np.sqrt(np.mean((test_pred - t_test)**2))
train_errors.append(train_err)
test_errors.append(test_err)
n_rvs.append(len(rvm.relevance_indices))
print(f"{'N':>5} {'训练RMSE':>10} {'测试RMSE':>10} {'相关向量':>8}")
print("-" * 40)
for N, tr, te, nr in zip(N_values, train_errors, test_errors, n_rvs):
print(f"{N:5d} {tr:10.4f} {te:10.4f} {nr:8d}")
model_complexity_analysis()四、本章小结
第7章系统介绍了稀疏核机的理论和应用：
稀疏性的动机：全核模型在测试阶段需要计算与所有训练点的核函数，计算代价高昂。稀疏核机通过选择少量相关向量来降低计算成本。
相关向量机（RVM）：通过ARD先验为每个权重分配独立的超参数。在边缘似然优化过程中，大部分超参数趋向无穷大，对应的权重精确为零，只有少量相关向量被保留。
RVM vs SVM：RVM基于贝叶斯框架，通常比SVM更稀疏，且提供预测概率。SVM基于结构风险最小化，有更强的理论保证。
证据近似：通过最大化边缘似然确定超参数，自动平衡拟合度和复杂度，不会过拟合。
稀疏性与泛化：稀疏模型参数少、自动选择复杂度，通常具有更好的泛化能力。
稀疏核机展示了贝叶斯方法在模型选择中的强大能力。ARD先验提供了一种自然的方式来实现稀疏性，不需要显式的L1正则化。这种思想在后续的稀疏贝叶斯学习、压缩感知等领域得到了广泛应用。
五、思考与练习
思考题
RVM的ARD先验和Lasso的L1正则化都能产生稀疏性，它们的机制有什么本质区别？
为什么RVM通常比SVM产生更稀疏的模型？
边缘似然为什么不会过拟合？请从偏差-方差分解的角度分析。
当核函数不满足Mercer条件时，SVM不能使用但RVM仍然可以使用，为什么？
RVM的预测不确定性在哪些实际应用中很重要？
编程练习
不同核函数的RVM：实现使用多项式核和Laplacian核的RVM，比较它们的稀疏度和预测精度。
增量RVM：实现增量版本的RVM，逐步添加和删除基函数，观察相关向量的变化过程。
RVM用于时间序列预测：使用RVM对一个时间序列进行预测，分析相关向量在时间轴上的分布。