使用python的线性支持向量机中的软边距

Posted 2023-02-23

【中文标题】使用python的线性支持向量机中的软边距

【英文标题】：Soft margin in linear support vector machine using python

【发布时间】：2018-07-26 00:43:38

【问题描述】：

我正在学习支持向量机并试图提出一个简单的线性分类的简单 python 实现（我知道 sklearn 包，只是为了帮助更好地理解这些概念）。 This 是我引用的主要材料。

我试图通过最小化这个来解决原始的 SVM：

J wrt w 的导数是（根据上面的参考）：

所以这是使用“铰链”损失，而 C 是惩罚参数。如果我理解正确，设置更大的 C 将迫使 SVM 有更硬的余量。

下面是我的代码：

import numpy
from scipy import optimize

class SVM2C(object):
    def __init__(self,xdata,ydata,c=200.,learning_rate=0.01,
            n_iter=5000,method='GD'):

        self.values=numpy.unique(ydata)
        self.xdata=xdata
        self.ydata=numpy.where(ydata==self.values[-1],1,-1)
        self.c=c
        self.lr=learning_rate
        self.n_iter=n_iter
        self.method=method

        self.m=len(xdata)
        self.theta=numpy.random.random(xdata.shape[1])-0.5

    def costFunc(self,theta,x,y):
        zs=numpy.dot(x,theta)
        j=numpy.maximum(0.,1.-y*zs).mean()*self.c+0.5*numpy.sum(theta**2)
        return j

    def jac(self,theta,x,y):
        '''Derivative of cost function'''
        zs=numpy.dot(x,theta)
        ee=numpy.where(y*zs>=1.,0.,-y)[:,None]
        # multiply rows by ee
        dj=(ee*x).mean(axis=0)*self.c+theta
        return dj

    def train(self):

        #----------Optimize using scipy.optimize----------
        if self.method=='optimize':
            opt=optimize.minimize(self.costFunc,self.theta,args=(self.xdata,self.ydata),\
                    jac=self.jac,method='BFGS')
            self.theta=opt.x

        #---------Optimize using Gradient descent---------
        elif self.method=='GD':
            costs=[]
            lr=self.lr

            for ii in range(self.n_iter):
                dj=self.jac(self.theta,self.xdata,self.ydata)
                self.theta=self.theta-lr*dj
                cii=self.costFunc(self.theta,self.xdata,self.ydata)
                costs.append(cii)

            self.costs=numpy.array(costs)

        return self


    def predict(self,xdata):

        yhats=[]
        for ii in range(len(xdata)):
            xii=xdata[ii]
            yhatii=xii.dot(self.theta)
            yhatii=1 if yhatii>=0 else 0
            yhats.append(yhatii)
        yhats=numpy.array(yhats)

        return yhats



#-------------Main---------------------------------
if __name__=='__main__':

    xdata = numpy.array([[-1, -1], [-2, -1], [1, 1], [2, 1]])
    ydata = numpy.array([1, 1, 2, 2])

    mysvm=SVM2C(xdata,ydata,method='GD')
    mysvm.train()

    from sklearn import svm
    clf=svm.SVC(C=50,kernel='linear')
    clf.fit(xdata,ydata)

    print mysvm.theta
    print clf.coef_

    #-------------------Plot------------------------
    import matplotlib.pyplot as plt
    figure=plt.figure(figsize=(12,10),dpi=100)
    ax=figure.add_subplot(111)

    cmap=plt.cm.jet
    nclasses=numpy.unique(ydata).tolist()
    colors=[cmap(float(ii)/len(nclasses)) for ii in nclasses]

    #----------------Plot training data----------------
    for ii in range(len(ydata)):
        xii=xdata[ii][0]
        yii=xdata[ii][1]
        colorii=colors[nclasses.index(ydata[ii])]
        ax.plot(xii,yii,color=colorii,marker='o')

    plt.show(block=False)

所以这确实是一个玩具示例，其中只有 4 个线性可分训练样本，我已经删除了偏差项 b，结果 w 预期为 [0.5, 0.5]（skimage 结果），而我的实现无论是使用梯度下降还是scipy.optimize，都会倾向于给出大于 0.5 的值（例如 [1.4650, 1.4650]）。这只发生在将 C 参数设置为 >1 时，当 C==1 时，它给我 [0.5, 0.5]。同样当 C>1 时，scipy.optimize 会失败（我尝试了几种不同的方法，例如 Newton-CG、BFGS），尽管最终结果接近梯度下降结果。

我有点困惑为什么w 向量停止缩小。我认为当所有数据都被正确分类时，松弛惩罚将停止对总成本函数的贡献，因此它只会通过减小 w 的大小来优化 J。我是不是弄错了导数？

