# 支持向量机

## 算法原理：

### 一、寻找最大间隔

$L\left(W,\alpha ,,b\right)=\frac{1}{2}{W}^{T}W+\sum _{i=1}^{n}{\alpha }_{i}\left(1-{y}^{i}\left({W}^{T}{X}^{i}+b\right)\right)$

$\frac{\mathrm{\partial }L}{\mathrm{\partial }W}=W-\sum _{i=1}^{n}{\alpha }_{i}{y}^{i}{X}^{i}=0$

$\frac{\mathrm{\partial }l}{\mathrm{\partial }b}=\sum _{i=1}^{n}{\alpha }_{i}{y}^{i}=0$

$W=\sum _{i=1}^{n}{\alpha }_{i}{y}^{i}{X}^{i}$

$L\left(W,\alpha ,b\right)=\frac{1}{2}{W}^{T}W+\sum _{i=1}^{n}{\alpha }_{i}-\sum _{i=1}^{n}{\alpha }_{i}{y}^{i}{W}^{T}{X}_{i}-b\sum _{i=1}^{n}{\alpha }_{i}{y}^{i}$

$\frac{1}{2}\left(\sum _{i=1}^{n}{\alpha }_{i}{y}^{i}{X}^{i}{\right)}^{T}\left(\sum _{i=1}^{n}{\alpha }_{i}{y}^{i}{X}^{i}\right)+\sum _{i=1}^{n}{\alpha }_{i}-\sum _{i=1}^{n}{\alpha }_{i}{y}^{i}\left(\sum _{i=1}^{n}{\alpha }_{i}{y}^{i}{X}^{i}{\right)}^{T}{X}_{i}$

$J\left(\alpha \right)=\sum _{i=1}^{n}{\alpha }_{i}-\frac{1}{2}\sum _{i=1}^{n}\sum _{j=1}^{n}{\alpha }_{i}{\alpha }_{j}{y}^{i}{y}^{j}\left({X}^{i}{\right)}^{T}{X}^{i}$

$\sum _{i=1}^{n}{\alpha }_{i}{y}^{i}=0$

${\alpha }_{i}\ge 0$

### 二、引入松弛变量

${y}^{i}\left({W}^{T}{X}^{i}+b\right)\ge 1-{\xi }_{i}$

${\xi }_{i}\ge 0$

$L\left(W,b,{\xi }_{i},\alpha ,\mu \right)=\frac{1}{2}{W}^{T}W+C\sum _{i=1}^{n}{\xi }_{i}+\sum _{i=1}^{n}{\alpha }_{i}\left(1-{\xi }_{i}-{y}^{i}\left({W}^{T}{X}^{i}+b\right)\right)-\sum _{i=1}^{n}{\mu }_{i}{\xi }_{i}$

(这里将${\xi }_{i}\ge 0$$\xi_i\geq0$变成$-{\xi }_{i}\le 0$$-\xi_i\leq0$的形式，以适用拉格朗日乘子法)

$W=\sum _{i=1}^{n}{\alpha }_{i}{y}^{i}{X}^{i}$

$\sum _{i=1}^{n}{\alpha }_{i}{y}^{i}=0$

$J\left(\alpha \right)=\sum _{i=1}^{n}{\alpha }_{i}-\frac{1}{2}\sum _{i=1}^{n}\sum _{j=1}^{n}{\alpha }_{i}{\alpha }_{j}{y}^{i}{y}^{j}\left({X}^{i}{\right)}^{T}{X}^{i}$

$C-{\alpha }_{i}-{\mu }_{i}=0$

$\sum _{i=1}^{n}{\alpha }_{i}=0$

${\alpha }_{i}\ge 0$

${\mu }_{i}\ge 0$

$C-{\alpha }_{i}-{\mu }_{i}=0$

$0\le {\alpha }_{i}\le C$

## 三、SMO求解算法(简化版)

1）该段代码进行需要的各种库的导入，最后一行%matplotlib inline是为了保证在jupyter notebook 的运行环境下能够在浏览器上输出绘制的图像。

