机器学习：支持向量机

1、基本概念

决策面：区分两类的超平面；

分类间隔：在保证决策面方向不变且不会出现错分样本的情况下移动决策面，会在原来的决策面两侧找到两个极限位置（越过该位置就会产生错分现象），俩个极限位置垂直距离就是这个决策面对应的分类间隔；

最优决策面：具有“最大间隔”的决策面；

支持样本点：也称支持向量，最优解对应的两侧虚线穿过的样本点；

2、线性SVM

支持向量机_基本数学模型推导：

3、凸优化问题求解方法

• 无约束优化问题：费马定理（求导数，令导数为零）；
• 有等式约束的优化问题：拉格朗日乘子法；
• 有不等式约束的优化问题：KKT条件；

4、求解基本模型

拉格朗日函数：构造一个函数，使得该函数在可行解区域内与原目标函数完全 一致，而在可行解区域外的数值非常大，那么这个没有约束条件的新目标函数的优化问题就与原来有约束条件的原始目标函数的优化问题等价。

求解步骤：

1. 将有约束的原始目标函数转换为无约束的新构造的拉格朗日目标函数；
2. 使用拉格朗日对偶性，将不易求解的优化问题转化为易求解的优化问题；

第一步：构造拉格朗日函数：

第二步：求新目标函数的最小值：

求内侧最小值：

求外侧最大值：

5、软间隔问题

前面推导的目标函数，假设数据线性可分。实际情况可能有噪声，导致线性不可分；

通过引入松弛变量ξ和惩罚参数C来允许数据点可以处于超平面的错误的一侧；

6、SMO算法

算法原理：根据约束条件随机给α赋值，每次选取两个α，调节两个α 使得目标函数最小，然后再选取别两个α，调节α使得目标函数最小，不断循环，直到达到目标函数最小值；

7、非线性SVM

在线性不可分的情况下，SVM通过某种事先选择的非线性映射（核函数）将输入变量映射到一个高维特征空间，将其变成在高维空间线性可分，在这个高维空间中构造最优分类超平面；

高维映射：方法2的优点是计算量小；

核函数：

8、sklearn中的SVM

sklearn.svm.SVC — scikit-learn 1.1.3 documentation

代码：

SMO算法：

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
data=pd.read_table("SVM_test.txt",sep='	',names=['x1','x2','y'])
dataMat=np.mat(data[['x1','x2']])
labelMat=np.mat(data['y'])

# 选择正样本
plt.scatter(data[data['y']==1]['x1'], data[data['y']==1]['x2'],c='r')   # 正样本散点图
# 选择负样本
plt.scatter(data[data['y']==-1]['x1'], data[data['y']==-1]['x2'],c='b') # 负样本散点图



# SMO实现
# -*- coding:UTF-8 -*-
from time import sleep
import matplotlib.pyplot as plt
import numpy as np
import random
import types

"""
函数说明:随机选择alpha

Parameters:
    i - alpha
    m - alpha参数个数
"""
def selectJrand(i, m):
    j = i                                 #选择一个不等于i的j
    while (j == i):
        j = int(random.uniform(0, m))
    return j

"""
函数说明:修剪alpha
Parameters:
    aj - alpha值
    H - alpha上限
    L - alpha下限
Returns:
    aj - alpah值
"""
def clipAlpha(aj,H,L):
    if aj > H:
        aj = H
    if L > aj:
        aj = L
    return aj

