机器学习之支持向量机（SVM）的求解方法

前言

支持向量机就是寻找一个超平面，将不同的样本分分隔开来，其中间隔分为硬间隔和软间隔，硬间隔就是不允许样本分错，而软间隔就是允许一定程度上样本存在偏差，后者更符合实际。

支持向量机思路简单但是求解过程还是比较复杂，需要将原函数通过拉格朗日乘子法并附上KKT条件是的问题有强对偶性，再使用SMO等算法进行高效的求解。
推导过程可以参考：
机器学习之支持向量机之线性可分型原理介绍及代码实现（SVM）
下面主要实现模型的求解方法。

梯度下降法

梯度下降法是一种比较普适的方法，当模型无法得出解析解，或者解析解求解困难的时候，都可以使用梯度下降法来近似求解。因为梯度下降法需要一轮轮迭代，也需要定义损失函数，因此一般而言，梯度下降法只能获得近似最优解。
SVM可以使用梯度下降法求解，不过得出的解大概率只是近似解，并且不一定满足SVM的公式里的约束条件。

# coding=utf8
import sys;
import random;
import numpy as np
import math

EPS = 0.000000001  # 很小的数字，用于判断浮点数是否等于0
def load_data(filename, dim):
    '''
    输入数据格式: label\tindex1:value1\tindex2:value2\tindex3:value3..., 其中index是特征的编号, 从1开始
    data的数据格式: [[label, sample],[label, sample], ...],  其中sample: [v0, v1, v2, v3, ..., v[dim]]
    '''
    label_ = []
    data_ = []
    for line in open(filename, 'rt'):
        sample = [0.0 for v in range(0, dim + 1)];
        line = line.rstrip("\r\n\t ");
        fields = line.split("\t");
        label = int(fields[0]);  # LABEL取值: 1 or -1
        sample[0] = 1.0;  # sample第一个元素用于存x0特征, 默认置为1.0[方便把 WX+b => WX]
        for field in fields[1:]:
            kv = field.split(":");
            idx = int(kv[0]);  # ensure idx >= 1
            val = float(kv[1]);
            sample[idx] = val;
        label_.append(label)
        data_.append(sample)
    
    label_ = np.array(label_)
    data_ = np.array(data_)
    return label_, data_
    


def svm_train(train_x, train_y, dim, iterations, lm, lr):
    '''
    data4train:数据集
    dim:样本特征维度
    W:SVM模型的权重
    iterations:迭代次数
    目标函数: obj(, W) = (对所有SUM{max{0, 1 - W*X*y}}) + lm / 2 * ||W||^2, 即：hinge+L2
    '''
    X = np.zeros(dim + 1)  #  => 
    grad = np.zeros(dim + 1) # 梯度
    num_train = len(train_x);
    global W
    for i in range(0, iterations):
        # 每次迭代随机选择一个训练样本
        index = random.randint(0, num_train - 1);
        y = train_y[index]  # y其实就是label
        X = train_x[index]

        # 计算梯度
        # for j in range(0, dim + 1):
        #	grad = lm * W[j];
        WX = 0.0
        WX += (W * X).sum() 

        if 1 - WX * y > 0:
            grad = lm * W - X * y
        else:  # 1-WX *y <= 0的时候，目标函数的前半部分恒等于0, 梯度也是0
            grad = lm * W - 0;
        # 更新权重, lr是学习速率
        W = W - lr * grad
        


def svm_predict(x, y, dim, W):
    num_test = len(x);
    num_correct = 0;
    for i in range(0, num_test):
        target = y[i]  # 即label
        X = x[i]  # 即sample
        sum = 0.0;
        sum += (X * W).sum()
        predict = -1;
        # print sum;
        if sum > 0:  # 权值>0，认为目标值为1
            predict = 1;
        if predict * target > 0:  # 预测值和目标值符号相同
            num_correct += 1;

    return num_correct * 1.0 / num_test;


