1. alpha2 的修剪
if y1 != y2 :
α1 - α2 = k # 不用算k的具体大小
if k > 0: # 上图的左 下这条线
α2 的区间 (0, c-k)
k < 0 : # 上图的左 下这条线
α2 的区间 （-k, C）
所以：
L = max(0, -k) # k>0 还是<0 都统一表达了
H = min(c, c-k)
else:
y1 = y2 右边的图 同理
2.参考
参考：https://zhuanlan.zhihu.com/p/49331510
非支撑向量：
alpha=0
间隔上的支撑向量：
0<alpha<c
间隔之内：
alpha=c and 0< ei <1
上面的三种都是正确分类的；
错误分类的：
alpha=c and ei > 1
yi(wxi+b) = 1-ei < 0 使负的 预测的与实际的符号相反
3. 代码实现
"""
算法思想：
创建一个alpha向量 并初始化为0
当迭代次数小于 最大迭代次数（外层循环）：
遍历α中的每一个αi（内层循环）：
如果该αi可以被优化：
随机选择另一个αj
同时优化 αi aj
如果αi aj优化完毕退出内层循环
"""
import numpy
import numpy as np
import random
import matplotlib.pyplot as plt
def load_data(filename):
X = []
y = []
fr = open(filename)
for line in fr.readlines():
x1, x2, y_true = line.strip().split("\t")
X.append([float(x1), float(x2)])
y.append(float(y_true))
fr.close()
return numpy.array(X), np.array(y)
def show_data(filename, line=None):
X, y = load_data("./testSet.txt")
class_1_index = np.where(y == -1) # -1 负样本的索引
X_class_1 = X[class_1_index, :].reshape(-1, 2) # 负样本对应的X取出
class_2_index = np.where(y == 1) # +1 正样本的索引
X_class_2 = X[class_2_index, :].reshape(-1, 2) # 正样本对应的X取出
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.scatter(X_class_1[:, 0], X_class_1[:, 1], c="r", label="$-1$")
ax.scatter(X_class_2[:, 0], X_class_2[:, 1], c="g", label="$+1$")
plt.legend()
if line is not None:
alphas, b = line
x = np.linspace(1, 8, 50)
w = np.sum(alphas*y[:, np.newaxis]*X, axis=0)
# w1*x1 + w2*x2 +b = 0 # 或者 +1 -1
y = np.array([(-b - w[0]*x[i]) / w[1]for i in range(50)])
y1 = np.array([(1 -b - w[0]*x[i]) / w[1]for i in range(50)])
y2 = np.array([(-1 -b - w[0]*x[i]) / w[1]for i in range(50)])
ax.plot(x, y, "b-")
ax.plot(x, y1, "b--")
ax.plot(x, y2, "b--")
plt.show()
def smo_simple(dataArr, yArr, C, toler, maxIter):
numSample, numDim = dataArr.shape # 100行 2列
# 初始化
b = 0
alphas = np.zeros((numSample, 1))
iterations = 0 # 迭代次数
while iterations < maxIter:
alphaPairsChanged = 0
for i in range(numSample):
# 针对对i个样本给出预测值 这里没有使用核函数变换
# fXi= wXi+b
# alphas * y * X * Xi +b Xi为一条样本
fXi = np.sum(alphas * yArr[:, np.newaxis] * dataArr * dataArr[i, :]) + b
# Ei 为误差
Ei = fXi - yArr[i]
"""
kkt 约束条件：
yi*(wXi+b) >= 1 且alpha=0 远离边界 正常分类的 ----条件1
yi*(wXi+b) = 1 且0<alpha<c 边界上面
yi*(wXi+b) < 1 且 alpha=c 两条边界之间 ----条件2
# 0 <= alphas[i] <= C 已经在边界上的 alpha不能再优化
Ei = fXi - yArr[i] 两边都 乘以 yArr[i]
yArr[i]*Ei = yArr[i]*fXi - 1 # 偏差
如果 yArr[i]*Ei < -toler alpha应该为C 如果<C 就需要优化 ----对应条件2
如果 yArr[i]*Ei > toler alpha应该为0 如果>0 就需要优化 ----对应条件1
"""
if ( (yArr[i]*Ei < -toler) and (alphas[i] < C) ) or ( (yArr[i]*Ei > toler) and (alphas[i] > 0)):
# 进入这里 说明需要优化alpha[i] 我们在随机取一条alpha[j] j!=i
j = selectJrand(i, numSample)
fXj = np.sum(alphas * yArr[:, np.newaxis] * dataArr * dataArr[j, :]) + b
Ej = fXj - yArr[j]
# 先复制一份old alphaI alphaJ
alphaIold = alphas[i].copy()
alphaJold = alphas[j].copy()
# 根据 yArr[i] yArr[j] 是否同号 计算 L H
if yArr[i] != yArr[j]:
L = max(0, alphas[j] - alphas[i])
H = min(C, C + alphas[j] - alphas[i])
else:
L = max(0, alphas[j] + alphas[i] - C)
H = min(C, alphas[j] + alphas[i])
if L == H:
print("L == H")
continue
# 计算 eta, eta是alpha[j]的最大修正量 eta == 0 则 continue
# eta = K(1,1)+K(2,2) - 2K(1,2)
eta = (np.sum(dataArr[i,:] * dataArr[i, :]) +
np.sum(dataArr[j, :] * dataArr[j, :]) -
2 * np.sum(dataArr[i,:] * dataArr[j, :]))
if eta == 0:
continue
# 计算新的alphas[j]
alphas[j] = alphaJold + yArr[j]*(Ei - Ej) / eta
# 修剪新的alphas[j]
alphas[j] = clipAlpha(alphas[j], H, L)
# 新的alphas[j] 与 alphaJold 如果改动量很小 则 contine 退出
if abs(alphaJold - alphas[j]) < 0.00001:
continue
# 更新 alphas[i]
alphas[i] = alphaIold + yArr[i] * yArr[j] * (alphaJold - alphas[j])
# 计算参数b
bi = b - Ei - yArr[i] * (alphas[i] - alphaIold) * np.sum(dataArr[i, :] * dataArr[i, :]) - \
yArr[j] * (alphas[j] - alphaJold) * np.sum(dataArr[i, :] * dataArr[j, :])
bj = b - Ej - yArr[i] * (alphas[i] - alphaIold) * np.sum(dataArr[i, :] * dataArr[j, :]) - \
yArr[j] * (alphas[j] - alphaJold) * np.sum(dataArr[j, :] * dataArr[j, :])
if 0<alphas[i]<C:
b = bi
elif 0<alphas[j]<C:
b = bj
else:
b = (bi+bj)/2
# 走到这里说明 alpha b 更新了
alphaPairsChanged += 1
# 输出
print(f"iter: {iterations}, i: {i}, pairs changed: {alphaPairsChanged}")
if alphaPairsChanged ==0:
iterations += 1
else:
iterations =0
return b, alphas
def selectJrand(i, numSample):
j = i
while j == i:
j = int(random.uniform(0, numSample))
return j
def clipAlpha(aj, H, L):
if aj > H:
return H
if aj < L:
return L
return aj
if __name__ == '__main__':
X, y = load_data("./testSet.txt")
C = 0.6
toler = 0.001
maxIter = 40
b, alphas = smo_simple(X, y, C, toler, maxIter)
print(b)
print(alphas)
show_data("./testSet.txt", line=(alphas, b))
