Python实现硬间隔SVM、软间隔SVM

SVM：从类别上理解可以将SVM分为硬间隔SVM(hard-margin SVM)、软间隔SVM(soft-margin SVM)、核SVM。

个人理解：在数据线性可分的前提下，硬间隔SVM是找到离分类平面较近的支持向量，再由支持向量找到最优超平面将数据进行分类。软间隔SVM是为了在线性不可分的数据中适用，对每个样本点引入一个松弛变量，即在约束条件中增加一个惩罚项。核技巧能够让svm从普通的特征空间映射到高纬空间，从而使得线性不可分的数据在高纬空间中可分。

硬间隔SVM可以转化为一个凸优化问题，并调用cvxopt库进行求解。

在做手写数字识别时，以784维度进行计算时，提示错误，说存在多余约束。。尝试多种方法无果后，决定降低数据的维度 使用图像的方向梯度直方图（hog）作为特征向量，将784纬的特征降至64纬，就可以了，对于凸优化问题，由于我的数学能力有限，理解的不是很透彻。

数据集在模板匹配算法中给出。

参考资料：

(145条消息) 手推支持向量机01-硬间隔SVM-模型定义（最大间隔分类器）_敲代码的乔帮主的博客-CSDN博客

import os

import cv2
import cvxopt
import numpy as np
import scipy.io as sio
from cvxopt import solvers

'''
准备train test数据和标签
训练集400张来自minist_train，第一测试集50张来自minist_test，第二测试集50张来自手写数字
'''

#转换train，test的minist数据
train_data = sio.loadmat('./课程数据集/mat格式的MNIST数据/train_images.mat')
key_train  = list(train_data.keys())[-1]
train      = train_data[key_train]

train_data2 = sio.loadmat('./课程数据集/mat格式的MNIST数据/train_labels.mat')
key_label   = list(train_data2.keys())[-1]
train_label = train_data2[key_label]

test_data  = sio.loadmat('./课程数据集/mat格式的MNIST数据/test_images.mat')
key_test   = list(test_data.keys())[-1]
test       = test_data[key_test]

test_data2 = sio.loadmat('./课程数据集/mat格式的MNIST数据/test_labels.mat')
test_key   = list(test_data2.keys())[-1]
test_label = test_data2[test_key]

#从数据中选取出[6,4],并将其转换为[1,-1]

train_num  = 400
test_num   = 50

j          = 0
q          = 0
number     = [6,4]

#从train中选取数据
train_img  = np.zeros((28,28,400),dtype=np.uint8)
train_lab  = np.zeros((1,400),dtype=np.float64)

while True:
    for i in range(60000):
        if train_label[:,i] == number[0] or train_label[:,i] == number[1]:
            train_img[:,:,j] = train[:,:,i]
            train_lab[:,j]   = train_label[:,i]
            if train_lab[:,j] == number[0]:
                train_lab[:,j] = 1
            elif train_lab[:,j] == number[1]:
                train_lab[:,j] = -1
            j +=1
            if j >= 400:
                break
    if j >=400:
        break

#从test中选取数据

test_img  = np.zeros((28,28,50),dtype=np.uint8)
test_lab  = np.zeros((1,50),dtype=np.float64)

while True:
    for i in range(10000):
        if test_label[:,i] == number[0] or test_label[:,i]==number[1]:
            test_img[:,:,q] = test[:,:,i]
            test_lab[:,q]   = test_label[:,i]
            if test_lab[:,q] == number[0]:
                test_lab[:,q] = 1
            elif test_lab[:,q] == number[1]:
                test_lab[:,q] = -1
            q += 1
            if q >=50:
                break
    if q >=50:
        break

#从数字中提取特征
SZ=20
bin_n = 16
affine_flags = cv2.WARP_INVERSE_MAP|cv2.INTER_LINEAR
#使用二阶矩矫正图像
def deskew(img):
    m = cv2.moments(img)
    if abs(m['mu02']) < 1e-2:
        return img.copy()
    skew = m['mu11']/m['mu02']
    M = np.float32([[1, skew, -0.5*SZ*skew], [0, 1, 0]])
    img = cv2.warpAffine(img,M,(SZ, SZ),flags=affine_flags)
    return img
#提取hog特征
def hog(img):
    gx = cv2.Sobel(img, cv2.CV_32F, 1, 0)
    gy = cv2.Sobel(img, cv2.CV_32F, 0, 1)
    mag, ang = cv2.cartToPolar(gx, gy)
    bins = np.int32(bin_n*ang/(2*np.pi))
    bin_cells = bins[:10,:10], bins[10:,:10], bins[:10,10:], bins[10:,10:]
    mag_cells = mag[:10,:10], mag[10:,:10], mag[:10,10:], mag[10:,10:]
    hists = [np.bincount(b.ravel(), m.ravel(), bin_n) for b, m in zip(bin_cells, mag_cells)]
    hist = np.hstack(hists)
    return hist