if __name__ == "__main__":
    # 设置参数
    epochs = 10;  # 迭代轮数
    iterations = 10;  # 每一轮中梯度下降迭代次数, 这个其实可以和epochs合并为一个参数
    lm = 0.0001;  # lambda, 对权值做正则化限制的权重
    lr = 0.01;  # lr, 是学习速率，用于调整训练收敛的速度
    dim = 1000;  # dim, 特征的最大维度, 所有样本不同特征总数
    W = np.zeros(dim + 1)  # 权值
    # 导入测试集&训练集
    train_y, train_x = load_data("train.txt", dim)
    test_y, test_x = load_data("test.txt", dim)

    # # 训练, 实际迭代次数=epochs * iterations
    for i in range(0, epochs):
        svm_train(train_x, train_y, dim, iterations, lm, lr);
        accuracy = svm_predict(test_x, test_y, dim, W);

        print("epoch:%d\taccuracy:%f" % (i, accuracy));
    # 输出结果权值
    for i in range(0, dim + 1):
        if math.fabs(W[i]) > EPS:
            print("权值W%d\t%f" % (i, W[i]));
    print(W)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109

这个是参考支持向量机SVM-手写笔记&手动实现这篇博客的，代码改成numpy进行运算了，进行改动的过程也是读懂代码的过程。

这个例子使用的梯度下降法，损失函数应该大概也许相当于用了max(0, x)。整体思路就是，每次就选一个样本点进行参数更新，如果这个样本点对于当前的参数能够正确分类，那么就不更新，如果不能正确分类，就更新。

运行结果如下：

因为每次随机选择样本，那么其实找到的那个超平面大概率不会将样本集完全分开的，但是事实证明梯度下降法还是有效果的，所有，应该可以勉强认为训练出来的是软间隔的SVM吧？

SMO算法

贴上b站视频，以便日后再去看，实话说目前也没有完全走通这个算法的流程。
快速理解SMO算法
SMO算法思路很简单，因为存在约束条件

所以每次更新两个α，剩余看成常量，α能够通过看成的常量导出，然后满足约束条件的情况下求出极值，更新α。每次都更新两个，当固定其他的α时，能够求出选取α的更新的解析解，所以就算起来非常快。
那么该如何选择呢？
第一个αi应该选择违反KKT条件最大的。
第二个αj应该选择于第一个αi差值最大的。
这样能够保证每次更新都是向最快的方向进行更新。
思路是很简单，但是实践起来还是很困难的，因为里面涉及到许多的约束条件，不同情况下需要分类讨论等等。
直接贴上别人的代码吧

# -*- coding: utf-8 -*-
from numpy import *
import matplotlib.pyplot as plt
import random
def loadDataSet(filename): 
    dataMat=[]
    labelMat=[]
    fr=open(filename)
    for line in fr.readlines():
        lineArr=line.strip().split('\t')
        dataMat.append([float(lineArr[0]),float(lineArr[1])])
        labelMat.append(float(lineArr[2]))
    return dataMat,labelMat 
class optStruct:
    def __init__(self,dataMatIn, classLabels, C, toler):  
        self.X = dataMatIn  # 数据的特征
        self.labelMat = classLabels # 数据的标签
        self.C = C # 用于调整间隔大小
        self.tol = toler # 容错率
        self.m = shape(dataMatIn)[0] # 样本个数 
        self.alphas = mat(zeros((self.m,1))) # 需要求解的α
        self.b = 0 # 偏置
        self.eCache = mat(zeros((self.m,2))) # 第一列表示是否是个有效标志位，第二列存误差值E
    def print_m(self):
        print("self.X", self.X)
        print("self.labelMat", self.labelMat)
        print("self.C", self.C)
        print("self.tol", self.tol)
        print("self.m", self.m)
        print("self.alphas", self.alphas)
        print("self.b", self.b)
        print("self.eCache", self.eCache)

# 随机选取一个J值
def selectJrand(i,m): 
    j=i
    while (j==i):
        j=int(random.uniform(0,m))
    return j
 
# 根据关于α_1与α_2的优化问题对应的约束问题分析，对α进行截取约束
# 保证取值范围符合约束条件 
def clipAlpha(aj,H,L):  
    if aj>H:
        aj=H
    if L>aj:
        aj=L
    return aj
        
