一、简介
支持向量机（SVM）是经典的监督学习模型，这篇博文中，将通过线性规划（cvxopt包）展示SVM求解时的过程。数据集采用鸢尾花数据集进行二分类任务。
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from cvxopt import matrix, solvers
from sklearn.datasets import load_iris
iris = load_iris()
iris_df = pd.DataFrame(data= np.c_[iris["data"], iris["target"]], columns= iris["feature_names"] + ["target"])
# 二分类
iris_df = iris_df[iris_df["target"].isin([0,1])] ##找出0-1类
iris_df["target"] = iris_df[["target"]].replace(0,-1) ##0替换为-1
#选两个特征
iris_df = iris_df[["petal length (cm)", "petal width (cm)", "target"]]
##将数据转为numpy 可视化
X = iris_df[["petal length (cm)", "petal width (cm)"]].to_numpy()
y = iris_df[["target"]].to_numpy()
x_min = 0
x_max = 5.5
y_min = 0
y_max = 2
plt.figure(figsize=(8, 8))
colors = ["steelblue", "orange"]
plt.scatter(X[:, 0], X[:, 1], c=y.ravel(), alpha=0.5, cmap=matplotlib.colors.ListedColormap(colors), edgecolors="black")
plt.xlim(x_min, x_max)
plt.ylim(y_min, y_max)
plt.show()
数学来了
-
SVM最终的表达式化简如下
-
由拉格朗日定理，引入拉格朗日乘子α
- 因为J(w,b,α)是凸的，我们对其中的变量求微分如下
- 两个微分结果带入可得：
-
还是由微分结果，第二项为0了，最终结果如下：
-
我们用线性规划重新定义上试，定义一个矩阵H，Hij = yiyjxixj, 最终结果可改写成如下形式：
#取n个训练样本，定义矩阵H
n = X.shape[0]
H = np.dot(y*X, (y*X).T)
##第二项我们定于-1的列向量
q = np.repeat([-1.0], n)[..., None]
##第一个约束条件，我们定义A是1×n的矩阵，种类只有2类，令b=0
A = y.reshape(1, -1)
b = 0.0
##第二个约束条件，我们构造一个n×n单位矩阵G，和零向量h
G = np.negative(np.eye(n))
h = np.zeros(n)
##所有条件转为CVXOPT对象
P = matrix(H)
q = matrix(q)
G = matrix(G)
h = matrix(h)
A = matrix(A)
b = matrix(b)
##调用线性规划求解器
sol = solvers.qp(P, q, G, h, A, b)
alphas = np.array(sol["x"]) ##求出α
微分求解结果
我们定义支持向量,加入松弛因子,松弛因子的作用就是让约束条件不要太过绝对化。。。。
两个公式代码如下：
w = np.dot((y * alphas).T, X)[0]
S = (alphas > 1e-5).flatten()
b = np.mean(y[S] - np.dot(X[S], w.reshape(-1,1)))
##打印下结果
print("W:", w)
print("b:", b)
##可视化一下
xx = np.linspace(x_min, x_max)
a = -w[0]/w[1]
yy = a*xx - (b)/w[1]
margin = 1 / np.sqrt(np.sum(w**2))
yy_neg = yy - np.sqrt(1 + a**2) * margin
yy_pos = yy + np.sqrt(1 + a**2) * margin
plt.figure(figsize=(8, 8))
plt.plot(xx, yy, "b-")
plt.plot(xx, yy_neg, "m--")
plt.plot(xx, yy_pos, "m--")
colors = ["steelblue", "orange"]
plt.scatter(X[:, 0], X[:, 1], c=y.ravel(), alpha=0.5, cmap=matplotlib.colors.ListedColormap(colors), edgecolors="black")
plt.xlim(x_min, x_max)
plt.ylim(y_min, y_max)
plt.show()
第二方法调Sklearn包就很简单了
from sklearn import svm
clf = svm.SVC(kernel="linear", C=10.0)
clf.fit(X, y.ravel())
w = clf.coef_[0]
b = clf.intercept_
print("W:", w)
print("b:", b)
#W: [1.29411744 0.82352928]
#b: [-3.78823471]
完整代码如下：
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from cvxopt import matrix, solvers
from sklearn.datasets import load_iris
iris = load_iris()
iris_df = pd.DataFrame(data= np.c_[iris["data"], iris["target"]], columns= iris["feature_names"] + ["target"])
# 二分类
iris_df = iris_df[iris_df["target"].isin([0,1])] ##找出0-1类
iris_df["target"] = iris_df[["target"]].replace(0,-1) ##0替换为-1
#选两个特征
iris_df = iris_df[["petal length (cm)", "petal width (cm)", "target"]]
##将数据转为numpy 可视化
X = iris_df[["petal length (cm)", "petal width (cm)"]].to_numpy()
y = iris_df[["target"]].to_numpy()
x_min = 0
x_max = 5.5
y_min = 0
y_max = 2
plt.figure(figsize=(8, 8))
colors = ["steelblue", "orange"]
plt.scatter(X[:, 0], X[:, 1], c=y.ravel(), alpha=0.5, cmap=matplotlib.colors.ListedColormap(colors), edgecolors="black")
plt.xlim(x_min, x_max)
plt.ylim(y_min, y_max)
plt.show()
#取n个训练样本，定义矩阵H
n = X.shape[0]
H = np.dot(y*X, (y*X).T)
##第二项我们定于-1的列向量
q = np.repeat([-1.0], n)[..., None]
##第一个约束条件，我们定义A是1×n的矩阵，种类只有2类，令b=0
A = y.reshape(1, -1)
b = 0.0
##第二个约束条件，我们构造一个n×n单位矩阵G，和零向量h
G = np.negative(np.eye(n))
h = np.zeros(n)
##所有条件转为CVXOPT对象
P = matrix(H)
q = matrix(q)
G = matrix(G)
h = matrix(h)
A = matrix(A)
b = matrix(b)
##调用线性规划求解器
sol = solvers.qp(P, q, G, h, A, b)
alphas = np.array(sol["x"]) ##求出α
w = np.dot((y * alphas).T, X)[0]
S = (alphas > 1e-5).flatten()
b = np.mean(y[S] - np.dot(X[S], w.reshape(-1,1)))
##打印下结果
print("W:", w)
print("b:", b)
##可视化一下
xx = np.linspace(x_min, x_max)
a = -w[0]/w[1]
yy = a*xx - (b)/w[1]
margin = 1 / np.sqrt(np.sum(w**2))
yy_neg = yy - np.sqrt(1 + a**2) * margin
yy_pos = yy + np.sqrt(1 + a**2) * margin
plt.figure(figsize=(8, 8))
plt.plot(xx, yy, "b-")
plt.plot(xx, yy_neg, "m--")
plt.plot(xx, yy_pos, "m--")
colors = ["steelblue", "orange"]
plt.scatter(X[:, 0], X[:, 1], c=y.ravel(), alpha=0.5, cmap=matplotlib.colors.ListedColormap(colors), edgecolors="black")
plt.xlim(x_min, x_max)
plt.ylim(y_min, y_max)
plt.show()
##sklearn方法
from sklearn import svm
clf = svm.SVC(kernel="linear", C=10.0)
clf.fit(X, y.ravel())
w = clf.coef_[0]
b = clf.intercept_
print("W:", w)
print("b:", b)
#W: [1.29411744 0.82352928]
#b: [-3.78823471]
