由于在MNIST上运行SVM耗时过久,所以这里使用了伪造的数据集,并使用线性核和多项式核进行实验。
#!/usr/bin/env python3and
# -*- coding: utf-8 -*-
import time
import random
import logging
import pandas as pd
from generate_dataset import *
class SVM(object):
def __init__(self, kernel='linear', epsilon=0.001):
self.kernel = kernel
self.epsilon = epsilon # P129,检验在epsilon范围内进行
def init_parameters(self, features, labels):
'''
初始化参数
'''
self.X = features
self.Y = labels
self.b = 0.0
self.n = len(features[0]) # 2,每个样本维度数
self.N = len(features) # 1334,样本数
self.alpha = [0.0] * self.N # list快速初始化
self.E = [self._E(i) for i in range(self.N)]
self.C = 1000
self.Max_Interation = 5000
def K(self, x1, x2):
'''
核函数
'''
if self.kernel == 'linear': # 线性核,等价于没有用核函数
return sum([x1[k] * x2[k] for k in range(self.n)])
elif self.kernel == 'poly': # 多项式核,这里取三次方
return (sum([x1[k] * x2[k] for k in range(self.n)]) + 1) ** 3
def G(self, i):
'''
公式(7.104)
'''
result = self.b
for j in range(N):
result += self.alpha[j] * self.Y[j] * self.K(self.X[i], self.X[j])
return result
def _E(self, i):
'''
公式(7.105)
E(i)为输入xi的预测值与真实输出Yi的误差
'''
return self.G(i) - self.Y[i]
def satisfy_KKT(self, i):
'''
验证是否满足KKT,以便选取alpha
'''
ygx = self.Y[i] * self.G(i)
if abs(self.alpha[i]) < self.epsilon: # alpha与0的误差,(7.111)
return ygx > 1 or ygx == 1
elif abs(self.alpha[i] - self.C) < self.epsilon: # (7.113)
return ygx < 1 or ygx == 1
else:
return abs(ygx - 1) < self.epsilon # (7.112)
def select_two_parameters(self):
'''
按照书7.4.2选择两个变量
'''
index_list = [i for i in range(self.N)]
# 只选择间隔边界上的支持向量点
i1_list_1 = filter(lambda i: self.alpha[i] > 0 and self.alpha[i] < self.C, index_list)
i1_list_1 = list(i1_list_1) # 注意Python3返回的是迭代器对象,Python2返回的是列表
i1_list_2 = list(set(index_list) - set(i1_list_1)) # 其余点
i1_list = i1_list_1
i1_list.extend(i1_list_2) # 先选取0<alpha<C的,再选择其它的
for i in i1_list: # 遍历alpha的下标
if self.satisfy_KKT(i): # 选择一个违反KKT条件的alpha1
continue
# 下面选择alpha2
E1 = self.E[i]
# 第一个值记录E最大差值,第二个值记录alpha2下标
max_parameters = (0, 0)
for j in index_list:
if i == j:
continue
E2 = self.E[j]
if abs(E1 - E2) > max_parameters[0]:
max_parameters = (abs(E1 - E2), j)
return i, max_parameters[1]
def train(self, features, labels):
self.init_parameters(features, labels)
for iteration in range(self.Max_Interation):
# logging.debug('iteration %d:' % iteration)
# 选出alpha1和alpha2的下标
i1, i2 = self.select_two_parameters()
# P126,给出alpha2new的上下界
if self.Y[i1] != self.Y[i2]:
L = max(0, self.alpha[i2] - self.alpha[i1])
H = min(self.C, self.C + self.alpha[i2] - self.alpha[i1])
else:
L = max(0, self.alpha[i2] - self.alpha[i1] - self.C)
H = min(self.C, self.alpha[i2] + self.alpha[i1])
E1 = self.E[i1]
E2 = self.E[i2]
# (7.107)
eta = self.K(self.X[i1], self.X[i1]) + self.K(self.X[i2], self.X[i2]) - self.K(self.X[i1], self.X[i2])
# (7.106),未裁剪的alpha2new
alpha2_new_unc = self.alpha[i2] + self.Y[i2] * (E1 - E2) / eta
# (7.108),进行裁剪
alpha2_new = 0
if alpha2_new_unc > H:
alpha2_new = H
elif alpha2_new_unc < L:
alpha2_new = L
else:
alpha2_new = alpha2_new_unc
# (7.109),求alpha1new
alpha1_new = self.alpha[i1] + self.Y[i1] * self.Y[i2] * (self.alpha[i2] - alpha2_new)
# (7.115)和(7.116),选择b1new和b2new从而选择bnew
b_new = 0
b1_new = -E1 - self.Y[i1] * self.K(self.X[i1], self.X[i1]) * (alpha1_new - self.alpha[i1]) - self.Y[i2] * self.K(self.X[i2], self.X[i1]) * (alpha2_new - self.alpha[i2]) + self.b
b2_new = -E2 - self.Y[i1] * self.K(self.X[i1], self.X[i2]) * (alpha1_new - self.alpha[i1]) - self.Y[i2] * self.K(self.X[i2], self.X[i2]) * (alpha2_new - self.alpha[i2]) + self.b
if alpha1_new > 0 and alpha1_new < self.C:
b_new = b1_new
elif alpha2_new > 0 and alpha2_new < self.C:
b_new = b2_new
else:
b_new = (b1_new + b2_new) / 2
self.alpha[i1] = alpha1_new
self.alpha[i2] = alpha2_new
self.b = b_new
# (7.117),还需要更新E,这里需要bnew的值
self.E[i1] = self._E(i1)
self.E[i2] = self._E(i2)
def predict_one(self, feature):
result = self.b
for i in range(self.N):
result += self.alpha[i] * self.Y[i] * self.K(feature, self.X[i])
if result > 0:
return 1
else:
return -1
def predict(self, features):
results = []
for feature in features:
results.append(self.predict_one(feature))
return results
if __name__ == '__main__':
logger = logging.getLogger()
logger.setLevel(logging.DEBUG)
print('Start reading data:')
time1 = time.time()
# 选取2/3为训练数据,1/3为测试数据
train_features, train_labels, test_features, test_labels = generate_dataset(2000, visualization=False)
# print(len(train_features)) # list 1334
# print(len(train_features[0])) # list 2
# print(len(train_labels)) # list 1334
# print(len(test_features)) # list 666
# print(len(test_labels)) # list 666
time2 = time.time()
print('read data cost %f seconds' % (time2 - time1))
print('Start training:')
svm = SVM('poly')
svm.train(train_features, train_labels)
time3 = time.time()
print('training cost %f seconds' % (time3 - time2))
print('Start predicting:')
test_predict = svm.predict(test_features)
time4 = time.time()
print('predicting cost %f seconds' % (time4 - time3))
accuracy = sum(test_labels[i] == test_predict[i] for i in range(len(test_labels))) / len(test_labels)
print('The accuracy is %f!' % accuracy)
'''
output:(linear)
Start reading data:
read data cost 0.013078 seconds
Start training:
training cost 4.556223 seconds
Start predicting:
predicting cost 3.825693 seconds
The accuracy is 0.887387!
'''
'''
output:(poly)
Start reading data:
read data cost 0.013761 seconds
Start training:
training cost 4.902359 seconds
Start predicting:
predicting cost 4.038038 seconds
The accuracy is 0.875375!
'''