核感知机与谱聚类

概述

核感知机（Kernel Perceptron）和谱聚类（Spectral Clustering）代表了核方法的两个重要应用方向：分类和聚类。核感知机将线性感知机扩展到非线性情形，而谱聚类利用图谱理论进行流形聚类。¹

共同主题：通过核函数将数据映射到高维特征空间，在那里线性方法变得强大。

核感知机

标准感知机回顾

感知机算法（Rosenblatt, 1957）：

def perceptron_train(X, y, max_iter=1000):
    """标准感知机训练"""
    n, d = X.shape
    w = np.zeros(d)  # 权重
    b = 0            # 偏置
    
    for _ in range(max_iter):
        mistakes = 0
        for i in range(n):
            if y[i] * (X[i] @ w + b) <= 0:
                w += y[i] * X[i]
                b += y[i]
                mistakes += 1
        if mistakes == 0:
            break
    
    return w, b
 
def perceptron_predict(X, w, b):
    """预测"""
    return np.sign(X @ w + b)

线性不可分问题

标准感知机只能在线性可分数据上收敛。当数据线性不可分时，感知机不会收敛。

线性可分：                  线性不可分：
    ○ ╲                          ○ ●
      ╲  ●                   ●       ○
      ● ╱                   ○   ●
    ○                       ● ○
    
    感知机会收敛              感知机不收敛

核感知机的动机

通过核函数将数据映射到高维特征空间 $ϕ (x)$ ，使得原本线性不可分的数据变得线性可分：

原始空间 (低维)          特征空间 (高维)
      │                        │
      │   φ: x → φ(x)         │
      ▼                        ▼
   非线性边界              线性边界

核感知机算法

核感知机将权重向量表示为特征空间中数据点的线性组合：

f (x) = w^{T} ϕ (x) + b = i = 1 \sum n α_{i} y_{i} K (x_{i}, x) + b

其中 $α_{i}$ 是组合系数。

核感知机更新规则：

当样本 $(x_{i}, y_{i})$ 被错误分类时：

α_{i} \leftarrow α_{i} + 1

这等价于在特征空间中的更新：

w \leftarrow w + y_{i} ϕ (x_{i})

核感知机实现

import numpy as np
 
class KernelPerceptron:
    """核感知机"""
    
    def __init__(self, kernel='rbf', gamma=0.1, **kernel_params):
        self.kernel = kernel
        self.gamma = gamma
        self.kernel_params = kernel_params
        self.alpha = None
        self.support_vectors = None
        self.y_sv = None
        self.b = 0.0
        
    def kernel_function(self, X1, X2):
        """核函数"""
        if self.kernel == 'linear':
            return X1 @ X2.T
        elif self.kernel == 'rbf':
            # 高斯核
            pairwise_sq = np.sum(X1**2, axis=1, keepdims=True) + \
                          np.sum(X2**2, axis=1) - \
                          2 * (X1 @ X2.T)
            return np.exp(-self.gamma * pairwise_sq)
        elif self.kernel == 'poly':
            p = self.kernel_params.get('p', 2)
            return (X1 @ X2.T + 1) ** p
        else:
            raise ValueError(f"Unknown kernel: {self.kernel}")
    
    def fit(self, X, y, max_iter=1000):
        """
        训练核感知机
        """
        n = X.shape[0]
        self.alpha = np.zeros(n)
        self.support_vectors = X.copy()
        self.y_sv = y.copy()
        
        # 计算核矩阵
        K = self.kernel_function(X, X)
        
        for iteration in range(max_iter):
            mistakes = 0
            
            for i in range(n):
                # 预测
                pred = np.sign(np.sum(self.alpha * y * K[:, i]) + self.b)
                
                if pred * y[i] <= 0:
                    # 更新：增加该样本的权重
                    self.alpha[i] += 1
                    mistakes += 1
            
            if mistakes == 0:
                print(f"在第 {iteration+1} 次迭代收敛")
                break
        
        # 只保留非零alpha的样本作为支持向量
        sv_idx = self.alpha > 1e-7
        self.support_vectors = X[sv_idx]
        self.y_sv = y[sv_idx]
        self.alpha = self.alpha[sv_idx]
        
        return self
    
    def predict(self, X):
        """预测"""
        K = self.kernel_function(X, self.support_vectors)
        return np.sign(np.sum(self.alpha * self.y_sv * K, axis=1))
 