#计算每个样本点k的Ek值，就是计算误差值=预测值-标签值        
def calcEk(oS, k): 
    fXk = float(multiply(oS.alphas,oS.labelMat).T*(oS.X*oS.X[k,:].T) + oS.b)
    Ek = fXk - float(oS.labelMat[k])
    return Ek
 

#内循环的启发式方法，获取最大差值|Ei-Ej|对应的Ej的索引J
def selectJ(i, oS, Ei):
    # 传入第一个α对应的索引i和误差值Ei

    maxK = -1 #用于保存临时最大索引
    maxDeltaE = 0 #用于保存临时最大差值--->|Ei-Ej|
    Ej = 0 #保存我们需要的Ej误差值
    oS.eCache[i] = [1,Ei]
    # 获取有效误差的下标
    validEcacheList = nonzero(oS.eCache[:,0].A)[0]  
    if (len(validEcacheList)) > 1: # 如果存在符合条件的αj
        for k in validEcacheList:
            if k == i:
                continue
            Ek = calcEk(oS, k)
            deltaE = abs(Ei - Ek)
            if (deltaE > maxDeltaE): 
                maxK = k
                maxDeltaE = deltaE
                Ej = Ek
        return maxK, Ej # 找到与αi差值最大的αj(最优更新方向)
    else: # 不存在就任意选取
        j = selectJrand(i, oS.m)
        Ej = calcEk(oS, j)
    return j, Ej

# 更新Ek
def updateEk(oS, k): 
    Ek = calcEk(oS, k)
    oS.eCache[k] = [1,Ek]
 
