SVM算法实现及可视化边界绘制

import numpy as np
from utils import common
import matplotlib.pyplot as plt
from tqdm import tqdm
import scipy.spatial.distance as ssd

# 导入所需的库

def RBF_kernel(sigma):
    def kernel(x, y):
        '''
        :param x: N x d
        :param y: M x d
        :return: N x M
        '''
        dist = ssd.cdist(x, y, 'euclidean')  # 计算欧氏距离
        return np.exp(-np.square(dist) * sigma)  # 计算RBF核函数的值

    return kernel

# 定义RBF核函数

def auto_scale(X):
    X_var = np.var(X,ddof=1)  # 计算X的方差
    sigma = 1.0 / (X.shape[1] * X_var) if X_var != 0 else 1.0  # 计算自动缩放的参数sigma
    return sigma

# 自动缩放数据

def linear_kernel():
    def kernel(x, y):
        '''
        :param x: N x d
        :param y: M x d
        :return: N x M
        '''
        return x @ y.T  # 计算线性核函数的值

    return kernel

# 定义线性核函数

class SVM:
    def __init__(self, C, kernel, max_iter=500,show_fitting_bar=True):
        self.kernel = kernel  # 核函数
        self.C = C  # 惩罚参数
        self.max_iter = max_iter  # 最大迭代次数
        self.biary_encoder = common.BinaryEncoder()  # 二进制编码器
        self.tqdm = tqdm if show_fitting_bar else lambda x:x  # 进度条显示

    def error(self, index):
        pred = (self.alpha * self._label).T @ self.Q[:, index] + self.b  # 计算预测值
        return pred - self._label[index]  # 计算误差

    def fit(self, X: np.ndarray, y: np.ndarray, tol: float = 1e-3):
        '''
        :param X: N x d
        :param y: N
        :param tol: scalar
        :return:
        '''
        N, d = X.shape
        self._X = X  # 训练数据
        self._label = self.biary_encoder.fit_transform(y)  # 标签编码
        self.Q = self.kernel(X, X)  # 计算核矩阵
        self.alpha = np.zeros(N)  # 初始化拉格朗日乘子
        self.b = 0.0  # 初始化截距
        satisfied_count = 0

        for iter_count in self.tqdm(range(self.max_iter)):  # 迭代训练
            alpha_changed_count = 0
            for i in range(N):
                Ei = self.error(i)  # 计算第i个样本的误差
                alphai_old = self.alpha[i].copy()  # 保存旧的拉格朗日乘子
                yi = self._label[i]  # 第i个样本的标签

                # 选择违反KKT条件的alpha i
                if (yi * Ei < -tol and alphai_old < self.C) or (yi * Ei > tol and alphai_old > 0):
                    # 随机选择alpha j
                    j = np.random.randint(0, N)
                    while j == i:
                        j = np.random.randint(0, N)
                    yj = self._label[j]  # 第j个样本的标签
                    Ej = self.error(j)  # 计算第j个样本的误差
                    alphaj_old = self.alpha[j].copy()  # 保存旧的拉格朗日乘子

                    # 计算边界L和H
                    if yi == yj:
                        L = max(0.0, alphai_old + alphaj_old - self.C)
                        H = min(self.C, alphai_old + alphaj_old)
                    else:
                        L = max(0.0, alphaj_old - alphai_old)
                        H = min(self.C, self.C + alphaj_old - alphai_old)

                    # 计算eta
                    eta = 2.0 * self.Q[i, j] - self.Q[i, i] - self.Q[j, j]
                    if eta >= 0:
                        continue

                    # 更新alpha j
                    self.alpha[j] =alphaj_old - (yj * (Ei - Ej) / eta)

                    # 限制alpha j的取值范围
                    self.alpha[j] = np.clip(self.alpha[j], L, H)

                    # 更新alpha i
                    delta_j = alphaj_old - self.alpha[j]
                    if np.abs(delta_j) < tol:
                        continue
                    self.alpha[i] =alphai_old + yi * yj * delta_j

                    # 更新截距b
                    bi = self.b - Ei \
                         - yi * (self.alpha[i] - alphai_old) * self.Q[i, i] \
                         - yj * (self.alpha[j] - alphaj_old) * self.Q[i, j]
                    bj = self.b - Ej \
                         - yi * (self.alpha[i] - alphai_old) * self.Q[i, j] \
                         - yj * (self.alpha[j] - alphaj_old) * self.Q[j, j]
                    if 0.0+tol < self.alpha[i] < self.C-tol:
                        self.b = bi
                    elif 0.0+tol < self.alpha[j] < self.C-tol:
                        self.b = bj
                    else:
                        self.b = (bi + bj) / 2.0

                    alpha_changed_count += 1
            if alpha_changed_count == 0:
                satisfied_count += 1
            else:
                satisfied_count = 0
            if satisfied_count >= 5:
                break

        mask = self.alpha > tol  # 选取大于阈值的拉格朗日乘子
        self.support_vectors = mask  # 支持向量
        self.X_sup = X[mask]  # 支持向量的特征
        self.Y_sup = self._label[mask]  # 支持向量的标签
        self.alpha_sup = self.alpha[mask]  # 支持向量的拉格朗日乘子

        del self._X
        del self._label
        del self.alpha
        del self.Q

    def predict_proba(self, X):
        sigma = self.kernel(X, self.X_sup)  # 计算测试样本与支持向量的核矩阵
        sigma = self.Y_sup[np.newaxis, :] * self.alpha_sup[np.newaxis, :] * sigma  # 计算预测值
        y_hat = np.sum(sigma, axis=1) + self.b  # 计算预测概率
        return y_hat

    def predict(self, X):
        r = self.predict_proba(X)  # 预测概率
        r[r>0]=1
        r[r<=0]=-1
        r = self.biary_encoder.inverse_transform(r)  # 将预测概率转换为标签
        return r

