使用Python的线性支持向量机中的软间隔

12
我正在学习支持向量机,并尝试编写一个简单的Python实现(我知道sklearn包,只是为了更好地理解概念),用于进行简单的线性分类。是我参考的主要材料。
我正在尝试从原始问题开始解决SVM,通过最小化以下内容:

enter image description here

J对w的导数是(根据上面的参考文献):

enter image description here

所以这里使用了“铰链”损失函数,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结果),而我的实现将倾向于给出大于0.5的值(例如[1.4650,1.4650]),无论是使用梯度下降还是scipy.optimize。只有当C>1时,才会发生这种情况,当C==1时,它会给我[0.5,0.5]。此外,当C>1时,scipy.optimize会失败(我尝试了几种不同的方法,例如Newton-CG、BFGS),尽管最终结果接近于梯度下降结果。
我有点困惑为什么w向量停止缩小。我认为当所有数据都被正确分类时,松弛惩罚将停止对总成本函数的贡献,因此它只会通过减少w的大小来优化J。我算错了吗?
我知道这可能是一个新手问题,我在粘贴一些脏代码,这已经困扰我几天了,我周围没有人可以提供帮助,所以任何支持都将不胜感激!
更新:
感谢所有的帮助。我正在更新代码以处理稍微复杂一些的样本。这次我包括了偏置项,并使用以下内容进行更新:

enter image description here

根据我收到的反馈,我尝试了scipy.optimize的Nelder-Mead方法,并尝试了两种自适应梯度下降方法。以下是代码:

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方法倾向于给出vtinf值,因为对于较大的Cvt增长得非常快。我是不是在计算梯度时出了问题?

提前无限感谢!


我正在尝试跟随你的代码流程,但是我很困惑。1)self.theta=aa.x只在其中一个if分支中调用。2)第二个aa=optimize(..)在if语句之外被调用。3)我看不到你调用mysvm.predict()。 希望这有所帮助。 - Marcel Flygare
1
现在我意识到[1.4650,1.4650]和[0.5,0.5]定义了相同的平面/线,只是缩放不同。但这没有纳入公式中,是吗?我仍然不明白为什么C>1会导致scipy.optimize失败而C=1不会。 - Jason
你可以通过打印“opt”获取错误,并通过设置容差,例如tol=0.1来消除错误。(由于您比我更了解问题,因此我没有进行进一步的深入分析。) - Marcel Flygare
关于基本的优化理论:您正在使两种方法的假设无效:平滑性(铰链损失)!此外,GD的收敛需要一些注意:例如线搜索。 - sascha
仔细查看参考资料表明,在使用次梯度下降时,参数w应该每个样本更新一次,而我却将所有样本的导数平均值取了,因此结果可能是错误的。但仅仅改变为每个样本的方案并不能解决所有问题,还有更多我不理解的地方。 - Jason
显示剩余3条评论
1个回答

11
关于scipy.optimize,您误用了其优化方法。 Newton-CG和BFGS都假定您的成本函数是平滑的,但事实并非如此。 如果您使用强大的无梯度方法,例如Nelder-Mead,则在大多数情况下会收敛到正确的点(我已经尝试过)。
理论上,您的问题可以通过梯度下降来解决,但只有在将其适应为非平滑函数时才能实现。 目前,您的算法快速接近最优解,但由于学习率较大且成本函数中的最大值从0变为正数的梯度发生了急剧变化,导致开始“跳跃”而不是收敛。

enter image description here

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

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

这个小修改会使你的代码收敛得更好:

enter image description here

在结果参数中,这些参数非常接近最优解,例如[0.50110433 0.50076661][0.50092616 0.5007394 ]
在现代应用程序(如神经网络)中,学习率的适应是通过高级梯度下降算法(如ADAM)实现的,该算法不断跟踪梯度的均值和方差。
更新。此答案的第二部分涉及代码的第二个版本。
关于ADAM。您得到了爆炸性的vt,因为有一行代码vt=vt/(1-beta2**t)。您应该仅对用于计算梯度步长的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方法的优化参数如xatolfatol来改善收敛性。
然而,即使您实现了完美的收敛(例如[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或截距,则成本会增加-因此,解决方案是最佳的。但是,从mysvm的角度来看,sklearn的解决方案不是最优的,因为代码{{code}}。

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

打印40.31527145374271!这意味着你已经到达了一个局部最小值,但是sklearn的SVM最小化了不同的东西。

如果你阅读sklearn的文档,你可以找到问题所在:他们最小化了sum(errors) * C + 0.5 * penalty,而你最小化了mean(errors) * C + 0.5 * penalty!这是差异的最可能原因。


漂亮的答案! - MaxU - stand with Ukraine
谢谢您的见解!您提供的答案有所帮助,但我仍然在解决一个稍微复杂一些的原始问题时遇到了问题(也是线性可分的),这次包括偏置项在梯度中。我也尝试过ADAM,但它往往会在“vt”项中给出“inf”。无论是使用梯度下降还是Nelder-Meader,它总是错误地分类1个样本,无论“C”有多大。我开始觉得对偶形式更加稳健。我应该更新问题以包括新的情况吗? - Jason
@Jason,你的vt出现了爆炸性的问题,是因为这一行代码vt=vt/(1-beta2**t)。你应该只对用于计算梯度步长的vt值进行归一化,而不是下一次迭代的值。 - David Dale
@DavidDale 确实,谢谢!现在它不会崩溃了,但是收敛速度有点慢(大约8k次迭代后),而且仍然存在单个错误分类。 - Jason
@Jason 收敛速度可能通过调整学习率来固定。在收敛之后,您可以检查您的解决方案是否确实是一个局部最优解(成本比邻域点低),如果是这种情况,那么您可能在成本函数中犯了某种愚蠢的错误。 - David Dale
显示剩余3条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接