python高斯核函数_支持向量机（SVM）模型python复现 - SMO算法；核技巧；高斯核函数...

'''支持向量机-------------------------------------数据集：Mnist训练集数量：10000测试集数量：1000---------------------------------训练结果：初始化时间：327s训练时间：285s正确率：0.99'''

import numpy as np

import time

from itertools import islice

import math

'''readlines方法读取训练数据--------输入：file：数据地址--------输出：data：转换后的数据label：转换后的标记'''

def LoadData(file):

# 打开数据文件

fr = open(file, 'r')

# 准备两个list存储data和label

data = []

label = []

# 逐行读取数据，使用islice可以直接跳过第一行表头进行逐行读取

for line in islice(fr,1,None):

# 对每一行消去空格，并根据','进行分割

splited = line.strip().split(',')

# 分割后的第一个元素是Label，跳过label遍历所有特征值

# 归一化

int_line = [int(num) / 255 for num in splited[1:]]

# 逐行存储数据

data.append(int_line)

# 将问题转换为二分类问题

if int(splited[0]) == 0:

label.append(1)

else:

label.append(-1)

# 转换成ndarray形式方便后续计算

data = np.array(data)

label = np.array(label)

# 返回数据的特征部分和标记部分

return data, label

'''定义SVM类'''

class SVM:

'''初始化参数'''

def __init__(self, train_data, train_label, sigma, C, toler, itertime):

self.train_data = train_data # 训练集数据

self.train_label = train_label # 训练集标记

self.m, self.n = np.shape(train_data) # self.m为训练集样本容量，self.n为特征数量

self.sigma = sigma # 高斯核分母上的超参数

self.KernalMatrix = self.CalKernalMatrix() # 高斯核矩阵

self.alpha = np.zeros(self.m) # 初始化拉格朗日向量，长度为训练集样本容量

self.b = 0 # 初始化参数b

self.C = C # 惩罚参数

self.toler = toler # 松弛变量

self.itertime = itertime # 迭代次数

self.E = [float(-1 * y) for y in self.train_label] # 初始化Elist，因为alpha和b初始值为0，因此E的初始值为训练集标记的值

'''计算高斯核矩阵'''

def CalKernalMatrix(self):

# 初始化高斯核矩阵，矩阵为m*m大小

KernalMatrix = [[0 for i in range(self.m)] for j in range(self.m)]

# 遍历每一个样本

for i in range(self.m):

# 首先得到一个样本的数据

X = self.train_data[i]

# 仅遍历从i到self.m的样本，因为Kij和Kji的值相同，所以只需要计算一次

for j in range(i, self.m):

# 得到另一个样本的数据

Z = self.train_data[j]

# 计算核函数

K = np.exp(-1 * np.dot(X - Z, X - Z) / (2 * np.square(self.sigma)))

# 存储在核矩阵中对称的位置

KernalMatrix[i][j] = K

KernalMatrix[j][i] = K

KernalMatrix = np.array(KernalMatrix)

return KernalMatrix

'''计算g(xi)'''

def Calgxi(self, i):

Index = [index for index, value in enumerate(self.alpha) if value != 0]

gxi = 0

for index in Index:

gxi += self.alpha[index] * self.train_label[index] * self.KernalMatrix[i][index]

gxi = gxi + self.b

return gxi

'''判断是否符合KKT条件'''

def isSatisfyKKT(self, i):

# 获得alpha[i]的值

alpha_i = self.alpha[i]

# 计算yi * g(xi)

gxi = self.Calgxi(i)

yi = self.train_label[i]

yi_gxi = yi * gxi

# 判断是否符合KKT条件

if -1 * self.toler < alpha_i < self.toler and yi_gxi >= 1:

return True

elif -1 * self.toler < alpha_i < self.C + self.toler and math.fabs(yi_gxi - 1) < self.toler:

return True

elif self.C - self.toler < alpha_i < self.C + self.toler and yi_gxi <= 1:

return True

return False

'''SMO算法'''

def SMO(self):

# 迭代

t = 0

parameterchanged = 1

while t < self.itertime and parameterchanged > 0:

t += 1

parameterchanged = 0

'''选择两个alpha'''

# 外层循环，选择第一个alpha

for i in range(self.m):

# 判断是否符合KKT条件，如果不满足，则选择该alpha为alpha1

# 如果满足，则继续外层循环

TorF = self.isSatisfyKKT(i)

if TorF == False:

alpha1 = self.alpha[i]

# 从Earray得到alpha1对应的E1

E1 = self.E[i]

# 复制一个EMatrix，并令E1的位置为nan

# 这样在接下来找最大值和最小值时将不会考虑E1

# 这里需要使用copy，如果不用copy，改变EM_temp也会同时改变EMatrix

EM_temp = np.copy(self.E)

EM_temp[i] = np.nan

# 我们需要使|E1-E2|的值最大，由此选择E2

# 首先初始化maxE1_E2和E2及E2的下标j

maxE1_E2 = -1

E2 = np.nan