我知道这可能是一个新手问题，我正在粘贴一些脏代码，这让我困惑了几天，我周围没有人可以提供帮助，所以任何支持都将不胜感激！

更新：

感谢所有帮助。我正在更新代码以处理稍微复杂的示例。这次我包含了偏差项并使用以下内容对其进行了更新：

根据我得到的反馈，我为scipy.optimize 尝试了 Nelder-Mead，并尝试了 2 种自适应梯度下降方法。代码如下：

import numpy
from scipy import optimize

class SVM2C(object):
    def __init__(self,xdata,ydata,c=9000,learning_rate=0.001,
            n_iter=600,method='GD'):

        self.values=numpy.unique(ydata)
        # Add 1 dimension for bias
        self.xdata=numpy.hstack([xdata,numpy.ones([xdata.shape[0],1])])
        self.ydata=numpy.where(ydata==self.values[-1],1,-1)
        self.c=c
        self.lr=learning_rate
        self.n_iter=n_iter
        self.method=method

        self.m=len(xdata)
        self.theta=numpy.random.random(self.xdata.shape[1])-0.5

    def costFunc(self,theta,x,y):
        zs=numpy.dot(x,theta)
        j=numpy.maximum(0.,1.-y*zs).mean()*self.c+0.5*numpy.sum(theta[:-1]**2)
        return j

    def jac(self,theta,x,y):
        '''Derivative of cost function'''
        zs=numpy.dot(x,theta)
        ee=numpy.where(y*zs>=1.,0.,-y)[:,None]
        dj=numpy.zeros(self.theta.shape)
        dj[:-1]=(ee*x[:,:-1]).mean(axis=0)*self.c+theta[:-1] # weights
        dj[-1]=(ee*self.c).mean(axis=0)                      # bias

        return dj

    def train(self):

        #----------Optimize using scipy.optimize----------
        if self.method=='optimize':
            opt=optimize.minimize(self.costFunc,self.theta,args=(self.xdata,self.ydata),\
                    jac=self.jac,method='Nelder-Mead')
            self.theta=opt.x

        #---------Optimize using Gradient descent---------
        elif self.method=='GD':

            costs=[]
            lr=self.lr
            # ADAM parameteres
            beta1=0.9
            beta2=0.999
            epsilon=1e-8

            mt_1=0
            vt_1=0
            for ii in range(self.n_iter):
                t=ii+1
                dj=self.jac(self.theta,self.xdata,self.ydata)
                '''
                mt=beta1*mt_1+(1-beta1)*dj
                vt=beta2*vt_1+(1-beta2)*dj**2
                mt=mt/(1-beta1**t)
                vt=vt/(1-beta2**t)
                self.theta=self.theta-lr*mt/(numpy.sqrt(vt)+epsilon)
                mt_1=mt
                vt_1=vt

                cii=self.costFunc(self.theta,self.xdata,self.ydata)
                '''
                old_theta=self.theta
                self.theta=self.theta-lr*dj
                if ii>0 and cii>costs[-1]:
                    lr=lr*0.9
                    self.theta=old_theta


                costs.append(cii)
            self.costs=numpy.array(costs)

        self.b=self.theta[-1]
        self.theta=self.theta[:-1]

        return self


    def predict(self,xdata):

        yhats=[]
        for ii in range(len(xdata)):
            xii=xdata[ii]
            yhatii=numpy.sign(xii.dot(self.theta)+self.b)
            yhatii=xii.dot(self.theta)+self.b
            yhatii=self.values[-1] if yhatii>=0 else self.values[0]
            yhats.append(yhatii)
        yhats=numpy.array(yhats)

        return yhats