# 测试
np.random.seed(42)
# 生成非线性可分的圆形数据
t1 = np.random.uniform(0, 2*np.pi, 100)
inner_circle = np.column_stack([np.cos(t1), np.sin(t1)]) + 0.1*np.random.randn(100, 2)
label1 = np.ones(100)
 
t2 = np.random.uniform(0, 2*np.pi, 100)
outer_circle = np.column_stack([3*np.cos(t2), 3*np.sin(t2)]) + 0.2*np.random.randn(100, 2)
label2 = -np.ones(100)
 
X = np.vstack([inner_circle, outer_circle])
y = np.concatenate([label1, label2])
 
# 训练核感知机
kp = KernelPerceptron(kernel='rbf', gamma=0.5)
kp.fit(X, y, max_iter=100)
 
# 评估
pred = kp.predict(X)
accuracy = np.mean(pred == y)
print(f"训练准确率: {accuracy:.2%}")
print(f"支持向量数量: {len(kp.alpha)}")

核感知机的理论性质

性质	描述
收敛性	在特征空间线性可分时收敛
支持向量	权重非零的样本构成支持向量
与SVM的关系	感知机找到分离超平面，不保证最大间隔
过拟合风险	线性不可分时可能过拟合

谱聚类（Spectral Clustering）

图论基础

无向加权图 $G = (V, E, W)$ ：

$V = {1, 2, \dots, n}$ ：顶点集（数据点）
$E$ ：边集
$W$ ：邻接权重矩阵， $W_{ij} \geq 0$

相似度图

将聚类问题转化为图分割问题：

def build_knn_graph(X, k=10, sigma=1.0):
    """
    构建k近邻相似度图
    
    使用ε-邻域或k-近邻方法
    """
    n = len(X)
    
    # 计算欧氏距离
    dist_sq = np.sum(X**2, axis=1, keepdims=True) + np.sum(X**2, axis=1) - 2*X@X.T
    
    # k-近邻
    W = np.zeros((n, n))
    for i in range(n):
        nearest = np.argsort(dist_sq[i])[1:k+1]  # 排除自身
        W[i, nearest] = np.exp(-dist_sq[i, nearest] / (2*sigma**2))
    
    # 对称化
    W = (W + W.T) / 2
    
    return W
 
def build_gaussian_graph(X, sigma=1.0):
    """全连接高斯相似度图"""
    dist_sq = np.sum(X**2, axis=1, keepdims=True) + np.sum(X**2, axis=1) - 2*X@X.T
    W = np.exp(-dist_sq / (2*sigma**2))
    np.fill_diagonal(W, 0)
    return W

图 Laplacian 矩阵

未归一化 Laplacian：

L = D - W

其中 $D$ 是度矩阵， $D_{ii} = \sum_{j} W_{ij}$ 。

归一化 Laplacian：

L_{sym} = I - D^{- 1/2} W D^{- 1/2}

或

L_{rw} = I - D^{- 1} W

Laplacian 的性质

def compute_laplacian(W, normalize='sym'):
    """
    计算图Laplacian矩阵
    """
    D = np.diag(np.sum(W, axis=1))
    L = D - W
    
    if normalize == 'sym':
        # 对称归一化 Laplacian
        D_inv_sqrt = np.diag(1.0 / np.sqrt(np.diag(D)))
        L_norm = D_inv_sqrt @ L @ D_inv_sqrt
        return L_norm
    elif normalize == 'rw':
        # 随机游走 Laplacian
        D_inv = np.diag(1.0 / np.diag(D))
        L_norm = D_inv @ L
        return L_norm
    else:
        return L
 
