Python机器学习4-线性规划手撕SVM

一、简介

支持向量机（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()

数学来了

1. SVM最终的表达式化简如下

   
   

2. 由拉格朗日定理，引入拉格朗日乘子α

   
   
3. 因为J(w,b,α)是凸的，我们对其中的变量求微分如下

4. 两个微分结果带入可得：

5. 还是由微分结果，第二项为0了，最终结果如下：

   
   

6. 我们用线性规划重新定义上试，定义一个矩阵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]