"""
函数说明:SMO算法

Parameters:
    dataMatIn - 数据矩阵
    classLabels - 数据标签
    C - 惩罚参数
    toler - 松弛变量
    maxIter - 最大迭代次数
"""
def smoSimple(dataMatIn, classLabels, C, toler, maxIter):
    #转换为numpy的mat存储
    dataMatrix = np.mat(dataMatIn); labelMat = np.mat(classLabels).transpose()
    #初始化b参数，统计dataMatrix的维度
    b = 0; m,n = np.shape(dataMatrix)
    #初始化alpha参数，设为0
    alphas = np.mat(np.zeros((m,1)))
    #初始化迭代次数
    iter_num = 0
    #最多迭代matIter次
    while (iter_num < maxIter):
        alphaPairsChanged = 0
        for i in range(m):
            #步骤1：计算误差Ei
            fXi = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[i,:].T)) + b
            Ei = fXi - float(labelMat[i])
            #优化alpha，更设定一定的容错率。
            if ((labelMat[i]*Ei < -toler) and (alphas[i] < C)) or ((labelMat[i]*Ei > toler) and (alphas[i] > 0)):
                #随机选择另一个与alpha_i成对优化的alpha_j
                j = selectJrand(i,m)
                #步骤1：计算误差Ej
                fXj = float(np.multiply(alphas,labelMat).T*(dataMatrix*dataMatrix[j,:].T)) + b
                Ej = fXj - float(labelMat[j])
                #保存更新前的aplpha值，使用深拷贝
                alphaIold = alphas[i].copy(); alphaJold = alphas[j].copy();
                #步骤2：计算上下界L和H
                if (labelMat[i] != labelMat[j]):
                    L = max(0, alphas[j] - alphas[i])
                    H = min(C, C + alphas[j] - alphas[i])
                else:
                    L = max(0, alphas[j] + alphas[i] - C)
                    H = min(C, alphas[j] + alphas[i])
                if L==H: print("L==H"); continue
                #步骤3：计算eta
                eta = 2.0 * dataMatrix[i,:]*dataMatrix[j,:].T - dataMatrix[i,:]*dataMatrix[i,:].T - dataMatrix[j,:]*dataMatrix[j,:].T
                if eta >= 0: print("eta>=0"); continue
                #步骤4：更新alpha_j
                alphas[j] -= labelMat[j]*(Ei - Ej)/eta
                #步骤5：修剪alpha_j
                alphas[j] = clipAlpha(alphas[j],H,L)
                if (abs(alphas[j] - alphaJold) < 0.00001): print("alpha_j变化太小"); continue
                #步骤6：更新alpha_i
                alphas[i] += labelMat[j]*labelMat[i]*(alphaJold - alphas[j])
                #步骤7：更新b_1和b_2
                b1 = b - Ei- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[i,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[i,:]*dataMatrix[j,:].T
                b2 = b - Ej- labelMat[i]*(alphas[i]-alphaIold)*dataMatrix[i,:]*dataMatrix[j,:].T - labelMat[j]*(alphas[j]-alphaJold)*dataMatrix[j,:]*dataMatrix[j,:].T
                #步骤8：根据b_1和b_2更新b
                if (0 < alphas[i]) and (C > alphas[i]): b = b1
                elif (0 < alphas[j]) and (C > alphas[j]): b = b2
                else: b = (b1 + b2)/2.0
                #统计优化次数
                alphaPairsChanged += 1
                #打印统计信息
                print("第%d次迭代 样本:%d, alpha优化次数:%d" % (iter_num,i,alphaPairsChanged))
        #更新迭代次数
        if (alphaPairsChanged == 0): 
            iter_num += 1
        else: iter_num = 0
        print("迭代次数: %d" % iter_num)
    return b,alphas

"""
函数说明:分类结果可视化
Parameters:
    dataMat - 数据矩阵
    w - 直线法向量
    b - 直线解决
"""
def showClassifer(data, w, b):
    #绘制样本点
    #选择正样本
    plt.scatter(data[data['y']==1]['x1'], data[data['y']==1]['x2'],c='r')   #正样本散点图
    #选择负样本
    plt.scatter(data[data['y']==-1]['x1'], data[data['y']==-1]['x2'],c='b') #负样本散点图
    #绘制直线
    dataMat=np.mat(data[['x1','x2']])
    dataMat=dataMat.tolist()
    x1 = max(dataMat)[0]
    x2 = min(dataMat)[0]
    a1, a2 = w
    b = float(b)
    a1 = float(a1[0])
    a2 = float(a2[0])
    y1, y2 = (-b- a1*x1)/a2, (-b - a1*x2)/a2
    plt.plot([x1, x2], [y1, y2])
    #找出支持向量点
    for i, alpha in enumerate(alphas):
        if alpha > 0:
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
    plt.show()


"""
函数说明:计算w
Parameters:
    dataMat - 数据矩阵
    labelMat - 数据标签
    alphas - alphas值
"""
def get_w(dataMat, labelMat, alphas):
    alphas, dataMat, labelMat = np.array(alphas), np.array(dataMat), np.array(labelMat)
    w = np.dot((np.tile(labelMat.reshape(1, -1).T, (1, 2)) * dataMat).T, alphas)
    return w.tolist()


if __name__ == '__main__':
    data=pd.read_table("SVM_test.txt",sep='	',names=['x1','x2','y'])
    dataMat=np.mat(data[['x1','x2']])
    labelMat=np.mat(data['y'])
    b,alphas = smoSimple(dataMat, labelMat, 0.6, 0.001, 40)
    w = get_w(dataMat, labelMat, alphas)
    showClassifer(data, w, b)

SVM低维映射到高维：

import matplotlib.pyplot as plt
from sklearn import datasets
import numpy as np
from mpl_toolkits.mplot3d import Axes3D  

x_data, y_data = datasets.make_circles(n_samples=500, factor=.3, noise=0.1)
plt.scatter(x_data[:,0], x_data[:,1], c=y_data)
plt.show()

z_data =  x_data[:,0]**2 + x_data[:,1]**2
ax = plt.figure().add_subplot(111, projection = '3d') 
ax.scatter(x_data[:,0], x_data[:,1], z_data, c = y_data, s = 10) #点为红色三角形   
#显示图像  
plt.show()  


# 3D效果
import mayavi.mlab as mlab
s=y_data
mlab.points3d(x_data[:,0], x_data[:,1], z_data,s,scale_mode='none' )
mlab.axes(xlabel='x', ylabel='y', zlabel='z')
mlab.colorbar()
mlab.show()

SVM线性分类：

import numpy as np
import matplotlib.pyplot as plt
from sklearn import svm

np.random.seed(10)
x_data = np.r_[np.random.randn(20,2)+[-2,-2],np.random.randn(20,2)+[2,2]] # 上下拼接
y_data = [0]*20 + [1]*20

plt.scatter(x_data[:,0],x_data[:,1],c=y_data)
plt.show()

model = svm.SVC(kernel='linear')
model.fit(x_data, y_data)
model.coef_             # 系数
model.intercept_        # 截距项

k = -model.coef_[0][0]/model.coef_[0][1]
d = -model.intercept_/model.coef_[0][1]
x_test = np.array([[-5],[5]])
y_test = x_test*k + d
plt.plot(x_test,y_test,'k')
plt.scatter(x_data[:,0],x_data[:,1],c=y_data)
plt.show()

model.support_vectors_  # 支持向量点
b1 = model.support_vectors_[0]
y_down = k*x_test + (b1[1] - k*b1[0])
b2 = model.support_vectors_[-1]
y_up = k*x_test + (b2[1] - k*b2[0])
k = -model.coef_[0][0]/model.coef_[0][1]
d = -model.intercept_/model.coef_[0][1]
x_test = np.array([[-5],[5]])
y_test = x_test*k + d
plt.plot(x_test,y_test,'k')
plt.plot(x_test,y_down,'b--')
plt.plot(x_test,y_up,'r--')
plt.scatter(x_data[:,0],x_data[:,1],c=y_data)
plt.show()

SVM非线性分类：

import matplotlib.pyplot as plt
import numpy as np
from sklearn.metrics import classification_report
from sklearn import svm
import os
# os.chdir('D:\CDA\File')

data = np.genfromtxt("svm_testset.txt", delimiter=",")  # 读取原始数据集

x0=data[data[:,-1]==0]       # 选择为0类的数据
x1=data[data[:,-1]==1]       # 选择为1类的数据
plt.scatter(x0[:,0], x0[:,1], c='b', marker='o')
plt.scatter(x1[:,0], x1[:,1], c='r', marker='x')
plt.show()

x_data = data[:,:-1]   # 数据
y_data = data[:,-1]    # 标签


model = svm.SVC(kernel='rbf',C=5,gamma=5)
model.fit(x_data, y_data) 
model.score(x_data,y_data)


# 获取数据值所在的范围
x_min, x_max = x_data[:, 0].min() - 0.2, x_data[:, 0].max() + 0.2
y_min, y_max = x_data[:, 1].min() - 0.2, x_data[:, 1].max() + 0.2

# 生成网格矩阵
xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.01),
                     np.arange(y_min, y_max, 0.01))

