在jupyternotebook上使用Python的cvxopt模块求解硬间隔的SVM,测试数据选择mnist数据集以及北京理工大学的手写数据集做泛化能力测试。对北理工的手写数据集进行预测时有进行预处理和不进行预处理两种情况。
时间: 2024-05-14 18:14:12 浏览: 230
首先,需要安装cvxopt模块。可以通过以下命令进行安装:
```
!pip install cvxopt
```
接下来,需要加载数据集。可以使用sklearn中的load_digits函数加载mnist数据集,代码如下:
```python
from sklearn.datasets import load_digits
digits = load_digits()
X = digits.data
y = digits.target
```
对于北理工的手写数据集,可以使用PIL库中的Image模块来读取图片,并使用numpy库将图片转换为数组。代码如下:
```python
from PIL import Image
import numpy as np
def read_image(file_path):
img = Image.open(file_path).convert('L')
img = img.resize((8, 8))
img_array = np.array(img)
img_array[img_array > 0] = 1
return img_array.reshape(1, -1)
X_test = read_image('test.png')
```
其中,'test.png'是北理工的手写数据集中的一张图片,需要将其放置在当前工作目录下。
接下来,可以使用cvxopt中的qp函数求解SVM的对偶问题。对于硬间隔SVM,其对偶问题可以表示为:
$$\max\limits_{\alpha}\sum_{i=1}^{n}\alpha_i-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}y_iy_j\alpha_i\alpha_j\langle x_i,x_j\rangle$$
$$s.t.\ \sum_{i=1}^{n}y_i\alpha_i=0,\ \alpha_i\geq 0$$
其中,$\alpha$为拉格朗日乘子向量,$n$为样本数量,$x_i$为第$i$个样本的特征向量,$y_i\in\{-1,1\}$为第$i$个样本的标签。
代码如下:
```python
from cvxopt import matrix, solvers
# 构造P矩阵
n_samples, n_features = X.shape
K = np.dot(X, X.T)
P = matrix(np.outer(y, y) * K)
# 构造q、G、h、A、b向量
q = matrix(-np.ones((n_samples, 1)))
G = matrix(np.vstack((-np.eye(n_samples), np.eye(n_samples))))
h = matrix(np.hstack((np.zeros(n_samples), np.ones(n_samples) * C)))
A = matrix(y.reshape(1, -1))
b = matrix(np.zeros(1))
# 求解QP问题
sol = solvers.qp(P, q, G, h, A, b)
alphas = np.array(sol['x'])
# 计算权重向量w和截距b
w = np.dot(X.T, y * alphas)
b = y[alphas > 1e-5][0] - np.dot(X[alphas > 1e-5], w)
```
在求解QP问题之前,需要先计算样本之间的内积矩阵$K$,并构造P矩阵。同时,还需要构造G、h、A、b向量,用于表示拉格朗日乘子的约束条件。求解QP问题可以使用cvxopt中的qp函数。
接下来,可以使用求解得到的权重向量和截距,对测试数据集进行预测。代码如下:
```python
y_pred = np.sign(np.dot(X_test, w) + b)
```
对于mnist数据集,可以使用交叉验证来评估模型的性能。代码如下:
```python
from sklearn.model_selection import KFold
kf = KFold(n_splits=5, shuffle=True, random_state=42)
accs = []
for train_index, test_index in kf.split(X):
X_train, y_train = X[train_index], y[train_index]
X_test, y_test = X[test_index], y[test_index]
# 训练模型
...
# 预测测试集
y_pred = np.sign(np.dot(X_test, w) + b)
# 计算准确率
acc = np.mean(y_pred == y_test)
accs.append(acc)
print('Accuracy: %.2f%% (+/- %.2f%%)' % (np.mean(accs) * 100, np.std(accs) * 100))
```
对于北理工的手写数据集,可以分别对预处理和未预处理的数据进行预测,并比较两种情况的性能。代码如下:
```python
# 未预处理
X_test = read_image('test.png')
y_pred = np.sign(np.dot(X_test, w) + b)
print('Prediction: %d' % y_pred[0])
# 预处理
X_test = read_image('test.png')
X_test_processed = np.dot(X_test - np.mean(X_train, axis=0), pca.components_.T)
y_pred = np.sign(np.dot(X_test_processed, w_pca) + b_pca)
print('Prediction: %d' % y_pred[0])
```
其中,未预处理的代码与之前相同,预处理的代码需要先对测试数据进行降维,再进行预测。在这里,使用PCA降维,代码如下:
```python
from sklearn.decomposition import PCA
# 对训练数据进行降维
pca = PCA(n_components=32)
X_train_processed = np.dot(X_train - np.mean(X_train, axis=0), pca.components_.T)
# 训练模型
...
# 对测试数据进行降维
X_test_processed = np.dot(X_test - np.mean(X_train, axis=0), pca.components_.T)
# 预测测试集
y_pred = np.sign(np.dot(X_test_processed, w) + b)
# 计算准确率
acc = np.mean(y_pred == y_test)
print('Accuracy: %.2f%%' % (acc * 100))
```
需要注意的是,在对测试数据进行降维时,需要使用与训练数据相同的PCA模型。因此,在训练模型时需要保存PCA模型,代码如下:
```python
# 对训练数据进行降维
pca = PCA(n_components=32)
X_train_processed = np.dot(X_train - np.mean(X_train, axis=0), pca.components_.T)
# 训练模型
...
# 保存PCA模型
w_pca = np.dot(pca.components_.T, w)
b_pca = b + np.dot(np.mean(X_train, axis=0), w_pca)
# 对测试数据进行降维
X_test_processed = np.dot(X_test - np.mean(X_train, axis=0), pca.components_.T)
# 预测测试集
y_pred = np.sign(np.dot(X_test_processed, w_pca) + b_pca)
# 计算准确率
acc = np.mean(y_pred == y_test)
print('Accuracy: %.2f%%' % (acc * 100))
```
阅读全文