# 定义支持向量机类

def plot_boundary():

    def generate_data(n_sample=50):
        x1 = np.random.rand(n_sample, 1) * 6 - 2
        y1 = x1 + np.random.normal(0, 0.7, (n_sample, 1))
        x2 = np.random.rand(n_sample, 1) * 6 - 3
        y2 = x2 + np.random.normal(0, 0.7, (n_sample, 1)) + 3
        c1 = np.hstack([x1, y1, np.ones_like(x1) * -1])
        c2 = np.hstack([x2, y2, np.ones_like(x2)])
        data_with_label = np.vstack([c1, c2])
        data, label = data_with_label[:, :-1], data_with_label[:, -1]

        return data, label

    import matplotlib.pyplot as plt
    np.random.seed(1)
    X,y = generate_data(50)
    kernels = [RBF_kernel(1.0),linear_kernel()]
    plt.figure(figsize=(10,5))
    for i,kernel in enumerate(kernels):
        plt.subplot(1,2,i+1)
        model = SVM(C=1, kernel=kernel,max_iter=1000)
        model.fit(X, y, tol=1e-5)
        sv = model.support_vectors

        plt.scatter(X[sv, 0], X[sv, 1], c=y[sv], marker='*',\n                    linewidths=0.5, edgecolors=(0, 0, 0, 1))
        plt.scatter(X[~sv, 0], X[~sv, 1], c=y[~sv],\n                    linewidths=0.5, edgecolors=(0, 0, 0, 1))
        xvals = np.linspace(-3, 6, 200)
        yvals = np.linspace(-3, 6, 200)
        xx, yy = np.meshgrid(xvals, yvals)

        pred = model.predict_proba(np.c_[xx.ravel(), yy.ravel()])
        zz = np.reshape(pred, xx.shape)
        plt.pcolormesh(xx, yy, zz, zorder=0)
        plt.contour(xx, yy, zz, levels=(-1, 0, 1), colors='w', linewidths=1.5, zorder=1, linestyles='solid')

        plt.xlim([-3, 6])
        plt.ylim([-3, 6])
    plt.show()


if __name__ == '__main__':
    plot_boundary()
    from sklearn.datasets import load_breast_cancer
    data, label = load_breast_cancer(return_X_y=True)

    sigma = auto_scale(data)
    svm = SVM(1.0, kernel=RBF_kernel(sigma))
    svm.fit(data, label)
    y_pred = svm.predict(data)
    acc = np.sum(y_pred == label) / len(y_pred)
    print('SVM Test Acc %.4f' % acc)

