原问题:
s.t.
SMO算法的基本思想是:
如果所有变量的解都满足最优化问题的KKT条件,那么这个最优化问题的解就得到了。否则,选择两个变量,固定其他变量,针对这两个变量构建一个二次规划问题。
这样做的目的是,通过求解两个变量的二次规划问题,能不断靠近原有凸二次规划问题的解,并且计算方法有解析方法。
子问题有两个变量,其中一个为违反KKT条件最严重的那一个,另一个由约束条件自动确定。由与约束条件的存在,子问题实际上只有一个自由变量。
整个SMO算法包含两个部分:求解两个变量二次规划的解析方法和选择变量的启发式算法。
两个变量的二次规划求解方法
假设
为变量,其余的量为固定量
设问题的原始可行解为
,最优解为
,并且在沿着约束方向未经剪辑时
的最优解为
由于约束条件存在,因此
的取值范围为:
其中
与
是
所在的对角线段端点的界。
如果
,则
如果
,则
下面先求未经剪辑的最优解
,记
令
定理:
最优化问题沿着约束方向未经剪辑的解时
其中,
其中,
是输入空间到特征空间的映射
经剪辑后
的解,
再由
求得
:
变量的选择方法
SMO算法在每个子问题中需要选择两个变量进行优化,并且至少一个变量是违反KKT条件的。
第一个变量的选择
外层循环,挑选在训练样本中选取违反KKT条件最严重的样本点,并将其对应的变量作为第一个变量。检验训练样本点
是否满足KKT条件,即
其中,
该检验是在
范围内进行的。在检验过程中,外层循环首先遍历所有满足条件
的样本点,如果都满足,则遍历整个数据集。
第二个变量选择
内层循环,假设已经寻找到了第一个变量
,现在需要寻找
,选择标准是希望能使
有足够大的变化。
思路1:由于
是依赖于
的,,我们只需选择
,使其对应的
最大。由于
已经确定,若
是正的,就选择最小的
作为
,如果
是负的,就选择最大的
作为
。
思路2:若思路1所选的
不能使目标函数有足够的下降,则遍历所有在间隔边界上的支持向量点,依次作为
试用,直到目标函数有足够的下降;若还是找不到,则遍历整个数据集;如何还是找不到,则抛弃
,重新寻找新的
。
计算阈值 与差值
每次完成两个变量的优化后,都需要重新计算阈值
,当
时,由KKT条件可知:
于是
用
进行表示,有
同理,如果
,那么
如果
同时满足条件
,那么
;如果
是0或者C,那么
以及他们之间的数都是符合KKT条件的阈值,选择他们的中点作为
。
在每次完成两个变量的优化后,还必须更新对应的
值
其中,S是所有支持向量
的集合。
代码如下:
from sklearn import datasets
import matplotlib.pyplot as plt # 画图工具
import numpy as np
X, y = datasets.make_classification(n_samples=100, n_features=2, n_informative=2, n_redundant=0, n_repeated=0, n_classes=2, n_clusters_per_class=1, random_state=2)
for i in range(X.shape[0]):
if y[i] == 0:
y[i] = -1
#plt.scatter(X[:, 0], X[:, 1], c=y)
# plt.show()
class SMO:
def __init__(self, X, y, C, toler, maxIter): # samples labels constance tolerate maxIteration
self.X = X
self.y = y
self.C = C
self.toler = toler
self.maxIter = maxIter
self.N = self.X.shape[0]
self.b = 0
self.a = np.zeros(self.N)
self.K = np.zeros((self.N, self.N))
self.o = 0.5
for i in range(self.N):
for j in range(self.N):
t = 0
for k in range(self.X.shape[1]):
t += (self.X[i][k] - self.X[j][k])**2
self.K[i][j] = np.exp(-t/(2*self.o**2))
def select_J(self, i): # random choose j which is not equal to i.
j = i
while j == i:
j = np.random.randint(0, self.N)
return j
def fix_alpha(self, a, H, L):
if a > H:
a = H
if a < L:
a = L
return a
def cal_Ei(self, j):
t = 0
for i in range(self.N):
t += self.a[i]*self.y[i]*self.K[i][j]
return t + self.b - self.y[j]
def update(self):
iter = 0
while iter < self.maxIter:
alphaPairsChanged = 0
for i in range(self.N):
Ei = self.cal_Ei(i)
if (Ei < -self.toler and self.a[i] < self.C) or (Ei > self.toler and self.a[i] > 0):
j = self.select_J(i) # choose j != i
Ej = self.cal_Ei(j)
aiold = self.a[i]
ajold = self.a[j]
if self.y[i] != self.y[j]:
L = max(0, self.a[j] - self.a[i])
H = min(self.C, self.C + self.a[j] - self.a[i])
else:
L = max(0, self.a[j] + self.a[i] - self.C)
H = min(self.C, self.a[j] + self.a[i])
if L == H:
print("L==H")
continue
eta = self.K[i][i] + self.K[j][j] - 2*self.K[i][j]
if eta <= 0 :
print('eta <= 0')
continue
ajnew = ajold + self.y[j]*(Ei - Ej)/eta
ajnew = self.fix_alpha(ajnew, H, L)
if abs(ajnew - ajold) < 0.00001:
print("j not moving enough")
continue
self.a[j] = self.fix_alpha(ajnew, H, L)
self.a[i] = aiold + self.y[i]*self.y[j]*(ajold - ajnew)
b1 = -Ei - self.y[i]*self.K[i][i]*(self.a[i] - aiold) - self.y[j]*self.K[i][j]*(ajnew - ajold) + self.b
b2 = -Ej - self.y[i]*self.K[i][j]*(self.a[i] - aiold) - self.y[j]*self.K[j][j]*(ajnew - ajold) + self.b
if 0 < self.a[i] < self.C:
self.b = b1
elif 0 < self.a[j] < self.C:
self.b = b2
else:
self.b = (b1 + b2)/2
alphaPairsChanged += 1
if alphaPairsChanged == 0:
iter += 1
else:
iter = 0
return self.b, self.a
def score(self, X, y):
K = np.zeros((X.shape[0], self.N))
for i in range(X.shape[0]):
for j in range(self.N):
t = 0
for k in range(self.X.shape[1]):
t += (self.X[i][k] - self.X[j][k]) ** 2
K[i][j] = np.exp(-t / (2 * self.o ** 2))
y_pre = []
for i in range(X.shape[0]):
t = 0
for j in range(self.N):
t += self.a[j]*self.y[j]*K[i][j]
y_pre.append(t + self.b)
s = 0
for i in range(X.shape[0]):
if y_pre[i] < 0 and y[i] == -1:
s += 1
elif y_pre[i] > 0 and y[i] == 1:
s += 1
return s/X.shape[0]
L = SMO(X, y, 1, 0.001, 100)
b, a = L.update()
score = L.score(X, y)
print('b = {}'.format(b))
print('a = {}'.format(a))
print('the accuracy of this SVM is {}'.format(score))
运行代码,所得到的结果为:
分别是对偶问题当中的变量
以及最终训练出的SVM模型的精确度为98%