当前位置：首页 > news >正文

支持向量机（SVM）python语言版本

news 2025/6/25 12:57:22

这个版本是AI写的，根据其他语言改编而来，跑了一下可以运行，但没有逐语句检查过：
可以通过gitee下载：】
git clone https://gitee.com/iamohenry/svm.git

import numpy as np
import matplotlib.pyplot as plt
from time import time

class SVM:
    def __init__(self, C=100, sigma=0.25, eps=1e-3, max_iter=100):
        self.C = C          # 惩罚系数
        self.sigma = sigma # 高斯核参数
        self.eps = eps      # 停止阈值
        self.max_iter = max_iter

    def fit(self, X, y):
        self.X = X
        self.y = y
        self.n_samples = X.shape[0]

        # 初始化参数
        self.a = np.zeros(self.n_samples)
        self.b = 0.0
        self.E = -y.copy() # E = f(x) - y，初始时f(x)=0

        # 预计算核矩阵（关键修正点）
        self.K = np.zeros((self.n_samples, self.n_samples))
        for i in range(self.n_samples):
            for j in range(i, self.n_samples):
                self.K[i,j] = self.kernel(X[i], X[j])
                self.K[j,i] = self.K[i,j]

        # SMO主循环
        for iter in range(self.max_iter):
            # 选择alpha对
            i1, i2 = self.select_alpha_pair()
            if i1 is None:
                break

            # 优化alpha对
            self.optimize_alpha_pair(i1, i2)

            # 提前停止检查
            if (iter+1) % 10 == 0:
                if self.check_stopping_condition():
                    break

        # 提取支持向量
        self.extract_support_vectors()
        print(f'找到 {len(self.sv_indices)} 个支持向量')

    def kernel(self, x1, x2):
        return np.exp(-np.sum((x1-x2)**2) / (2*self.sigma**2))

    def select_alpha_pair(self):
        # 第一层选择：违反KKT条件的样本
        violate_max = 0
        i1 = None

        # 遍历所有样本
        for i in range(self.n_samples):
            ai = self.a[i]
            yi = self.y[i]
            Ei = self.E[i]

            # KKT条件判断
            if (ai < self.C and yi*Ei < -self.eps) or (ai > 0 and yi*Ei > self.eps):
                # 计算违反程度
                violate = abs(yi*Ei)
                if violate > violate_max:
                    violate_max = violate
                    i1 = i

        if i1 is None:
            return None, None

        # 选择第二个alpha：最大|E1-E2|
        i2 = np.argmax(np.abs(self.E - self.E[i1]))
        if i2 == i1:
            i2 = (i1 + 1) % self.n_samples
        return i1, i2

    def optimize_alpha_pair(self, i1, i2):
        if i1 == i2:
            return

        # 获取相关参数
        y1, y2 = self.y[i1], self.y[i2]
        a1_old, a2_old = self.a[i1], self.a[i2]
        E1, E2 = self.E[i1], self.E[i2]
        K11 = self.K[i1,i1]
        K12 = self.K[i1,i2]
        K22 = self.K[i2,i2]

        # 计算上下界
        if y1 != y2:
            L = max(0, a2_old - a1_old)
            H = min(self.C, self.C + a2_old - a1_old)
        else:
            L = max(0, a2_old + a1_old - self.C)
            H = min(self.C, a2_old + a1_old)
        if L >= H:
            return

        # 计算eta
        eta = K11 + K22 - 2*K12
        if eta <= 0:
            return

        # 更新a2
        a2_new = a2_old + y2*(E1 - E2)/eta
        a2_new = np.clip(a2_new, L, H)

        # 更新a1
        a1_new = a1_old + y1*y2*(a2_old - a2_new)

        # 更新偏置b
        b1 = E1 + y1*(a1_new - a1_old)*K11 + y2*(a2_new - a2_old)*K12 + self.b
        b2 = E2 + y1*(a1_new - a1_old)*K12 + y2*(a2_new - a2_old)*K22 + self.b
        if 0 < a1_new < self.C:
            new_b = b1
        elif 0 < a2_new < self.C:
            new_b = b2
        else:
            new_b = (b1 + b2)/2

        # 更新参数
        self.a[i1] = a1_new
        self.a[i2] = a2_new
        self.b = new_b

        # 更新误差缓存（关键修正点）
        for i in range(self.n_samples):
            self.E[i] = self._fx(i) - self.y[i]

    def _fx(self, i):
        return np.dot(self.a * self.y, self.K[:,i]) - self.b

    def check_stopping_condition(self):
        # 检查对偶间隙
        dual_gap = np.sum(self.a) - 0.5*np.sum(np.outer(self.a*self.y, self.a*self.y)*self.K)
        if dual_gap < self.eps:
            return True
        return False

    def extract_support_vectors(self):
        self.sv_indices = np.where(self.a > 1e-5)[0]
        self.sv_X = self.X[self.sv_indices]
        self.sv_y = self.y[self.sv_indices]
        self.sv_a = self.a[self.sv_indices]

    def predict(self, X):
        n = X.shape[0]
        y_pred = np.zeros(n)
        for i in range(n):
            k = np.array([self.kernel(X[i], sv) for sv in self.sv_X])
            y_pred[i] = np.sign(np.dot(self.sv_a * self.sv_y, k) - self.b)
        return y_pred