j = -1

# 内层循环

# 遍历EM_temp中的每一个Ei，得到使|E1-E2|最大的E和它的下标

for j_temp, Ej in enumerate(EM_temp):

if math.fabs(E1 - Ej) > maxE1_E2:

maxE1_E2 = math.fabs(E1 - Ej)

E2 = Ej

j = j_temp

# alpha2为E2对应的alpha

alpha2 = self.alpha[j]

'''求最优alpha1和alpha2'''

y1 = self.train_label[i]

y2 = self.train_label[j]

# 计算η

K11 = self.KernalMatrix[i][i]

K22 = self.KernalMatrix[j][j]

K12 = self.KernalMatrix[i][j]

eta = K11 + K22 - 2 * K12

# 计算alpha2_new

alpha2_new = alpha2 + y2 * (E1 - E2) / eta

# 计算上限H和下限L

if y1 != y2:

L = max([0, alpha2 - alpha1])

H = min([self.C, self.C + alpha2 - alpha1])

else:

L = max([0, alpha2 + alpha1 - self.C])

H = min([self.C, alpha2 + alpha1])

# 剪切alpha2_new

if alpha2_new > H:

alpha2_new = H

elif alpha2_new < L:

alpha2_new = L

# 得到alpha1_new

alpha1_new = alpha1 + y1 * y2 * (alpha2 - alpha2_new)

'''更新b'''

# 计算b1_new和b2_new

b1_new = -1 * E1 - y1 * K11 * (alpha1_new - alpha1) \

- y2 * K12 * (alpha2_new - alpha2) + self.b

b2_new = -1 * E2 - y1 * K12 * (alpha1_new - alpha1) \

- y2 * K22 * (alpha2_new - alpha2) + self.b

# 根据alpha1和alpha2的范围确定b_new

if 0 < alpha1_new < self.C and 0 < alpha2_new < self.C:

b_new = b1_new

else:

b_new = (b1_new + b2_new) / 2

'''更新E'''

# 首先需要更新两个alpha和b

self.alpha[i] = alpha1_new

self.alpha[j] = alpha2_new

self.b = b_new

# 计算Ei_new和Ej_new

E1_new = self.Calgxi(i) - y1

E2_new = self.Calgxi(j) - y2

# 更新E

self.E[i] = E1_new

self.E[j] = E2_new

if math.fabs(alpha2_new - alpha2) >= 0.000001:

parameterchanged += 1

print('itertime: %a parameterchanged: %a' % (t,parameterchanged))

# 最后遍历一遍alpha，大于0的下标即对应支持向量

VecIndex = [index for index, value in enumerate(self.alpha) if value > 0]

# 返回支持向量的下标，之后在预测时还需要用到

return VecIndex

'''计算b'''

def OptimizeB(self):

for j, a in enumerate(self.alpha):

if 0 < a < self.C:

break

yj = self.train_label[j]

summary = 0

for i in range(self.alpha):

summary += self.alpha[i] * self.train_label[i] * self.KernalMatrix[i][j]

optimiezedB = yj - summary

self.b = optimiezedB

'''计算单个核函数'''

def CalSingleKernal(self, x, z):

SingleKernal = np.exp(-1 * np.dot(x - z, x - z) / (2 * np.square(self.sigma)))

return SingleKernal

'''单个新输入实例的预测'''

def predict(self, x, VecIndex):

# 决策函数计算

# 求和项初始化

summary = 0

# Index中存储着不为0的alpha的下标

for i in VecIndex:

alphai = self.alpha[i]

yi = self.train_label[i]

Kernali = self.CalSingleKernal(x, self.train_data[i])

summary += alphai * yi * Kernali

# 最后+b

# np.sign得到符号

result = np.sign(summary + self.b)

return result

'''测试模型'''

def test(self, test_data, test_label, VecIndex):

# 测试集实例数量

TestNum = len(test_label)

errorCnt = 0

# 对每一个实例进行预测

for i in range(TestNum):

result = self.predict(test_data[i], VecIndex)

if result != test_label[i]:

errorCnt += 1

Acc = 1 - errorCnt / TestNum

return Acc

'''测试模型'''

if __name__ == "__main__":

print('start loading')

# 输入训练集

train_data, train_label = LoadData('E:/我的学习笔记/统计学习方法/实战/训练数据/mnist_train.csv')

# 输入测试集

test_data, test_label = LoadData('E:/我的学习笔记/统计学习方法/实战/训练数据/mnist_test.csv')

print('end loading')

# 初始化

print('start initiating')

start = time.time()

svm = SVM(train_data[0:10000], train_label[0:10000], 10, 200, 0.001, 10)

end = time.time()

print('end initiating')

print('initiating time: ', end - start)

# 开始SMO算法

print('start training')

start = time.time()

VecIndex = svm.SMO()

end = time.time()

print('end training')

print('training time: ', end - start)

# 开始测试模型

Acc = svm.test(test_data[0:1000], test_label[0:1000], VecIndex)

print('Accurate: ', Acc)