import numpy as np
import random
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline    

file_name = r"D:\python_data\MachineLearningInaction\machinelearninginaction\Ch06\testSet.txt"
file.head()

3）可以看到所读取的txt文件的前几行如下图，每个样本有2个特征，类标签分别为+1或-1。

4）接下来再绘制一张样本点分布的散点图，蓝色的是正例，即+1类的点；红色的是负例，即-1类的点。横坐标为“factor1”，纵坐标为“factor2”。

positive = file[file["class"] == 1]
negative = file[file["class"] == -1]
fig, ax = plt.subplots(figsize=(10,5))
ax.scatter(positive["factor1"], positive["factor2"], s=30, c="b", marker="o", label="class 1")
ax.scatter(negative["factor1"], negative["factor2"], s=30, c="r", marker="x", label="class -1")
ax.legend()
ax.set_xlabel("factor1")
ax.set_ylabel("factor2")

5）运行后，在浏览器上输出了如下的图像。

6）现在对数据已经有了清晰直观的了解，接下来需要把pandas读取出来的整个表格进行切分，前2列作为样本点的特征矩阵，最后一列作为类标签向量。下面函数输入的是pandas都取出来的文件file，调用as_matrix方法将其转化为矩阵orig_data，将这个矩阵去除最后一列的部分作为样本点的特征矩阵data_mat，而最后一列作为类标签向量label_mat。

def load_data_set(file):
orig_data = file.as_matrix()
cols = orig_data.shape[1]
data_mat = orig_data[:,0:cols-1]
label_mat = orig_data[:,cols-1:cols]
return data_mat, label_mat

7）接下来调用该函数，得到样本点的特征矩阵data_mat和对应的类标签向量label_mat。

data_mat, label_mat = load_data_set(file) 

8）还需要定义一个select_jrand函数，该函数的功能是给定输入i和m，随机输出一个0到m之间的与i不同的整数。这个函数在后面将会用到。

def select_jrand(i, m):
j = i
while(j==i):
j = int(random.uniform(0,m))
return j

9）此外还需要定义一个clip_alpha函数，该函数的功能是输入将aj限制在L和H之间

def clip_alpha(aj, H, L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj

10）接下来是主体功能的实现函数，由于这个函数太长，为了方便理解，将以下代码分成一小段一小段的，这样就可以理解每一小段的功能。

10.1）循环外的初始化工作：首先该函数的输入为特征矩阵、类标签向量，常数C，容错率，最大迭代次数。首先将输入的特征矩阵和类标签向量都转化为矩阵(若特征矩阵和类标签向量已经是矩阵，则不需要这个步骤，若输入的是列表的形式，则需要该步骤)，设置b为0，m和n分别为特征矩阵的行数和列数。创建一个m行1列的$\alpha$$\alpha$向量(即有多少个样本点就有多少个$\alpha$$\alpha$)，迭代次数初始化为0(每次迭代在$\alpha$$\alpha$向量中寻找一对可以改变的$\alpha$$\alpha$，若没有找到则迭代次数+1，达到最大迭代次数还没有找到可以改变的$\alpha$$\alpha$，说明已经$\alpha$$\alpha$已经优化的差不多了，则退出循环)。

10.2）内循环的初始化工作：现在开始第一轮迭代，将改变的一对$\alpha$$\alpha$的数量alpha_pairs_changed初始化为0，然后遍历m个$\alpha$$\alpha$依次作为第一个$\alpha$$\alpha$，再在(0,m)范围内随机的找到第二个$\alpha$$\alpha$(第二个$\alpha$$\alpha$与第一个不能是同一个，这就用到了之前定义的select_jrand函数)。

10.3)第一小段代码：首先从i等于0开始，已知公式$W=\sum _{i}^{n}{\alpha }_{i}{y}^{i}{X}^{i}$$W=\sum_i^n\alpha_iy^iX^i$，若第0号数据点是支持向量，则该数据点必然在虚线${W}^{T}X+b+1=0$$W^TX+b+1=0$${W}^{T}X+b-1=0$$W^TX+b-1=0$上，即${W}^{T}X+b=-1$$W^TX+b=-1$${W}^{T}X+b=+1$$W^TX+b=+1$