z = model.predict(np.c_[xx.ravel(), yy.ravel()]) #ravel 降维 左右拼接
z = z.reshape(xx.shape)

# 等高线图
cs = plt.contourf(xx, yy, z) #等高线图
plt.scatter(x0[:,0], x0[:,1], c='b', marker='o')
plt.scatter(x1[:,0], x1[:,1], c='r', marker='x')
plt.show()

9、案例

手写数字的识别：

import os
os.chdir('D:\\CDA\\File')
import numpy as np
from sklearn.svm import SVC
from os import listdir
from sklearn.neighbors import KNeighborsClassifier as kNN

# 代码功能：读取文件，并把文件转换为1*1024    
def img2vector(filename):
    returnVect = np.zeros((1, 1024))
    fr = open(filename)
    for i in range(32):
        lineStr = fr.readline()
        for j in range(32):
            returnVect[0, 32*i+j] = int(lineStr[j])
    return returnVect

# 训练集
hwLabels = []          # 新建空的列表，存放训练集标签
trainingFileList = listdir('trainingDigits')        # listdir()是返回目录下的文件名
m = len(trainingFileList)                           # 样本数量 1934       
trainingMat = np.zeros((m, 1024))      # x矩阵 初始化训练的Mat矩阵,测试集，就是批量处理
for i in range(m):                                  # 从文件名中解析出训练集的类别
    fileNameStr = trainingFileList[i]               # 获得文件的名字；是为了下一步获取每个文件的数字类别
    classNumber = int(fileNameStr.split('_')[0])    # 获得分类的数字；因为分隔后的第一个是代表类别的数字
    hwLabels.append(classNumber)                    # 将获得的类别添加到hwLabels中
    trainingMat[i,:] = img2vector('trainingDigits/%s' % (fileNameStr))     #将每一个文件的1x1024数据存储到trainingMat矩阵中，trainingMat最后是m行1024列的矩阵

#测试集   
testFileList = listdir('testDigits')                # 返回testDigits目录下的文件列表   同样也要对测试集进行相同的处理
testLabels = [] 
mTest = len(testFileList)                           # 测试数据的数量
testMat = np.zeros((mTest, 1024)) 
for i in range(mTest):                              # 从文件中解析出测试集的类别并进行分类测试  
    fileNameStr = testFileList[i]                   # 获得文件的名字
    classNumber = int(fileNameStr.split('_')[0])    # 获得分类的数字
    testLabels.append(classNumber) 
    testMat[i,:] = img2vector('testDigits/%s' % (fileNameStr))     # 获得测试集的1x1024向量,用于训练   

# 参数优化
from sklearn.model_selection import GridSearchCV
params={'C':[0.01,0.1,1,10,100],'gamma':np.linspace(0.001,0.1,10),'class_weight':['','balanced'],'decision_function_shape':['ovo','ovr']}
clf = SVC(kernel='rbf')
grid_search=GridSearchCV(clf,param_grid=params,cv=10,n_jobs=-1,verbose=2)
grid_search.fit(trainingMat,hwLabels)
grid_search.best_params_