#-------------Main---------------------------------
if __name__=='__main__':

    #------------------Sample case 1------------------
    #xdata = numpy.array([[-1, -1], [-2, -1], [1, 1], [2, 1]])
    #ydata = numpy.array([1, 1, 2, 2])

    #------------------Sample case 2------------------
    from sklearn import datasets
    iris=datasets.load_iris()
    xdata=iris.data[20:,:2]
    ydata=numpy.where(iris.target[20:]>0,1,0)

    #----------------------Train----------------------
    mysvm=SVM2C(xdata,ydata,method='GD')
    mysvm.train()

    ntest=20
    xtest=2*(numpy.random.random([ntest,2])-0.5)+xdata.mean(axis=0)

    from sklearn import svm
    clf=svm.SVC(C=50,kernel='linear')
    clf.fit(xdata,ydata)

    yhats=mysvm.predict(xtest)
    yhats2=clf.predict(xtest)

    print 'mysvm weights:', mysvm.theta, 'intercept:', mysvm.b
    print 'sklearn weights:', clf.coef_, 'intercept:', clf.intercept_
    print 'mysvm predict:',yhats
    print 'sklearn predict:',yhats2

    #-------------------Plot------------------------
    import matplotlib.pyplot as plt
    figure=plt.figure(figsize=(12,10),dpi=100)
    ax=figure.add_subplot(111)

    cmap=plt.cm.jet
    nclasses=numpy.unique(ydata).tolist()
    colors=[cmap(float(ii)/len(nclasses)) for ii in nclasses]

    #----------------Plot training data----------------
    for ii in range(len(ydata)):
        xii=xdata[ii][0]
        yii=xdata[ii][1]
        colorii=colors[nclasses.index(ydata[ii])]
        ax.plot(xii,yii,color=colorii,marker='o',markersize=15)

    #------------------Plot test data------------------
    for ii in range(ntest):
        colorii=colors[nclasses.index(yhats2[ii])]
        ax.plot(xtest[ii][0],xtest[ii][1],color=colorii,marker='^',markersize=5)

    #--------------------Plot line--------------------
    x1=xdata[:,0].min()
    x2=xdata[:,0].max()

    y1=(-clf.intercept_-clf.coef_[0][0]*x1)/clf.coef_[0][1]
    y2=(-clf.intercept_-clf.coef_[0][0]*x2)/clf.coef_[0][1]

    y3=(-mysvm.b-mysvm.theta[0]*x1)/mysvm.theta[1]
    y4=(-mysvm.b-mysvm.theta[0]*x2)/mysvm.theta[1]

    ax.plot([x1,x2],[y1,y2],'-k',label='sklearn line')
    ax.plot([x1,x2],[y3,y4],':k',label='mysvm line')
    ax.legend(loc=0)
    plt.show(block=False)

我遇到的新问题：

它是不稳定的，取决于初始随机参数是什么，结果可能会有很大的不同。大约一半的时间，即使我将C 设置为一个相当大的值，它也会对训练集中的 1 个样本进行错误分类。 scipy.optimize 和 GD 都会出现这种情况。 ADAM 方法倾向于为vt 提供infs，至于大C，vt 增长非常快。我弄错渐变了吗？

提前致谢！

【问题讨论】：

我正在尝试遵循您的代码流，但是我很挣扎。 1) self.theta=aa.x 仅在 if 分支之一中调用 2) 第二个 aa=optimize(..) 在 if 语句之外调用。 3) 我看不到你调用 mysvm.predict()。希望对你有帮助

@MarcelFlygare，感谢您的指出。 1）这是 scipy.optimize 的正确调用。 2）这是我调试代码并已被删除。 3) predict() 没有被调用，因为我目前关心返回参数。预测值与 sklearn 结果相似，但并非始终相似。

我想你现在只需要稍微抑制一下。使用给定的参数，您的结果高度依赖于 theta 的初始随机状态。我在 (C=20, lr=0.01) 或 (C=200, lr=0.001) 得到了更好的数字

是的，它相当不稳定，这让我很担心，因为 svm 的成本函数据说是凸的，所以应该有一个全局最小值。我想大多数人都是从对偶形式来解决 svm 的，我仍在努力完全理解。

现在我意识到 [1.4650, 1.4650] 和 [0.5, 0.5] 定义了相同的平面/线，只是缩放不同。但这并没有纳入公式，是吗？而且我仍然不明白 C>1 如何使 scipy.optimize 失败而 C=1 不会。

【参考方案1】：

至于scipy.optimize，你误用了它的优化方法。 Newton-CG 和 BFGS 假设您的成本函数是平滑的，但事实并非如此。如果您使用稳健的 无梯度 方法，例如 Nelder-Mead，您将在大多数情况下收敛到正确的点（我已经尝试过）。

您的问题理论上可以通过梯度下降来解决，但前提是您将其调整为非平滑函数。目前，您的算法快速接近最优值，但开始跳跃而不是收敛，这是由于学习率大以及梯度的急剧变化，其中成本函数的最大值从 0 变为正值：

当您的成本相对于前一次迭代没有降低时，您可以通过降低学习率来平息这些振荡

def train(self):

    #----------Optimize using scipy.optimize----------
    if self.method=='optimize':
        opt=optimize.minimize(self.costFunc,self.theta,args=(self.xdata,self.ydata),\
                jac=self.jac,method='BFGS')
        self.theta=opt.x

    #---------Optimize using Gradient descent---------
    elif self.method=='GD':
        costs=[]
        lr=self.lr

        for ii in range(self.n_iter):
            dj=self.jac(self.theta,self.xdata,self.ydata)
            old_theta = self.theta.copy()
            self.theta=self.theta-lr*dj
            cii=self.costFunc(self.theta,self.xdata,self.ydata)

            # if cost goes up, decrease learning rate and restore theta
            if len(costs) > 0 and cii > costs[-1]:
                lr *= 0.9
                self.theta = old_theta
            costs.append(cii)

        self.costs=numpy.array(costs)

    return self