10.4）第二小段代码：由于有一个限制条件即$\sum _{i}^{n}{\alpha }_{i}{y}^{i}=0$$\sum_i^n\alpha_iy^i=0$，且任意的$\alpha$$\alpha$在(0,C)之间，所以当调整其中一个${\alpha }_{i}$$\alpha_i$的时候，另一个${\alpha }_{j}$$\alpha_j$的改变必须有一定限制，即${y}^{i}{\alpha }_{i}+{y}^{j}{\alpha }_{j}=k$$y^i\alpha_i+y^j\alpha_j=k$。即如下图：${\alpha }_{i}$$\alpha_i$${\alpha }_{j}$$\alpha_j$必须在两条线段上。

10.5）第三小段代码：下面为${\alpha }_{j}$$\alpha_j$的更新公式：

${\alpha }_{j}={\alpha }_{j}-\frac{{y}^{j}\left({E}_{i}-{E}_{j}\right)}{\eta }$

10.6）第四小段代码：对于样本点i，已经更新了其${\alpha }_{i}$$\alpha_i$值，若b值也是正确的，则该点应该从更新之前的不在${W}^{T}X+b±1=0$$W^TX+b\pm1=0$上，更新到位于${W}^{T}X+b±1=0$$W^TX+b\pm1=0$上，亦即将样本点i带入，得到${E}_{i}^{new}={{W}_{new}}^{T}{X}^{i}+{b}^{new}-{y}^{i}=0$$E_i^{new}={W_{new}}^TX^i+b^{new}-{y^i}=0$

${b}^{new}={b}^{old}-{E}_{i}^{old}+{{W}_{old}}^{T}{X}^{i}-{{W}_{new}}^{T}{X}^{i}$

${y}^{i}\left(\alpha {{}_{i}}^{old}-\alpha {{}_{i}}^{new}\right)<{X}^{i},{X}^{i}>+{y}^{j}\left({\alpha }_{j}^{old}-{\alpha }_{j}^{new}\right)<{X}^{j},{X}^{i}>$

10.7)收尾工作：若程序运行到此处，则已经更新了一对$\alpha$$\alpha$，并且也更新了$b$$b$值，则将alpha_pairs_changed加上1，代表更新了一对，继续更新下一对。若程序没有执行更新$\alpha$$\alpha$的操作，即alpha_pairs_changed=0，则将iter+1，代表迭代了一次，否则将iter置为1，重新开始迭代。若程序迭代的次数达到了最大迭代次数，还没有$\alpha$$\alpha$被更新，说明已经更新的很好了，则可以退出while循环，返回b值和$\alpha$$\alpha$向量。

def smo_simple(data_mat, class_label, C, toler, max_iter):
# 循环外的初始化工作
data_mat = np.mat(data_mat); label_mat = np.mat(class_label)
b = 0
m,n = np.shape(data_mat)
alphas = np.zeros((m,1))
iter = 0
while iter < max_iter:
# 内循环的初始化工作
alpha_pairs_changed = 0
for  i in range(m):
# 第一小段代码
WT_i = np.dot(np.multiply(alphas, label_mat).T, data_mat)
f_xi = float(np.dot(WT_i, data_mat[i,:].T)) + b
Ei = f_xi - float(label_mat[i])
if((label_mat[i]*Ei < -toler) and  (alphas[i] < C)) or \
((label_mat[i]*Ei > toler) and (alphas[i] > 0)):
j = select_jrand(i, m)
WT_j = np.dot(np.multiply(alphas, label_mat).T, data_mat)
f_xj = float(np.dot(WT_j, data_mat[j,:].T)) + b
Ej = f_xj - float(label_mat[j])
alpha_iold = alphas[i].copy()
alpha_jold = alphas[j].copy()

