svm c++实现_《统计学习方法》：SVM

今天终于实现了SVM，实现只用了一个小时，但是debug花了一整个晚上，下面我将列举一下实现中值得注意的点。实现过程十分艰辛，不过也乐趣非凡。另外，我还会将简要分析其他人的实现策略，以及存在的小bug等。

SVM

SVM作为本书中篇幅最长的一章，其基本思路是使得优化支撑向量之间的间隔最大化。

目标函数等价于最小化

，加入软间隔后，则在目标函数中加入了误分类的惩罚。

svm问题一般取其对偶形式进行优化，在对偶问题中，出现了原数据的内积

，可以通过核技巧，将内积转化为一个再生核希尔伯特空间上的内积

。从而引入了非线性，使得svm能适应更为复杂的分类情况。

核

由于核函数定义为内积，其具有对称的形式，实现也很简单，不过rbf核实现后发现分类效果极差，可能是数据集的问题，在此就只展示linear核以及poly核。

class Kernel:
    def __init__(self, name='linear', para=None):
        self.name = name
        self.para = para

    def __call__(self, data1, data2):
        if self.name == 'linear':
            return data1 @ data2.T
        if self.name == 'poly':
            if self.para == None:
                print('please assign a integer to p')
            return (data1 @ data2.T+1)**self.para

SMO算法

文中大部分时间都用在介绍SMO算法，这也是实现svm的重中之重。将svm转化为对偶问题后，优化变量变为了系数向量

, SMO算法的核心思想就是每次只优化两个变量，其他的保持不变，而由于svm问题中

，且对偶问题的约束条件保证这一子问题能转化为单变量优化问题，从而极大地提高了求解速度。

KKT条件

KKT 条件（7.111-113）中，

• 表示没有贡献，即处在支撑边界外的点
• 表示约束是积极的，即处在约束边界上的点
• 时，则点处在支撑边界内。

在我的实现中，由（7.105）有

，则

，故比较均和0进行对比，这一点和SmirkCao的实现思路类似。不同的是，他直接给出了违背KKT条件的判断方法。

def KKT(self, i):
    ye = self.y[i]*self.calE(i)
    ai = self.a[i]
    return np.isclose(ai, 0.0, atol=self.eps) and ye >= 0
        or self.eps < ai < self.C and np.isclose(ye, 0.0, atol=self.eps)
        or np.isclose(ai, self.C, atol=self.eps) and ye <= 0

这里必须指出，初始化优化变量

时应该从以0初始化，因为最后模型中大多数点都是没有贡献的，若使用随机初始化或者初始化为1，则会带来实践上的问题。比如wzyonggege的实现中是以1初始化的，优化完成后，所有的

均为1，这是实现中的一个错误，在此指出。 @yonggege

两个待优化变量的选取技巧

1. 筛选支撑向量集合，将其放置在其他点前进行循环，选取违背KKT的点1。
2. 首先在支撑向量集合中计算每个点的增益，代码中为delta，并选取增益最大的的点作为点2。
3. 若支撑向量集合中的点的增益均小于阈值，则继续在其他点中寻找。
4. 若其他点中也不能有足够的增益，则重新寻找点1。

def select_ij(self):
    active = [] # active set
    others = []
    for i in range(self.X.shape[0]):
        if 0 < self.a[i] < self.C:
            active.append(i)
        else:
            others.append(i)
    # 1. select i,j
    # active set first, if failed, try others
    for i in active+others:
        if not self.KKT(i):
            Ei = self.calE(i)
            okj, okEj, maxv, oka2 = 0, 0, 0, 0
            # active set first, if failed, try others
            for idx, j in enumerate(active+others):
                # 7.107
                eta = self.K[i, i]+self.K[j, j]-2*self.K[i, j]
                # page 126
                L, H = max(0, self.a[j]-self.a[i]), min(self.C, self.C+self.a[j]-self.a[i])
                if self.y[i] == self.y[j]:
                    L, H = max(0, self.a[j]+self.a[i]-self.C), min(self.C, self.a[j]+self.a[i])
                # odd case
                if L == H or eta <= 0:
                    continue
                Ej = self.calE(j)
                # 7.106
                a2 = self.a[j]+self.y[j]*(Ei-Ej)/eta
                a2 = np.clip(a2, L, H)
                # get the max update value in the active set
                delta = a2-self.a[j]
                if abs(delta) > maxv:
                    okj, okEj, maxv, oka2 = j, Ej, abs(delta), a2
                # if good step found in active, break
                if idx+1 == len(active) and not np.isclose(maxv, 0, self.eps):
                    break
            if np.isclose(maxv, 0, self.eps):
                continue
            return i, okj, Ei, okEj, oka2
    return None

另外，值得一提的是，SmirkCao的实现中，第一个变量是直接遍历1->n，第二个变量则采用随机的方法生成，同样具有较好的效果。（找不到对应知乎用户，暂且艾特一下一位提到了作者的 @红色石头 ）而wzyonggege的实现中没有针对第二个变量进行选取，直接按照书上的选取最大或最小的

，这样更新的点过少，实际操作中会使得优化难以进行下去。 @yonggege

更新变量

这个没什么好说的，直接从上一步优化变量后根据公式（7.115-116）更新偏置

即可。

def fit(self):
    it = 0
    while it < self.maxiter:
        res = self.select_ij()
        if res == None:
            break
        i, j, Ei, Ej, a2 = res
        a1 = self.a[i]+self.y[i]*self.y[j]*(self.a[j]-a2)
        # 3. cal b (7.115-116)
        b1 = self.b-Ei-self.y[i]*self.K[i, i] * (a1-self.a[i])
                      -self.y[j]*self.K[j, i]*(a2-self.a[j])
        b2 = self.b-Ej-self.y[i]*self.K[i, j] * (a1-self.a[i])
                      -self.y[j]*self.K[j, j]*(a2-self.a[j])
        # select b (under 7.116)
        self.b = b1 if 0 < a1 < self.C else (
                 b2 if 0 < a2 < self.C else (b1+b2)/2)
        self.a[i] = a1
        self.a[j] = a2
        it += 1

    if it == self.maxiter:
        print("Max iter number exceeded!")
    else:
        print("Done!")
    return it

其他部分的代码可以参见Github，希望看官老爷们能喜欢~

https://github.com/cherichy/statistical_learning/blob/master/python/svm.py​github.com