# 实现内循环函数
def innerL(i, oS):
    # 确定了第一个αi
    # 计算出Ei
    Ei = calcEk(oS, i) 
    # 如果违背了KKT条件
    if ((oS.labelMat[i]*Ei < -oS.tol) and (oS.alphas[i] < oS.C)) or ((oS.labelMat[i]*Ei > oS.tol) and (oS.alphas[i] > 0)): 
        # 选择i对应确定的αj
        j,Ej = selectJ(i, oS, Ei)
        # 因为后面还需要用到原始值，因此需要记录下来
        alphaIold = oS.alphas[i].copy()
        alphaJold = oS.alphas[j].copy()

        # 确定αj的取值范围
        if (oS.labelMat[i] != oS.labelMat[j]): # 标签异号时
            L = max(0, oS.alphas[j] - oS.alphas[i])
            H = min(oS.C, oS.C + oS.alphas[j] - oS.alphas[i])
        else: # 标签同号
            L = max(0, oS.alphas[j] + oS.alphas[i] - oS.C)
            H = min(oS.C, oS.alphas[j] + oS.alphas[i])
        if L==H:
            print("L==H")
            return 0
        # 2*xi*xj - xi^2 - xj^2 > 0 <=>
        # xi^2 + xj^2 - 2xi*xj < 0  <=>
        # （xi - xj）^2 < 0 <=>
        # 没什么好算的  
        eta = 2.0 * oS.X[i,:]*oS.X[j,:].T-oS.X[i,:]*oS.X[i,:].T-oS.X[j,:]*oS.X[j,:].T
        if eta >= 0:
            print("eta>=0")
            return 0
        # 满足条件后更新αj
        oS.alphas[j] -= oS.labelMat[j]*(Ei - Ej)/eta 
        oS.alphas[j] = clipAlpha(oS.alphas[j],H,L)
        # 更新Ek 
        updateEk(oS, j)
        # 检测变化量是否显著
        if (abs(oS.alphas[j] - alphaJold) < oS.tol): 
            print("j not moving enough")
            return 0
        # αj 的变化量显著则更新 αi
        oS.alphas[i] += oS.labelMat[j]*oS.labelMat[i]*(alphaJold - oS.alphas[j])
        # 更新Ek
        updateEk(oS, i)
        # 更新b 
        b1 = oS.b - Ei- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[i,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[i,:]*oS.X[j,:].T
        b2 = oS.b - Ej- oS.labelMat[i]*(oS.alphas[i]-alphaIold)*oS.X[i,:]*oS.X[j,:].T - oS.labelMat[j]*(oS.alphas[j]-alphaJold)*oS.X[j,:]*oS.X[j,:].T
        #根据统计学习方法中阈值b在每一步中都会进行更新，
        #1.当新值alpha_1不在界上时(0
        #2.当新值alpha_2不在界上时(0 < alpha_2 < C)，b_new的计算规则为：b_new = b2
        #3.否则当alpha_1和alpha_2都不在界上时，b_new = 1/2(b1+b2)
        if (0 < oS.alphas[i]<oS.C):
            oS.b = b1
        elif (0 < oS.alphas[j]<oS.C):
            oS.b = b2
        else:
            oS.b = (b1 + b2)/2.0
        # 更新完毕
        return 1
    else: # 没有违反KKT条件
        return 0
 
 
# 根据西瓜书6.37计算W参数
def calcWs(dataMat, labelMat, alphas):
    alphas, dataMat, labelMat = array(alphas), array(dataMat), array(labelMat)
    w = dot((tile(labelMat.reshape(1, -1).T, (1, 2)) * dataMat).T, alphas)
    return w.tolist()
 
 
def smoP(dataMatIn, classLabels, C, toler, maxIter): 
    # 初始化数据
    oS = optStruct(mat(dataMatIn),mat(classLabels).transpose(),C,toler)
    # oS.print_m()
    iter = 0 # 当前迭代次数
    entireSet = True # 标志是否应该遍历整个数据集
    alphaPairsChanged = 0  # 标志一次循环中α更新的次数

    # 当迭代次数没有达到要求并且（有α被更新或者应该遍历整个数据集）则继续迭代
    while (iter < maxIter) and ((alphaPairsChanged > 0) or (entireSet)):
        alphaPairsChanged = 0
        if entireSet:
            for i in range(oS.m):
                alphaPairsChanged += innerL(i,oS) # innerL成功更新return 1
                print("fullSet, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged)) 
            iter += 1
        else:
            # 获取非边界值得索引
            nonBoundIs = nonzero((oS.alphas.A > 0) * (oS.alphas.A < C))[0]
            for i in nonBoundIs: 
                alphaPairsChanged += innerL(i,oS)
                print("non-bound, iter: %d i:%d, pairs changed %d" % (iter,i,alphaPairsChanged))
            iter += 1

        # 实现交替遍历
        if entireSet:
            entireSet = False
        elif (alphaPairsChanged == 0): 
            entireSet = True
        print("iteration number: %d" % iter)
    # oS.print_m()
    return oS.b,oS.alphas,oS
 
 
def showClassifer(dataMat,labelMat,alphas, w, b):
    data_plus = []                                  
    data_minus = []
    for i in range(len(dataMat)):
        if labelMat[i] > 0:
            data_plus.append(dataMat[i])
        else:
            data_minus.append(dataMat[i])
    data_plus_np = array(data_plus)              
    data_minus_np = array(data_minus)            
    plt.scatter(transpose(data_plus_np)[0], transpose(data_plus_np)[1], s=30, alpha=0.7)   
    plt.scatter(transpose(data_minus_np)[0], transpose(data_minus_np)[1], s=30, alpha=0.7) 
    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 0.6>abs(alpha) > 0:
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='red')
        if 50==abs(alpha) :
            x, y = dataMat[i]
            plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolor='yellow')
    plt.show()

if __name__ == "__main__":
    x = [[1,8],[3,20],[1,15],[3,35],[5,35],[4,40],[7,80],[6,49],[1.5,25],[3.5,45],[4.5,50],[6.5,15],[5.5,20],[5.8,74],[2.5,5]]
    y = [1,1,-1,-1,1,-1,-1,1,-1,-1,-1,1,1,-1,1]
    b, alphas, oS = smoP(dataMatIn=x,classLabels=y,C=50, toler=0.001,maxIter=400)
    w = calcWs(x,y,alphas)
    showClassifer(x,y,alphas, w, b)
    # oS.print_m() # 打印参数
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227

运行结果：


其中画圈的就是支持向量机的支持向量，也就是α不为0的样本是对超平面位置有影响的样本点。

SMO算法实现细节可能还是没完全明白，以后懂了再来补吧。

参考

支持向量机SVM-手写笔记&手动实现
https://www.cnblogs.com/ssyfj/p/13363526.html
SVM SMO算法代码详细剖析