# 验证性质
def laplacian_properties(W, L):
    n = len(W)
    D = np.diag(np.sum(W, axis=1))
    
    print("Laplacian 性质验证：")
    print(f"1. 对称性: {np.allclose(L, L.T)}")
    print(f"2. 行和为0: {np.allclose(np.sum(L, axis=1), 0)}")
    
    # 特征值
    eigvals = np.linalg.eigvalsh(L)
    print(f"3. 特征值(前5个): {eigvals[:5]}")
    print(f"   最小特征值: {eigvals[0]:.6f} (应为0)")

谱聚类算法

未归一化谱聚类（Shi-Malik算法）：

输入：相似度矩阵 W，聚类数 k

1. 构建 Laplacian 矩阵 L = D - W
2. 计算 L 的前 k 个最小特征值对应的特征向量
3. 构造矩阵 U ∈ ℝ^{n×k}（特征向量按列）
4. 对 U 的行向量进行 k-means 聚类
5. 输出聚类结果

输出：每个数据点的聚类标签

def spectral_clustering(X, n_clusters, sigma=1.0, mode='knn'):
    """
    谱聚类算法
    
    Args:
        X: 数据矩阵 (n × d)
        n_clusters: 聚类数
        sigma: 高斯核参数
        mode: 'knn' 或 'full'
    """
    n = len(X)
    
    # 构建相似度图
    if mode == 'knn':
        W = build_knn_graph(X, k=10, sigma=sigma)
    else:
        W = build_gaussian_graph(X, sigma=sigma)
    
    # 计算 Laplacian
    D = np.diag(np.sum(W, axis=1))
    L = D - W
    
    # 标准化（可选）
    D_inv_sqrt = np.diag(1.0 / np.sqrt(np.diag(D) + 1e-10))
    L = D_inv_sqrt @ L @ D_inv_sqrt
    
    # 计算特征值和特征向量
    eigvals, eigvecs = np.linalg.eigh(L)
    
    # 选择前 k 个特征向量
    U = eigvecs[:, :n_clusters]
    
    # 对 U 的行进行 k-means
    from sklearn.cluster import KMeans
    
    # 归一化行向量（推荐）
    U_normalized = U / (np.linalg.norm(U, axis=1, keepdims=True) + 1e-10)
    
    kmeans = KMeans(n_clusters=n_clusters, random_state=42, n_init=10)
    labels = kmeans.fit_predict(U_normalized)
    
    return labels, eigvals[:n_clusters]
 
# 示例：非凸形状聚类
np.random.seed(42)
 
# 生成两个半月形数据
t1 = np.linspace(0, np.pi, 100)
moon1 = np.column_stack([np.cos(t1), np.sin(t1)]) + 0.1*np.random.randn(100, 2)
 
t2 = np.linspace(0, np.pi, 100)
moon2 = np.column_stack([1-np.cos(t2), 0.5-np.sin(t2)]) + 0.1*np.random.randn(100, 2)
 
X = np.vstack([moon1, moon2])
true_labels = np.array([0]*100 + [1]*100)
 
labels, eigvals = spectral_clustering(X, n_clusters=2, sigma=0.5)
accuracy = np.mean(labels == true_labels)
print(f"谱聚类准确率: {accuracy:.2%}")
print(f"前2个特征值: {eigvals}")

归一化 vs 未归一化谱聚类

类型	Laplacian	适用场景
Shi-Malik	$L_{rw}$	聚类大小相近
Ng-Jordan-Weiss	$L_{sym}$	聚类大小差异大
未归一化	$L$	一般情况

谱聚类的理论保证

定理（谱聚类收敛性）：如果数据点来自 $k$ 个分离良好的簇，则：

Laplacian 的前 $k$ 个特征值接近零
对应的特征向量趋近于指示向量
k-means 聚类正确率高

直观理解：

当簇内连接紧密、簇间连接稀疏时
$L$ 的谱结构反映图的连通性
小特征值对应低频模式（簇结构）