代码解释

导入所需库：

• numpy：用于数值计算
• utils.common：包含一些常用的工具函数，例如二进制编码器
• matplotlib.pyplot：用于绘制图形
• tqdm：用于显示进度条
• scipy.spatial.distance：用于计算距离

RBF核函数：

• RBF_kernel(sigma)：定义RBF核函数，接受参数sigma，返回一个函数kernel。
• kernel(x, y)：计算两个向量x和y之间的RBF核函数值。
  • ssd.cdist(x, y, 'euclidean')：计算x和y之间的欧氏距离
  • np.exp(-np.square(dist) * sigma)：根据欧氏距离和sigma计算RBF核函数值

自动缩放数据：

• auto_scale(X)：自动缩放数据，接受参数X，返回一个sigma值。
  • np.var(X,ddof=1)：计算X的方差
  • sigma = 1.0 / (X.shape[1] * X_var) if X_var != 0 else 1.0：根据方差计算sigma值，避免出现除零错误

线性核函数：

• linear_kernel()：定义线性核函数，返回一个函数kernel。
• kernel(x, y)：计算两个向量x和y之间的线性核函数值，即内积x @ y.T。

支持向量机类：

• SVM：定义支持向量机类，接受参数C、kernel、max_iter、show_fitting_bar，包含以下方法：
  • __init__(self, C, kernel, max_iter=500,show_fitting_bar=True)：初始化SVM类，设置参数
  • error(self, index)：计算第index个样本的误差
  • fit(self, X: np.ndarray, y: np.ndarray, tol: float = 1e-3)：训练SVM模型，接受训练数据X、标签y和容差tol，返回训练好的模型
    • 迭代训练：使用tqdm显示进度条
    • 选择违反KKT条件的alpha i：根据KKT条件选择需要更新的拉格朗日乘子
    • 随机选择alpha j：随机选择另一个拉格朗日乘子进行更新
    • 计算边界L和H：限制alpha j的取值范围
    • 计算eta：用于更新alpha j
    • 更新alpha j：根据eta和误差更新alpha j
    • 限制alpha j的取值范围：将alpha j限制在L和H之间
    • 更新alpha i：根据alpha j的更新量更新alpha i
    • 更新截距b：根据更新后的alpha和误差更新截距
    • 选取支持向量：将大于阈值的拉格朗日乘子对应的样本标记为支持向量
  • predict_proba(self, X)：预测测试样本的概率，接受测试数据X，返回预测概率
  • predict(self, X)：预测测试样本的标签，接受测试数据X，返回预测标签

绘制决策边界：

• plot_boundary()：绘制决策边界，包含以下步骤：
  • 生成数据：使用generate_data()生成测试数据
  • 创建模型：使用SVM类创建模型
  • 训练模型：使用fit()方法训练模型
  • 绘制支持向量：使用plt.scatter()绘制支持向量
  • 绘制其他样本：使用plt.scatter()绘制其他样本
  • 计算预测概率：使用predict_proba()方法计算测试数据的预测概率
  • 绘制等高线：使用plt.contour()绘制等高线
  • 绘制填充区域：使用plt.pcolormesh()绘制填充区域

测试SVM模型：

• 使用load_breast_cancer()加载乳腺癌数据集
• 使用auto_scale()计算sigma值
• 创建SVM模型并训练
• 使用predict()方法预测测试数据
• 计算准确率并打印

总结

本代码实现了一个简单的SVM算法，包含RBF核函数和线性核函数，并可视化决策边界。代码中包含详细的注释，解释了每行代码的作用。你可以根据需要修改参数，例如惩罚参数C、核函数类型、迭代次数等，来调整模型的性能。

注意：

• 此代码仅供参考，可能需要根据实际情况进行调整。
• 由于代码较为复杂，建议你仔细阅读注释，理解代码的逻辑。
• 如果你对SVM算法有更深入的了解，可以尝试使用其他优化方法，例如SMO算法，来提高模型的效率。