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

作者: 西萌XXX | 来源:发表于2021-07-21 16:25 被阅读0次

    一、简介

    支持向量机(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,α)是凸的,我们对其中的变量求微分如下
    1. 两个微分结果带入可得:
    1. 还是由微分结果,第二项为0了,最终结果如下:


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

    相关文章

      网友评论

        本文标题:Python机器学习4-线性规划手撕SVM

        本文链接:https://www.haomeiwen.com/subject/emvlmltx.html