# 第二小段代码
if (label_mat[i] != label_mat[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 H == L :continue

# 第三小段代码
eta = 2.0 * data_mat[i,:]*data_mat[j,:].T - data_mat[i,:]*data_mat[i,:].T - \
data_mat[j,:]*data_mat[j,:].T
if eta >= 0: continue
alphas[j] = (alphas[j] - label_mat[j]*(Ei - Ej))/eta
alphas[j] = clip_alpha(alphas[j], H, L)
if (abs(alphas[j] - alpha_jold) < 0.00001):
continue

alphas[i] = alphas[i] + label_mat[j]*label_mat[i]*(alpha_jold - alphas[j])

# 第四小段代码
b1 = b - Ei + label_mat[i]*(alpha_iold - alphas[i])*np.dot(data_mat[i,:], data_mat[i,:].T) +\
label_mat[j]*(alpha_jold - alphas[j])*np.dot(data_mat[i,:], data_mat[j,:].T)
b2 = b - Ej + label_mat[i]*(alpha_iold - alphas[i])*np.dot(data_mat[i,:], data_mat[j,:].T) +\
label_mat[j]*(alpha_jold - alphas[j])*np.dot(data_mat[j,:], data_mat[j,:].T)
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
alpha_pairs_changed += 1
if (alpha_pairs_changed == 0): iter += 1
else: iter = 0
return b, alphas

11）接下来调用该函数，得到b值和$\alpha$$\alpha$向量，由于$\alpha$$\alpha$的值大多为0，所以只打印大于零的$\alpha$$\alpha$值。

b, alphas = smo_simple(data_mat, label_mat, 0.6, 0.001, 10)
print(b, alphas[alphas>0])

12）辛辛苦苦走到这里，就看到几个所谓的$\alpha$$\alpha$和一个b值，真的不够。“锦瑟无端五十弦，一弦一柱思华年”，每一个算法的推导，每一段代码的理解，都给我留下了深刻的印象。过程的享受是一种润物无声的细腻，但我想要的更是成功之后那会“当凌绝顶，一览众山小”的快意，更要那“冲天香阵透长安，满城尽戴黄金甲”的震撼。于是，我又敲下了如下代码，得到了，那虽不太完美，却凝聚了我心血的分、割、线

support_x = []
support_y = []
class1_x = []
class1_y = []
class01_x = []
class01_y = []
for i in range(100):
if alphas[i] > 0.0:
support_x.append(data_mat[i,0])
support_y.append(data_mat[i,1])
for i in range(100):
if label_mat[i] == 1:
class1_x.append(data_mat[i,0])
class1_y.append(data_mat[i,1])
else:
class01_x.append(data_mat[i,0])
class01_y.append(data_mat[i,1])
w_best = np.dot(np.multiply(alphas, label_mat).T, data_mat)
fig, ax = plt.subplots(figsize=(10,5))
ax.scatter(support_x, support_y, s=100, c="y", marker="v", label="support_v")
ax.scatter(class1_x, class1_y, s=30, c="b", marker="o", label="class 1")
ax.scatter(class01_x, class01_y, s=30, c="r", marker="x", label="class -1")
lin_x = np.linspace(3,6)
lin_y = (-float(b) - w_best[0,0]*lin_x) / w_best[0,1]
plt.plot(lin_x, lin_y, color="black")
ax.legend()
ax.set_xlabel("factor1")
ax.set_ylabel("factor2")

+++++++++++刚好在这一刻，2018年的蓝色月全食发生了+++++++++++++++
++++++向月全食致敬，祝我一年好运气，2018.1.31，19：41月全食++++++++

09-20 2606

10-15 1万+

01-09

09-08 4046

05-29 1179

05-06

#### Python可以这样学（第九季 机器学习案例与实战）

01-24
©️2020 CSDN 皮肤主题: 大白 设计师: CSDN官方博客

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、C币套餐、付费专栏及课程。