占据统计学习方法大半江山的SVM，终于在昨天以我的代码写完，画上了一个句号，不得不说，支持向量机不是一个简单的东西，尽管花了很大的力气去做，在最后，也没有将所有的细节都理解了，在这，我准备总结一些支持向量机，也是对前面几篇理解不正确的地方的一个修正。
0X01:算法
从大体上看，给我们一组数据，我们怎么才能通过支持向量机训练模型来拟合这些数据呢？
确定使用什么核技巧，是线性核，还是多项式核，这个是取决于你数据是不是非线性可分的。
确定最大迭代次数，确定惩罚参数C的大小，与软间隔硬间隔宽窄有关。
使用SMO算法确定要优化的两个拉格朗日乘子，通过检验样本点是否满足KKT条件确定第一个乘子，通过等式约束条件确定第二个乘子。
根据选定的两个乘子对应的label是否相同，求出拉格朗日乘子的上下界。
通过预测函数
求出误差值
:
再通过data，计算出,也就是说：
最终得到的未修正值:
-
通过上面得到的上下界，修剪
:
def compare(self, alpha_not_clip, L, H): if alpha_not_clip < L: return L elif alpha_not_clip > H: return H else: return alpha_not_clip 通过
:
这个公式是这么来的：
-
根据
得到
:
同理
这个公式是从这里来的：
再用E的定义式替换，就得到了最终的b 更新E的值：
直到迭代结束
0x02 不能理解的地方
尽管看了很多内容，仍然对一些地方理解不了：
- 判断样本点是否满足KKT条件的方式
0x03 补充
见好多人都提到了坐标上升法，我也在这里说一说：
坐标上升法确实与SMO确实很相似，但是有一些很细微，也很重要的区别：
-
坐标上升法只优化一个变量
SMO优化两个变量。
这是为什么，SMO为什么不能优化一个变量？
因为：
SMO比坐标上升法多出来一些约束条件，所以，如果只优化一个变量，就可以根据其他的“常量”将那个变量推断出来。
还有KKT条件的一些理解：
KKT条件其实是在有约束的优化问题上（不等式约束上）推到出来的，最优解应该满足的条件。
具体参考这篇blog:浅谈最优化问题的KKT条件
0x04 代码
代码参考了小伙伴Tenshine的代码，但也对他的代码做了一些改进，将他代码中的循环全部以矩阵运算的方式展开了，加快了训练速度。（矩阵运算与循环的差异可以参考这个回答：在机器学习的算法中， 矩阵运算的效率总是优于for循环吗？)
import numpy as np
import pandas as pd
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
def create_data():
iris = load_iris()
df = pd.DataFrame(iris.data, columns=iris.feature_names)
df['label'] = iris.target
df.columns = ['sepal length', 'sepal width',
'petal length', 'petal width', 'label']
data = np.array(df.iloc[:100, [0, 1, -1]])
for i in range(len(data)):
if data[i, -1] == 0:
data[i, -1] = -1
return data[:, :2], data[:, -1]
class SVM:
def __init__(self, data, max_iters=1000, kernal="linear", C=0.0001):
self.max_iters = max_iters
self.Kernal = kernal
self.load_data(data[0], data[1])
self.C = C
def load_data(self, X, y):
self.X = X
self.Y = y
def init_args(self):
self.data_num = self.X.shape[0]
self.dim = self.X.shape[1]
self.alpha = np.zeros(self.data_num)
self.b = 0
self.f = [self.com_f(i) for i in range(self.data_num)]
self.E = [self.com_E(i) for i in range(self.data_num)]
def com_f(self, idx):
return self.kernal(np.multiply(self.alpha, self.Y), np.dot(self.X[idx, :], self.X.T)) + self.b
def com_E(self, idx):
return self.f[idx] - self.Y[idx]
def com_eta(self, i, j):
return self.kernal(self.X[i], self.X[i]) + self.kernal(self.X[j], self.X[j]) - 2*self.kernal(self.X[i], self.X[j])
def search_sample(self):
l1 = [i for i in range(self.data_num) if 0 < self.alpha[i] < self.C]
l2 = [i for i in range(self.data_num) if i not in l1]
l1.extend(l2)
for i in l1:
if self.kkt(i):
continue
E = self.E[i]
if E > 0:
j = min(range(self.data_num), key=lambda x: self.E[x])
else:
j = max(range(self.data_num), key=lambda x: self.E[x])
return i, j
def kkt(self, idx):
judge = self.Y[idx]*self.f[idx]
if self.alpha[idx] == 0:
return judge >= 1
elif self.alpha[idx] > 0 and self.alpha[idx] < self.C:
return judge == 1
else:
return judge <= 1
def compare(self, alpha_not_clip, L, H):
if alpha_not_clip < L:
return L
elif alpha_not_clip > H:
return H
else:
return alpha_not_clip
def kernal(self, x_i, x_j):
if self.Kernal == "linear":
return np.dot(x_i, x_j)
elif self.Kernal == "poly":
return pow(np.dot(x_i, x_j)+1, 2)
def train(self):
for _ in range(self.max_iters):
i, j = self.search_sample()
print(i,",",j)
eta = self.com_eta(i, j)
# print(eta)
if self.Y[i] == self.Y[j]:
L = max(0, self.alpha[i]+self.alpha[j] - self.C)
H = min(self.C, self.alpha[i]+self.alpha[j])
else:
L = max(0, self.alpha[j]-self.alpha[i])
H = min(self.C, self.C + self.alpha[j]-self.alpha[i])
E1 = self.E[i]
E2 = self.E[j]
if eta == 0:
continue
alpha2_new_unclip = self.alpha[j] + self.Y[j]*(E1 - E2) / eta
alpha2 = self.compare(alpha2_new_unclip, L, H)
alpha1 = self.alpha[i] + self.Y[i] * self.Y[j] * (self.alpha[j] - alpha2)
b1_new = -E1 - self.Y[i] * self.kernal(self.X[i], self.X[i]) * (
alpha1-self.alpha[i]) - self.Y[j] * self.kernal(self.X[j], self.X[i]) * (alpha2-self.alpha[j]) + self.b
b2_new = -E2 - self.Y[i] * self.kernal(self.X[i], self.X[j]) * (
alpha1-self.alpha[i]) - self.Y[j] * self.kernal(self.X[j], self.X[j]) * (alpha2-self.alpha[j]) + self.b
if 0 < alpha1 < self.C:
b_new = b1_new
elif 0 < alpha2 < self.C:
b_new = b2_new
else:
b_new = (b1_new + b2_new) / 2
self.alpha[i] = alpha1
self.alpha[j] = alpha2
self.b = b_new
self.E[i] = self.com_E(i)
self.E[j] = self.com_E(j)
print('train done!')
def predict(self, x):
r = self.kernal(np.multiply(self.alpha, self.Y), np.dot(x, self.X.T)) + self.b
return 1 if r > 0 else -1
def score(self, X_test, y_test):
r = 0
for i in range(len(X_test)):
x = X_test[i]
result = self.predict(x)
if result == y_test[i]:
r += 1
return r / len(X_test)
def weight(self):
yx = self.Y.reshape(-1, 1)*self.X
self.w = np.dot(yx.T, self.alpha)
return self.w
def main():
X, y = create_data()
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25)
plt.scatter(X[:50, 0], X[:50, 1], label='0')
plt.scatter(X[50:, 0], X[50:, 1], label='1')
svm = SVM((X_train, y_train))
svm.init_args()
svm.train()
score = svm.score(X_test, y_test)
print('score', score)
a1, a2 = svm.weight()
b = svm.b
x_min = min(svm.X, key=lambda x: x[0])[0]
x_max = max(svm.X, key=lambda x: x[0])[0]
y1, y2 = (-b - a1 * x_min)/a2, (-b - a1 * x_max)/a2
plt.plot([x_min, x_max], [y1, y2])
plt.show()
if __name__ == "__main__":
main()
0x05 结果
result
