2.完整版的SMO算法
在几百个点组成的小规模数据集上,简化版SMO算法的运行是没有什么问题的,但是在更大的数据集上的运行速度就会变慢。
刚才已经讨论了简化版SMO算法,下面我们就讨论完整版的Platt SMO算法。
在这两个版本中,实现alpha的更改和代数运算的优化环节一模一样。在优化过程中,唯一的不同就是选择alpha的方式。
完整版的Platt SMO算法应用了一些能够提速的启发方法。上一节的例子在执行时存在一定的时间提升空间。
Platt SMO算法是通过一个外循环来选择第一个alpha值的,并且其选择过程会在两种方式之间进行交替:一种方式是在所有数据集上进行单遍扫描,另一种方式则是在非边界alpha中实现单遍扫描。
而所谓非边界alpha指的就是那些不等于边界0或C的alpha值。
对整个数据集的扫描相当容易,而实现非边界alpha值的扫描时,首先需要建立这些alpha值的列表,然后再对这个表进行遍历。同时,该步骤会跳过那些已知的不会改变的alpha值。
在选择第一个alpha值后,算法会通过一个内循环来选择第二个alpha值。在优化过程中,会通过最大化步长的方式来获得第二个alpha值。
在简化版SMO算法中,我们会在选择j之后计算错误率Ej。但在这里,我们会建立一个全局的缓存用于保存误差值,并从中选择使得步长或者说Ei-Ej最大的alpha值。
辅助函数和上面的步骤一样,即下面的三个函数:
"""
函数说明:读取数据
Parameters:
file_name - 文件名
Returns:
data_mat - 数据矩阵
label_mat - 数据标签
"""
def load_dataset(file_name):
# 数据矩阵
data_mat = []
# 标签向量
label_mat = []
# 打开文件
fr = open(file_name)
# 逐行读取
for line in fr.readlines():
# 去掉每一行首尾的空白符,例如'\n','\r','\t',' '
# 将每一行内容根据'\t'符进行切片
line_array = line.strip().split('\t')
# 添加数据(100个元素排成一行)
data_mat.append([float(line_array[0]), float(line_array[1])])
# 添加标签(100个元素排成一行)
label_mat.append(float(line_array[2]))
return data_mat, label_mat
"""
函数说明:随机选择alpha_j
Parameters:
i - alpha_i的索引值
m - alpha参数个数
Returns:
j - alpha_j的索引值
"""
def select_j_random(i, m):
j = i
while(j == i):
# uniform()方法将随机生成一个实数,它在[x, y)范围内
j = int(random.uniform(0, m))
return j
"""
函数说明:修剪alpha_j
Parameters:
aj - alpha_j值
H - alpha上限
L - alpha下限
Returns:
aj - alpha_j值
"""
def clip_alpha(aj, H, L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
我们首先定义一个数据结构:
"""
类说明:维护所有需要操作的值
Parameters:
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
Returns:
None
"""
"""
定义一个新的数据结构
"""
class opt_struct:
def __init__(self, data_mat_in, class_labels, C, toler):
# 数据矩阵
self.X = data_mat_in #传进来的数据
# 数据标签
self.label_mat = class_labels
# 松弛变量
self.C = C
# 容错率
self.tol = toler
# 矩阵的行数
self.m = np.shape(data_mat_in)[0]
# 根据矩阵行数初始化alphas矩阵,一个m行1列的全零列向量
self.alphas = np.mat(np.zeros((self.m, 1)))
# 初始化b参数为0
self.b = 0
# 根据矩阵行数初始化误差缓存矩阵,第一列为是否有效标志位,其中0无效,1有效;第二列为实际的误差Ei的值
"""
我们之前的定义为:Ei=gxi-yi
yi是标签的实际值。
gx=alpha_i*y_i*x_i.x,就相当于是w.x+b
因为误差值经常用到,所以希望每次计算后放到一个缓存当中,将ecache一分为二,第一列是标志位,取值为0或者1,为1时表示已经算出来
"""
self.ecache = np.mat(np.zeros((self.m, 2)))
首要的事情就是建立一个数据结构来保存所有的重要值,而这个过程可以通过一个对象来完成。只是作为一个数据结构来使用对象。
在将值传给函数时,我们可以通过将所有数据移到一个结构中来实现,这样就可以省掉手工输入的麻烦了。而此时,数据就可以通过一个对象来进行传递。实际上,当完成其实现时,可以很容易通过Python的字典来完成。
该方法可以实现其成员变量的填充。除了增加了一个m×2的矩阵成
员变量ecache之外,这些做法和简化版SMO一模一样。ecache的第一列给出的是ecache是否有效的标志位,而第二列给出的是实际的E值。
对于给定的alpha值,第一个辅助函数cal-Ek()能够计算E值并返回。以前,该过程是采用内嵌的方式来完成的,但是由于该过程在这个版本的SMO算法中出现频繁,这里必须要将其单独拎出来。
"""
函数说明:计算误差
Parameters:
os - 数据结构
k - 标号为k的数据
Returns:
Ek - 标号为k的数据误差
"""
def cal_Ek(os, k):
# multiply(a,b)就是个乘法,如果a,b是两个数组,那么对应元素相乘
# .T为转置
fXk = float(np.multiply(os.alphas, os.label_mat).T * (os.X * os.X[k, :].T) + os.b)
# 计算误差项
Ek = fXk - float(os.label_mat[k])
# 返回误差项
return Ek
函数select-j()用于选择第二个alpha或者说内循环的alpha值。回想一下,这里的目标是选择合适的第二个alpha值以保证在每次优化中采用最大步长。该函数的误差值与第一个alpha值Ei和下标i有关。
首先将输入值Ei在缓存中设置成为有效的。这里的有效意味着它已经计算好了。在ecache中,代码nonzero(os.ecache[:,0].A)[0]构建出了一个非零表。
NumPy函数nonzero()返回了一个列表,而这个列表中包含以输入列表为目录的列表值,当然这里的值并非零。nonzero()语句返回的是非零E值所对应的alpha值,而不是E值本身。
程序会在所有的值上进行循环并选择其中使得改变最大的那个值。如果这是第一次循环的话,那么就随机选择一个alpha值。当然也存在有许多更复杂的方式来处理第一次循环的情况,而上述做法就能够满足我们的目的。
"""
函数说明:内循环启发方式2
选择第二个待优化的alpha_j,选择一个误差最大的alpha_j
即,我们在选择alpha_2的时候做了改进,选择误差最大的
Parameters:
i - 标号为i的数据的索引值
oS - 数据结构
Ei - 标号为i的数据误差
Returns:
j - 标号为j的数据的索引值
maxK - 标号为maxK的数据的索引值
Ej - 标号为j的数据误差
"""
def select_j(i, os, Ei):
# 初始化
max_K = -1 #下标的索引值
max_delta_E = 0
Ej = 0
# 根据Ei更新误差缓存,即先计算alpha_1以及E1值
os.ecache[i] = [1, Ei] #放入缓存当中,设为有效
# 对一个矩阵.A转换为Array类型
# 返回误差不为0的数据的索引值
valid_ecache_list = np.nonzero(os.ecache[:, 0].A)[0] #找出缓存中不为0
# 有不为0的误差
if(len(valid_ecache_list) > 1):
# 遍历,找到最大的Ek
for k in valid_ecache_list: #迭代所有有效的缓存,找到误差最大的E
# 不计算k==i节省时间
if k == i: #不选择和i相等的值
continue
# 计算Ek
Ek = cal_Ek(os, k)
# 计算|Ei - Ek|
delta_E = abs(Ei - Ek)
# 找到maxDeltaE
if(delta_E > max_delta_E):
max_K = k
max_delta_E = delta_E
Ej = Ek
# 返回max_K,Ej
return max_K, Ej #这样我们就得到误差最大的索引值和误差最大的值
# 没有不为0的误差
else: #第一次循环时是没有有效的缓存值的,所以随机选一个(仅会执行一次)
# 随机选择alpha_j的索引值
j = select_j_random(i, os.m)
# 计算Ej
Ej = cal_Ek(os, j)
# 返回j,Ej
return j, Ej
最后一个辅助函数是update-Ek(),它会计算误差值并存入缓存当中。在对alpha值进行优化之后会用到这个值。
"""
函数说明:计算Ek,并更新误差缓存
Parameters:
os - 数据结构
k - 标号为k的数据的索引值
Returns:
None
"""
def update_Ek(os, k):
# 计算Ek
Ek = cal_Ek(os, k)
# 更新误差缓存
os.ecache[k] = [1, Ek]
接下来将简单介绍一下用于寻找决策边界的优化例程:
"""
函数说明:优化的SMO算法
Parameters:
i - 标号为i的数据的索引值
os - 数据结构
Returns:
1 - 有任意一对alpha值发生变化
0 - 没有任意一对alpha值发生变化或变化太小
"""
def innerL(i, os):
# 步骤1:计算误差Ei
Ei = cal_Ek(os, i)
# 优化alpha,设定一定的容错率
if((os.label_mat[i] * Ei < -os.tol) and (os.alphas[i] < os.C)) or ((os.label_mat[i] * Ei > os.tol) and (os.alphas[i] > 0)):
# 使用内循环启发方式2选择alpha_j,并计算Ej
j, Ej = select_j(i, os, Ei) #这里不再是随机选取了
# 保存更新前的alpha值,使用深层拷贝
alpha_i_old = os.alphas[i].copy()
alpha_j_old = os.alphas[j].copy()
# 步骤2:计算上界H和下界L
if(os.label_mat[i] != os.label_mat[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
# 步骤3:计算eta
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
# 步骤4:更新alpha_j
os.alphas[j] -= os.label_mat[j] * (Ei - Ej) / eta
# 步骤5:修剪alpha_j
os.alphas[j] = clip_alpha(os.alphas[j], H, L)
# 更新Ej至误差缓存
update_Ek(os, j)
if(abs(os.alphas[j] - alpha_j_old) < 0.00001):
print("alpha_j变化太小")
return 0
# 步骤6:更新alpha_i
os.alphas[i] += os.label_mat[i] * os.label_mat[j] * (alpha_j_old - os.alphas[j])
# 更新Ei至误差缓存
update_Ek(os, i)
# 步骤7:更新b_1和b_2:
b1 = os.b - Ei - os.label_mat[i] * (os.alphas[i] - alpha_i_old) * os.X[i, :] * os.X[i, :].T - os.label_mat[j] * (os.alphas[j] - alpha_j_old) * os.X[j, :] * os.X[i, :].T
b2 = os.b - Ej - os.label_mat[i] * (os.alphas[i] - alpha_i_old) * os.X[i, :] * os.X[j, :].T - os.label_mat[j] * (os.alphas[j] - alpha_j_old) * os.X[j, :] * os.X[j, :].T
# 步骤8:根据b_1和b_2更新b
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:
return 0 #表示没有更新
代码几乎与上一节中给出的smo-simple函数一摸一样,但是这里的代码已经使用了自己的数据结构。该结构在参数os中传递。
完整版SMO算法的外循环代码:
"""
函数说明:完整的线性SMO算法
Parameters:
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
maxIter - 最大迭代次数
Returns:
oS.b - SMO算法计算的b
oS.alphas - SMO算法计算的alphas
"""
def smo_p(data_mat_in, class_labels, C, toler, max_iter):
# 初始化数据结构
os = opt_struct(np.mat(data_mat_in), np.mat(class_labels).transpose(), C, toler)
# 初始化当前迭代次数
iter = 0
entrie_set = True #是否在全部数据集上迭代
alpha_pairs_changed = 0
# 遍历整个数据集alpha都没有更新或者超过最大迭代次数,则退出循环
while(iter < max_iter) and ((alpha_pairs_changed > 0) or (entrie_set)):
alpha_pairs_changed = 0
if entrie_set: # 遍历整个数据集
for i in range(os.m):
# 使用优化的SMO算法
alpha_pairs_changed += innerL(i, os) #innerL返回的值是0或者1
print("全样本遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter, i, alpha_pairs_changed))
iter += 1
# 遍历非边界值
else:
# 遍历不在边界0和C的alpha
non_bound_i_s = np.nonzero((os.alphas.A > 0) * (os.alphas.A < C))[0]
for i in non_bound_i_s:
alpha_pairs_changed += innerL(i, os)
print("非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter, i, alpha_pairs_changed))
iter += 1
# 遍历一次后改为非边界遍历
if entrie_set:
entrie_set = False #进行切换,遍历非边界数据集
# 如果alpha没有更新,计算全样本遍历
elif(alpha_pairs_changed == 0):
entrie_set = True
print("迭代次数:%d" % iter)
# 返回SMO算法计算的b和alphas
return os.b, os.alphas
其输入和函数smo-simple()完全一样。函数一开始构建一个数据结构来容纳所有的数据,然后需要对控制函数退出的一些变量进行初始化。整个代码的主体是while循环,这与smo-simple()有些类似,但是这里的循环退出条件更多一些。
当迭代次数超过指定的最大值,或者遍历整个集合都未对任意alpha对进行修改时,就退出循环。
这里的max-iter变量和函数smo-simple()中的作用有一点不同,后者当没有任何alpha发生改变时会将整个集合的一次遍历过程计成一次迭代,而这里的一次迭代定义为一次循环过程,而不管该循环具体做了什么事。此时,如果在优化过程中存在波动就会停止,因此这里的做法优于smo-simple()函数中的计数方法。
while循环的内部与smo-simple()中有所不同,一开始的for循环在数据集上遍历任意可能的alpha。
我们通过调用innerL()来选择第二个alpha,并在可能时对其进行优化处理。如果有任意一对alpha值发生改变,那么会返回
1。第二个for循环遍历所有的非边界alpha值,也就是不在边界0或C上的值。
接下来,我们对for循环在非边界循环和完整遍历之间进行切换,并打印出迭代次数。最后程序将会返回常数b和alpha值。
我们的部分结果如下:
全样本遍历:第0次迭代 样本:0, alpha优化次数:1
全样本遍历:第0次迭代 样本:1, alpha优化次数:1
全样本遍历:第0次迭代 样本:2, alpha优化次数:1
全样本遍历:第0次迭代 样本:3, alpha优化次数:1
全样本遍历:第0次迭代 样本:4, alpha优化次数:1
全样本遍历:第0次迭代 样本:5, alpha优化次数:1
全样本遍历:第0次迭代 样本:6, alpha优化次数:1
全样本遍历:第0次迭代 样本:7, alpha优化次数:1
全样本遍历:第0次迭代 样本:8, alpha优化次数:2
全样本遍历:第0次迭代 样本:9, alpha优化次数:2
......
非边界遍历:第1次迭代 样本:26, alpha优化次数:3
alpha_j变化太小
非边界遍历:第1次迭代 样本:29, alpha优化次数:3
alpha_j变化太小
非边界遍历:第1次迭代 样本:46, alpha优化次数:3
alpha_j变化太小
非边界遍历:第1次迭代 样本:54, alpha优化次数:3
非边界遍历:第1次迭代 样本:55, alpha优化次数:3
迭代次数:2
alpha_j变化太小
非边界遍历:第2次迭代 样本:8, alpha优化次数:0
alpha_j变化太小
非边界遍历:第2次迭代 样本:10, alpha优化次数:0
非边界遍历:第2次迭代 样本:23, alpha优化次数:0
alpha_j变化太小
非边界遍历:第2次迭代 样本:29, alpha优化次数:0
alpha_j变化太小
非边界遍历:第2次迭代 样本:46, alpha优化次数:0
alpha_j变化太小
非边界遍历:第2次迭代 样本:54, alpha优化次数:0
非边界遍历:第2次迭代 样本:55, alpha优化次数:0
......
全样本遍历:第3次迭代 样本:97, alpha优化次数:0
全样本遍历:第3次迭代 样本:98, alpha优化次数:0
全样本遍历:第3次迭代 样本:99, alpha优化次数:0
迭代次数:4
我们最终的决策边界和支持向量可视化为:
b= [[-3.06127512]]
w= [[ 0.69141217]
[-0.3411598 ]]
alphas= [[0.06158536 0.08744873 0.09509032 0.07285032 0.01603246 0.15524162]]
[3.542485, 1.977398] -1.0
[3.125951, 0.293251] -1.0
[4.658191, 3.507396] -1.0
[8.197181, 1.545132] 1.0
[5.286862, -2.358286] 1.0
[6.080573, 0.418886] 1.0
完整版代码如下:
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import random
import matplotlib as mpl
"""
改进SMO算法以加快我们的SVM运行速度!
"""
"""
函数说明:读取数据
Parameters:
file_name - 文件名
Returns:
data_mat - 数据矩阵
label_mat - 数据标签
"""
def load_dataset(file_name):
# 数据矩阵
data_mat = []
# 标签向量
label_mat = []
# 打开文件
fr = open(file_name)
# 逐行读取
for line in fr.readlines():
# 去掉每一行首尾的空白符,例如'\n','\r','\t',' '
# 将每一行内容根据'\t'符进行切片
line_array = line.strip().split('\t')
# 添加数据(100个元素排成一行)
data_mat.append([float(line_array[0]), float(line_array[1])])
# 添加标签(100个元素排成一行)
label_mat.append(float(line_array[2]))
return data_mat, label_mat
"""
函数说明:随机选择alpha_j
Parameters:
i - alpha_i的索引值
m - alpha参数个数
Returns:
j - alpha_j的索引值
"""
def select_j_random(i, m):
j = i
while(j == i):
# uniform()方法将随机生成一个实数,它在[x, y)范围内
j = int(random.uniform(0, m))
return j
"""
函数说明:修剪alpha_j
Parameters:
aj - alpha_j值
H - alpha上限
L - alpha下限
Returns:
aj - alpha_j值
"""
def clip_alpha(aj, H, L):
if aj > H:
aj = H
if L > aj:
aj = L
return aj
"""
类说明:维护所有需要操作的值
Parameters:
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
Returns:
None
"""
"""
定义一个新的数据结构
"""
class opt_struct:
def __init__(self, data_mat_in, class_labels, C, toler):
# 数据矩阵
self.X = data_mat_in #传进来的数据
# 数据标签
self.label_mat = class_labels
# 松弛变量
self.C = C
# 容错率
self.tol = toler
# 矩阵的行数
self.m = np.shape(data_mat_in)[0]
# 根据矩阵行数初始化alphas矩阵,一个m行1列的全零列向量
self.alphas = np.mat(np.zeros((self.m, 1)))
# 初始化b参数为0
self.b = 0
# 根据矩阵行数初始化误差缓存矩阵,第一列为是否有效标志位,其中0无效,1有效;第二列为实际的误差Ei的值
"""
我们之前的定义为:Ei=gxi-yi
yi是标签的实际值。
gx=alpha_i*y_i*x_i.x,就相当于是w.x+b
因为误差值经常用到,所以希望每次计算后放到一个缓存当中,将ecache一分为二,第一列是标志位,取值为0或者1,为1时表示已经算出来
"""
self.ecache = np.mat(np.zeros((self.m, 2)))
"""
函数说明:计算误差
Parameters:
os - 数据结构
k - 标号为k的数据
Returns:
Ek - 标号为k的数据误差
"""
def cal_Ek(os, k):
# multiply(a,b)就是个乘法,如果a,b是两个数组,那么对应元素相乘
# .T为转置
fXk = float(np.multiply(os.alphas, os.label_mat).T * (os.X * os.X[k, :].T) + os.b)
# 计算误差项
Ek = fXk - float(os.label_mat[k])
# 返回误差项
return Ek
"""
函数说明:内循环启发方式2
选择第二个待优化的alpha_j,选择一个误差最大的alpha_j
即,我们在选择alpha_2的时候做了改进,选择误差最大的
Parameters:
i - 标号为i的数据的索引值
oS - 数据结构
Ei - 标号为i的数据误差
Returns:
j - 标号为j的数据的索引值
maxK - 标号为maxK的数据的索引值
Ej - 标号为j的数据误差
"""
def select_j(i, os, Ei):
# 初始化
max_K = -1 #下标的索引值
max_delta_E = 0
Ej = 0
# 根据Ei更新误差缓存,即先计算alpha_1以及E1值
os.ecache[i] = [1, Ei] #放入缓存当中,设为有效
# 对一个矩阵.A转换为Array类型
# 返回误差不为0的数据的索引值
valid_ecache_list = np.nonzero(os.ecache[:, 0].A)[0] #找出缓存中不为0
# 有不为0的误差
if(len(valid_ecache_list) > 1):
# 遍历,找到最大的Ek
for k in valid_ecache_list: #迭代所有有效的缓存,找到误差最大的E
# 不计算k==i节省时间
if k == i: #不选择和i相等的值
continue
# 计算Ek
Ek = cal_Ek(os, k)
# 计算|Ei - Ek|
delta_E = abs(Ei - Ek)
# 找到maxDeltaE
if(delta_E > max_delta_E):
max_K = k
max_delta_E = delta_E
Ej = Ek
# 返回max_K,Ej
return max_K, Ej #这样我们就得到误差最大的索引值和误差最大的值
# 没有不为0的误差
else: #第一次循环时是没有有效的缓存值的,所以随机选一个(仅会执行一次)
# 随机选择alpha_j的索引值
j = select_j_random(i, os.m)
# 计算Ej
Ej = cal_Ek(os, j)
# 返回j,Ej
return j, Ej
"""
函数说明:计算Ek,并更新误差缓存
Parameters:
os - 数据结构
k - 标号为k的数据的索引值
Returns:
None
"""
def update_Ek(os, k):
# 计算Ek
Ek = cal_Ek(os, k)
# 更新误差缓存
os.ecache[k] = [1, Ek]
"""
函数说明:优化的SMO算法
Parameters:
i - 标号为i的数据的索引值
os - 数据结构
Returns:
1 - 有任意一对alpha值发生变化
0 - 没有任意一对alpha值发生变化或变化太小
"""
def innerL(i, os):
# 步骤1:计算误差Ei
Ei = cal_Ek(os, i)
# 优化alpha,设定一定的容错率
if((os.label_mat[i] * Ei < -os.tol) and (os.alphas[i] < os.C)) or ((os.label_mat[i] * Ei > os.tol) and (os.alphas[i] > 0)):
# 使用内循环启发方式2选择alpha_j,并计算Ej
j, Ej = select_j(i, os, Ei) #这里不再是随机选取了
# 保存更新前的alpha值,使用深层拷贝
alpha_i_old = os.alphas[i].copy()
alpha_j_old = os.alphas[j].copy()
# 步骤2:计算上界H和下界L
if(os.label_mat[i] != os.label_mat[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
# 步骤3:计算eta
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
# 步骤4:更新alpha_j
os.alphas[j] -= os.label_mat[j] * (Ei - Ej) / eta
# 步骤5:修剪alpha_j
os.alphas[j] = clip_alpha(os.alphas[j], H, L)
# 更新Ej至误差缓存
update_Ek(os, j)
if(abs(os.alphas[j] - alpha_j_old) < 0.00001):
print("alpha_j变化太小")
return 0
# 步骤6:更新alpha_i
os.alphas[i] += os.label_mat[i] * os.label_mat[j] * (alpha_j_old - os.alphas[j])
# 更新Ei至误差缓存
update_Ek(os, i)
# 步骤7:更新b_1和b_2:
b1 = os.b - Ei - os.label_mat[i] * (os.alphas[i] - alpha_i_old) * os.X[i, :] * os.X[i, :].T - os.label_mat[j] * (os.alphas[j] - alpha_j_old) * os.X[j, :] * os.X[i, :].T
b2 = os.b - Ej - os.label_mat[i] * (os.alphas[i] - alpha_i_old) * os.X[i, :] * os.X[j, :].T - os.label_mat[j] * (os.alphas[j] - alpha_j_old) * os.X[j, :] * os.X[j, :].T
# 步骤8:根据b_1和b_2更新b
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:
return 0 #表示没有更新
"""
函数说明:完整的线性SMO算法
Parameters:
dataMatIn - 数据矩阵
classLabels - 数据标签
C - 松弛变量
toler - 容错率
maxIter - 最大迭代次数
Returns:
oS.b - SMO算法计算的b
oS.alphas - SMO算法计算的alphas
"""
def smo_p(data_mat_in, class_labels, C, toler, max_iter):
# 初始化数据结构
os = opt_struct(np.mat(data_mat_in), np.mat(class_labels).transpose(), C, toler)
# 初始化当前迭代次数
iter = 0
entrie_set = True #是否在全部数据集上迭代
alpha_pairs_changed = 0
# 遍历整个数据集alpha都没有更新或者超过最大迭代次数,则退出循环
while(iter < max_iter) and ((alpha_pairs_changed > 0) or (entrie_set)):
alpha_pairs_changed = 0
if entrie_set: # 遍历整个数据集
for i in range(os.m):
# 使用优化的SMO算法
alpha_pairs_changed += innerL(i, os) #innerL返回的值是0或者1
print("全样本遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter, i, alpha_pairs_changed))
iter += 1
# 遍历非边界值
else:
# 遍历不在边界0和C的alpha
non_bound_i_s = np.nonzero((os.alphas.A > 0) * (os.alphas.A < C))[0]
for i in non_bound_i_s:
alpha_pairs_changed += innerL(i, os)
print("非边界遍历:第%d次迭代 样本:%d, alpha优化次数:%d" % (iter, i, alpha_pairs_changed))
iter += 1
# 遍历一次后改为非边界遍历
if entrie_set:
entrie_set = False #进行切换,遍历非边界数据集
# 如果alpha没有更新,计算全样本遍历
elif(alpha_pairs_changed == 0):
entrie_set = True
print("迭代次数:%d" % iter)
# 返回SMO算法计算的b和alphas
return os.b, os.alphas
"""
函数说明:分类结果可视化
Returns:
dataMat - 数据矩阵
classLabels - 数据标签
w - 直线法向量
b - 直线截距
Returns:
None
"""
def show_classifer(data_mat, class_labels, w, b):
data_mat, label_mat = load_dataset('testSet.txt')
fig = plt.figure()
ax = fig.add_subplot(111)
cm = mpl.colors.ListedColormap(['g', 'b'])
ax.scatter(np.array(data_mat)[:, 0], np.array(data_mat)[:, 1], c=np.array(label_mat), cmap=cm, s=20)
# 绘制直线
x1 = max(data_mat)[0]
x2 = min(data_mat)[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])
# 找出支持向量点
# enumerate在字典上是枚举、列举的意思
for i, alpha in enumerate(alphas):
# 支持向量机的点
if (abs(alpha) > 0):
x, y = data_mat[i]
plt.scatter([x], [y], s=150, c='none', alpha=0.7, linewidth=1.5, edgecolors='red')
plt.show()
"""
函数说明:计算w
Returns:
dataArr - 数据矩阵
classLabels - 数据标签
alphas - alphas值
Returns:
w - 直线法向量
"""
def cal_w_s(alphas, data_array, class_labels):
X = np.mat(data_array)
label_mat = np.mat(class_labels).transpose()
m, n = np.shape(X)
w = np.zeros((n, 1))
for i in range(m):
w += np.multiply(alphas[i] * label_mat[i], X[i, :].T)
return w
if __name__ == '__main__':
data_array, class_labels = load_dataset('testSet.txt')
b, alphas = smo_p(data_array, class_labels, 0.6, 0.001, 40)
w = cal_w_s(alphas, data_array, class_labels)
print('b=', b)
print('w=', w)
print('alphas=', alphas[alphas > 0])
for i in range(100):
if alphas[i] > 0.0:
print(data_array[i], class_labels[i])
show_classifer(data_array, class_labels, w, b)