# 数据生成函数
def generate_moon_data(n_train=300, n_val=2000):
    r, wth, dis = 10, 6, -6.5
    total = n_train + n_val
    X = np.zeros((total, 2))
    y = np.zeros(total)

    for i in range(total):
        if np.random.rand() > 0.5:
            theta = np.pi * np.random.rand()
            length = r + wth*(np.random.rand()-0.5)
            X[i] = [length*np.cos(theta), length*np.sin(theta)]
            y[i] = -1
        else:
            beta = -np.pi * np.random.rand()
            length = r + wth*(np.random.rand()-0.5)
            X[i] = [length*np.cos(beta)+r, length*np.sin(beta)-dis]
            y[i] = 1
    return X[:n_train], y[:n_train], X[n_train:], y[n_train:]

# 数据预处理
def normalize(X_train, X_val):
    # 统一归一化参数
    miu = np.mean(np.vstack([X_train, X_val]), axis=0)
    max_val = np.max(np.abs(np.vstack([X_train, X_val])), axis=0)

    X_train_norm = (X_train - miu) / max_val
    X_val_norm = (X_val - miu) / max_val
    return X_train_norm, X_val_norm

# 可视化函数
def visualize(model, X, y):
    plt.figure(figsize=(10,8))

    # 绘制训练数据
    pos = y == 1
    neg = y == -1
    plt.scatter(X[pos,0], X[pos,1], c='b', label='Class 1', alpha=0.5)
    plt.scatter(X[neg,0], X[neg,1], c='g', label='Class -1', alpha=0.5)

    if len(model.sv_indices) > 0:
        # 绘制支持向量
        plt.scatter(model.sv_X[:,0], model.sv_X[:,1], s=80,
                    facecolors='none', edgecolors='r', linewidths=2, label='Support Vectors')

        # 绘制决策边界
        x_min, x_max = X[:,0].min()-0.1, X[:,0].max()+0.1
        y_min, y_max = X[:,1].min()-0.1, X[:,1].max()+0.1
        xx, yy = np.meshgrid(np.linspace(x_min, x_max, 200),
                             np.linspace(y_min, y_max, 200))
        Z = model.predict(np.c_[xx.ravel(), yy.ravel()]).reshape(xx.shape)

        # 仅在有支持向量时绘制
        plt.contour(xx, yy, Z, colors='r', levels=[0], linewidths=2)

    plt.title('SVM Classification Result')
    plt.xlabel('Feature 1')
    plt.ylabel('Feature 2')
    plt.legend()
    plt.show()

if __name__ == "__main__":
    # 生成数据
    X_train, y_train, X_val, y_val = generate_moon_data(n_train=300, n_val=2000)

    # 数据归一化
    X_train, X_val = normalize(X_train, X_val)

    # 训练模型（调整参数）
    start = time()
    svm = SVM(C=10, sigma=0.2, max_iter=200) # 调整参数
    svm.fit(X_train, y_train)

    # 验证模型
    y_pred = svm.predict(X_val)
    error_rate = np.mean(y_pred != y_val)

    print(f'运行时间: {time()-start:.2f}秒')
    print(f'验证集错误率: {error_rate*100:.2f}%')

    # 可视化结果
    visualize(svm, X_train, y_train)