使用 NumPy 和 SMO 算法实现 SVM，并使用稀疏矩阵存储 Gram 矩阵

使用 NumPy 和 SMO 算法实现 SVM，并使用稀疏矩阵存储 Gram 矩阵

本代码示例展示了如何使用 NumPy 和 SMO 算法来实现 SVM，并使用稀疏矩阵存储 Gram 矩阵，以减少内存使用。代码中使用了 MNIST 数据集进行测试和训练，并提供了详细的注释。

import numpy as np
from scipy.sparse import csr_matrix
from sklearn.metrics import accuracy_score

class SVM:
    def __init__(self, C=1.0, tol=0.001, max_iter=100):
        self.C = C
        self.tol = tol
        self.max_iter = max_iter

    def fit(self, X, y):
        n_samples, n_features = X.shape

        # 初始化参数
        self.alpha = np.zeros(n_samples)
        self.b = 0
        self.errors = np.zeros(n_samples)

        # 保存训练数据
        self.X = csr_matrix(X)
        self.y = y

        # 计算 Gram 矩阵
        gram = self._kernel(self.X)

        # 迭代优化参数
        num_changed_alphas = 0
        iter_count = 0
        while iter_count < self.max_iter and num_changed_alphas > 0:
            num_changed_alphas = 0
            for i in range(n_samples):
                E_i = self._decision_function(gram, i) - y[i]
                if (y[i] * E_i < -self.tol and self.alpha[i] < self.C) or \
                        (y[i] * E_i > self.tol and self.alpha[i] > 0):

                    j = self._select_second_alpha(i, n_samples)
                    E_j = self._decision_function(gram, j) - y[j]

                    alpha_i_old = self.alpha[i]
                    alpha_j_old = self.alpha[j]

                    if y[i] != y[j]:
                        L = max(0, self.alpha[j] - self.alpha[i])
                        H = min(self.C, self.C + self.alpha[j] - self.alpha[i])
                    else:
                        L = max(0, self.alpha[i] + self.alpha[j] - self.C)
                        H = min(self.C, self.alpha[i] + self.alpha[j])

                    if L == H:
                        continue

                    eta = 2 * gram[i, j] - gram[i, i] - gram[j, j]
                    if eta >= 0:
                        continue

                    self.alpha[j] -= (y[j] * (E_i - E_j)) / eta
                    self.alpha[j] = np.clip(self.alpha[j], L, H)

                    if abs(self.alpha[j] - alpha_j_old) < 1e-5:
                        continue

                    self.alpha[i] += y[i] * y[j] * (alpha_j_old - self.alpha[j])

                    self.errors[i] = self._decision_function(gram, i) - y[i]
                    self.errors[j] = self._decision_function(gram, j) - y[j]

                    num_changed_alphas += 1

            iter_count += 1

    def predict(self, X):
        decisions = self._decision_function(self._kernel(self.X), X)
        return np.sign(decisions)

    def _decision_function(self, gram, X):
        decisions = np.dot(self.alpha * self.y, gram.T) + self.b
        return decisions

    def _select_second_alpha(self, first_alpha, n_samples):
        second_alpha = first_alpha
        while second_alpha == first_alpha:
            second_alpha = np.random.randint(n_samples)
        return second_alpha

    def _kernel(self, X1, X2=None):
        if X2 is None:
            X2 = X1
        return X1.dot(X2.T)


# 加载 MNIST 数据集
train_images = 'train-images-idx3-ubyte'  # 训练图像数据文件路径
train_labels = 'train-labels-idx1-ubyte'  # 训练标签数据文件路径
test_images = 't10k-images-idx3-ubyte'  # 测试图像数据文件路径
test_labels = 't10k-labels-idx1-ubyte'  # 测试标签数据文件路径

# 加载训练数据集
with open(train_images, 'rb') as f:
    train_images_data = np.frombuffer(f.read(), dtype=np.uint8, offset=16).reshape(-1, 28*28) / 255.0
with open(train_labels, 'rb') as f:
    train_labels_data = np.frombuffer(f.read(), dtype=np.uint8, offset=8)

# 加载测试数据集
with open(test_images, 'rb') as f:
    test_images_data = np.frombuffer(f.read(), dtype=np.uint8, offset=16).reshape(-1, 28*28) / 255.0
with open(test_labels, 'rb') as f:
    test_labels_data = np.frombuffer(f.read(), dtype=np.uint8, offset=8)

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_images_data, test_images_data, train_labels_data, test_labels_data

# 创建 SVM 分类器并进行训练
svm = SVM()
svm.fit(X_train, y_train)

# 在训练集上进行预测
y_pred_train = svm.predict(X_train)

# 在测试集上进行预测
y_pred_test = svm.predict(X_test)

# 计算训练集和测试集的准确率
accuracy_train = accuracy_score(y_train, y_pred_train)
accuracy_test = accuracy_score(y_test, y_pred_test)

print('训练集准确率:', accuracy_train)
print('测试集准确率:', accuracy_test)

代码说明：

1. SVM 类

   • __init__ 方法：初始化 SVM 模型的参数，包括正则化参数 C、容忍度 tol 和最大迭代次数 max_iter。
   • fit 方法：训练 SVM 模型，包括初始化参数、计算 Gram 矩阵、迭代优化参数。
   • predict 方法：使用训练好的模型进行预测。
   • _decision_function 方法：计算决策函数，用于判断样本所属类别。
   • _select_second_alpha 方法：选择第二个 alpha 值进行优化。
   • _kernel 方法：计算核函数，这里使用线性核函数。

2. 稀疏矩阵存储 Gram 矩阵

   • 使用 scipy.sparse.csr_matrix 创建稀疏矩阵，以减少存储 Gram 矩阵所需的内存。

3. 加载 MNIST 数据集

   • 从本地文件加载 MNIST 训练集和测试集数据。

4. 训练和评估

   • 创建 SVM 对象并使用 fit 方法进行训练。
   • 使用 predict 方法对训练集和测试集进行预测。
   • 使用 accuracy_score 计算训练集和测试集的准确率。

注意：

• 确保 MNIST 数据集文件 train-images-idx3-ubyte、train-labels-idx1-ubyte、t10k-images-idx3-ubyte、t10k-labels-idx1-ubyte 位于当前目录下。
• 可以根据需要调整 SVM 类的参数，例如正则化参数 C 和容忍度 tol。
• 可以尝试使用其他核函数，例如 RBF 核函数，以提高模型的性能。

本代码示例提供了一个简单的 SVM 实现，使用稀疏矩阵存储 Gram 矩阵，以减少内存使用。希望这能帮助你理解 SVM 的实现过程。