#训练数据和测试数据
new_train = np.zeros((400,64),dtype=np.float64)
new_test  = np.zeros((50,64),dtype=np.float64)
train_lab = train_lab.transpose().flatten()
test_lab  = test_lab.transpose().flatten()

for i in range(400):
    img = train_img[:,:,i]
    out_des = deskew(img)
    out_hog = hog(out_des)
    # print(out_hog.shape)
    new_train[i,:] = out_hog
for j in range(50):
    img2 = test_img[:,:,j]
    out_des2 = deskew(img2)
    out_hog2 = hog(out_des2)
    new_test[j,:] = out_hog2

write_num_test = np.zeros((50,64),dtype=np.float64)
write_num_label = np.zeros((50,1),dtype=np.float32)

for i in range(1,26):
    path = './课程数据集/手写数字/{}/{}.bmp'.format(str(4),i)
    img  = cv2.imdecode(np.fromfile(path,np.uint8),cv2.IMREAD_GRAYSCALE)
    hog_img    = hog(img)
    write_num_test[i-1,:] = hog_img

    write_num_label[i-1,:] = -1
for j in range(26,51):
    path_2 = './课程数据集/手写数字/{}/{}.bmp'.format(str(6),j)
    img_2  = cv2.imdecode(np.fromfile(path_2,np.uint8),cv2.IMREAD_GRAYSCALE)
    hog_img_2  = hog(img_2)
    write_num_test[j-1,:] = hog_img_2
    write_num_label[j-1,:] = 1
write_num_label = write_num_label.transpose().flatten()



def linear_kernel(x1, x2):
    return np.dot(x1, x2)


class mySVM():
    def __init__(self,kernel=linear_kernel,C=None):
        self.kernel = kernel
        self.C      = C

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

        # 需要求的参数  P q G h A b

        K = np.zeros((n_samples, n_samples))
        for i in range(n_samples):
            for j in range(n_samples):
                K[i][j] = self.kernel(X[i], X[j])
        P = cvxopt.matrix(np.outer(y, y) * K)   #np.outer求外积

        q = cvxopt.matrix(np.ones(n_samples) * -1)  # 为列向量

        A = cvxopt.matrix(y, (1, n_samples))  # y的转置 为行向量

        b = cvxopt.matrix(0.0)
        if self.C is None:
            G = cvxopt.matrix(np.diag(np.ones(n_samples) * -1))

            h = cvxopt.matrix(np.zeros(n_samples))
        else:
            # 软间隔与硬间隔的区别在于对每个样本点引入了一个松弛变量
            arg1 = np.diag(np.ones(n_samples) * -1)
            arg2 = np.diag(np.ones(n_samples))
            G = cvxopt.matrix(np.vstack((arg1, arg2)))
            arg1 = np.zeros(n_samples)
            arg2 = np.ones(n_samples) * self.C
            h = cvxopt.matrix(np.hstack((arg1, arg2)))

        solution = solvers.qp(P, q, G, h, A, b)

        a = np.ravel(solution['x'])
        
        if self.C is None:
            
            sv = a > 1e-9
        else:
            sv = (a > 1e-9) * (a < self.C)

        index = np.arange(len(a))[sv]  # a>1e-9  返回true or false的数组   取出对应为true的下标

        self.a = a[sv]  # 大于0对应的a的值index
        self.sv_X = X[sv]  # 支持向量的X index
        self.sv_y = y[sv]  # 支持向量的y  index

        print('%d 点中有 %d 个支持向量' % (n_samples, len(self.a)))

        # 求w
        self.w = np.zeros(n_feature)
        for i in range(len(self.a)):
            self.w += self.a[i] * self.sv_y[i] * self.sv_X[i]
        # 求b
        self.b = 0
        for i in range(len(self.a)):
            self.b += self.sv_y[i]
            self.b -= np.sum(self.a * self.sv_y * K[index[i]][index])
        self.b /= len(self.a)

    def predict(self, X_test):
        return np.sign(np.dot(X_test, self.w) + self.b)

def test(type,X_train,X_test,train_lab,test_lab):
    clf = mySVM()
    clf.fit(X_train,train_lab)

    test_pridict = clf.predict(X_test)
    correct      = np.sum(test_pridict==test_lab)
    print('{}个{}测试样本中预测正确个数为：{}'.format(len(test_lab),type,correct))

def test_2(type,X_train,X_test,train_lab,test_lab):
    clf = mySVM(C=1.0)
    clf.fit(X_train,train_lab)

    test_pridict = clf.predict(X_test)
    correct      = np.sum(test_pridict==test_lab)
    print('{}个{}测试样本中预测正确个数为：{}'.format(len(test_lab),type,correct))
    
if __name__ == "__main__":
    test("minist",new_train,new_test,train_lab,test_lab)
    test("手写数字",new_train,write_num_test,train_lab,write_num_label)
    
    #软间隔测试
    test_2("minist",new_train,new_test,train_lab,test_lab)
    test_2("手写数字",new_train,write_num_test,train_lab,write_num_label)