对代码的这个小修改会带来更好的收敛性：

得到的参数非常接近最优值 - 例如 [0.50110433 0.50076661] 或 [0.50092616 0.5007394 ]。

在现代应用程序（如神经网络）中，这种学习率适应是在 ADAM 等高级梯度下降算法中实现的，该算法不断跟踪梯度均值和方差的变化。

更新。答案的第二部分涉及代码的第二个版本。

关于亚当。因为vt=vt/(1-beta2**t) 这条线，你的vt 爆炸了。您应该只对用于计算梯度步长的 vt 值进行归一化，而不是对进入下一次迭代的值进行归一化。喜欢这里：

...
mt=beta1*mt_1+(1-beta1)*dj
vt=beta2*vt_1+(1-beta2)*dj**2
mt_temp=mt/(1-beta1**t)
vt_temp=vt/(1-beta2**t)
old_theta=self.theta
self.theta=self.theta-lr*mt_temp/(numpy.sqrt(vt_temp)+epsilon)
mt_1=mt
vt_1=vt
...

关于不稳定。 Nelder-Mead 方法和梯度下降都依赖于参数的初始值，这是可悲的事实。您可以尝试通过以更明智的方式进行更多的 GD 迭代和衰减学习率，或者通过减少 Nelder-Mead 方法的 xatol 和 fatol 等优化参数来提高收敛性。

但是，即使您实现了完美的收敛（在您的情况下为 [ 1.81818459 -1.81817712 -4.09093887] 之类的参数值），您也会遇到问题。收敛性可以通过以下代码粗略检查：

print(mysvm.costFunc(numpy.concatenate([mysvm.theta, [mysvm.b]]), mysvm.xdata, mysvm.ydata))
print(mysvm.costFunc(numpy.concatenate([mysvm.theta, [mysvm.b+1e-3]]), mysvm.xdata, mysvm.ydata))
print(mysvm.costFunc(numpy.concatenate([mysvm.theta, [mysvm.b-1e-3]]), mysvm.xdata, mysvm.ydata))
print(mysvm.costFunc(numpy.concatenate([mysvm.theta-1e-3, [mysvm.b]]), mysvm.xdata, mysvm.ydata))
print(mysvm.costFunc(numpy.concatenate([mysvm.theta+1e-3, [mysvm.b]]), mysvm.xdata, mysvm.ydata))

导致

6.7323592305075515
6.7335116664996
6.733895813394582
6.745819882839341
6.741974212439457

如果您更改 theta 或任何方向的截距，您的成本会增加 - 因此，解决方案是最佳的。但是sklearn 的解决方案不是最优的（从mysvm 的角度来看），因为代码

print(mysvm.costFunc(numpy.concatenate([clf.coef_[0], clf.intercept_]), mysvm.xdata, mysvm.ydata))

打印40.31527145374271！这意味着，您已经达到了局部最小值，但 sklearn 的 SVM 已经最小化了一些不同的东西。

如果你是read the documentation 或sklearn，你会发现问题所在：他们最小化sum(errors) * C + 0.5 * penalty，而你最小化mean(errors) * C + 0.5 * penalty！！！这是造成差异的最可能原因。

【讨论】：

漂亮的答案！

感谢您的见解！您建议的答案有所帮助，但是对于稍微复杂的情况（也是线性可分的），我仍然遇到解决原始问题的问题，这次是梯度中包含了偏差项。我也尝试过 ADAM，但它倾向于在 vt 术语中给出 infs。而且无论是使用梯度下降还是 Nelder-Meader，无论C 有多大，它总是错误分类 1 个样本。我开始觉得双重形式更强大。我应该更新问题以包含较新的案例吗？

@Jason，因为vt=vt/(1-beta2**t) 这行，你的vt 爆炸了。您应该只对用于计算梯度步长的 vt 的值进行归一化，而不是对进入下一次迭代的值进行归一化。

@DavidDale 确实，谢谢！现在它没有爆炸，但收敛速度有点慢（大约 8k 次迭代后），并且仍然存在单个错误分类。

@Jason 收敛速度可以通过学习率来固定。收敛后，您可以检查您的解决方案是否确实是局部最优的（成本低于邻点），如果是这种情况，那么您可能在成本函数中犯了一种愚蠢的错误。