Support Vector Machine实现（二）

soft margin的优化目标：    

                                          

最终转为：

                                                       

该等式优化问题使用SMO算法（sequential minimal optimization）进行训练

面向对象风格代码如下：

import numpy as np


class SupportVectorMachine:
    def __init__(self, maxIter=1e3, C=1, tol=1e-2, kernel=('lin', 0)):
        '''

        Parameters
        ----------
        maxIter：最大迭代次数默认1000次
        C：比例常数C
        tol：间隔容差
        kernel：核函数类型，默认为线性核
        '''
        self.X = None
        self.labels = None
        self.maxIter = maxIter
        self.C = C
        self.tol = tol
        self.kernel = kernel
        self.KernerMatrix = None   #样本数据输入核函数得到的Gram矩阵
        self.w = None   #线性分隔超平面的系数
        self.b = 0      #线性分隔超平面的截距
        self.alphas = None   #拉格朗日系数
        self.errCache = None

    def __initParams(self, X, labels):
        self.X = np.mat(X);
        self.labels = np.mat(labels.reshape((-1, 1)))
        m, n = X.shape
        self.alphas = np.mat(np.zeros((m, 1)))
        self.errCache = np.mat(np.zeros((m, 2)))
        self.w = np.zeros((n, 1))
        self.KernelMatrix = np.mat(np.zeros((m, m)))  # 计算Gram矩阵
        for i in range(m):
            self.KernelMatrix[:, i] = self.calKernelMatrix(X[i, :])

    def calKernelMatrix(self, sample):
        m, n = self.X.shape
        kernel = np.mat(np.zeros((m, 1)))
        if self.kernel[0] == 'lin':
            kernel = self.X * sample.T
        elif self.kernel[0] == 'rbf':
            for j in range(m):
                delta = self.X[j, :] - sample
                kernel[j] = delta * delta.T
            kernel = np.exp(kernel / (-2 * self.kernel[1] ** 2))
        else:
            raise NameError('not exists kernel function')
        return kernel

    def train(self, X, y):
        self.__initParams(X, y)
        self.__SequentialMinimalOptimization()

    def __calcEk(self, k):
        Fx_k = np.float(np.multiply(self.alphas, self.labels).T * self.KernelMatrix[:, k]) + self.b
        E_k = Fx_k - float(self.labels[k])
        return E_k

    def __updateE_k(self, k):
        E_k = self.__calcEk(k)
        self.errCache[k] = [1, E_k]

    @staticmethod
    def randSelectE_J(i, m):
        j = i
        while j == i:
            j = int(np.random.uniform(0, m))
        return j

    def __selectE_J(self, i, E_i):
        m = self.X.shape[0]
        maxK = -1
        maxDeltaE = 0
        E_j = 0
        self.errCache[i] = [1, E_i]
        validErrCacheList = np.nonzero(self.errCache[:, 0].A)[0]
        if len(validErrCacheList) > 1:
            for k in validErrCacheList:
                if k == i:
                    continue
                E_k = self.__calcEk(k)
                deltaE = abs(E_i - E_k)
                if deltaE > maxDeltaE:
                    maxK = k
                    maxDeltaE = deltaE
                    E_j = E_k
            return maxK, E_j
        else:
            j = self.randSelectE_J(i, m)
            E_j = self.__calcEk(j)
        return j, E_j

    @staticmethod
    def clipAlpha(aj, L, H):
        if aj > H:
            aj = H
        if aj < L:
            aj = L
        return aj

    def __calW(self):
        m, n = self.X.shape
        for i in range(m):
            self.w += np.multiply(self.alphas[i] * self.labels[i], self.X[i, :].T)

    def __innerL(self, i):
        E_i = self.__calcEk(i)
        if (self.labels[i] * E_i < -self.tol and self.alphas[i] < self.C) or (
                self.labels[i] * E_i > self.tol and self.alphas[i] > 0):
            j, E_j = self.__selectE_J(i, E_i)
            alpha_i_old = self.alphas[i].copy()
            alpha_j_old = self.alphas[j].copy()
            if self.labels[i] != self.labels[j]:
                L = max(0, self.alphas[j] - self.alphas[i])
                H = min(self.C, self.C + self.alphas[j] - self.alphas[i])
            else:
                L = max(0, self.alphas[j] + self.alphas[i] - self.C)
                H = min(self.C, self.alphas[j] + self.alphas[i])
            if L == H:
                print("L==H")
                return 0
            eta = self.KernelMatrix[i, i] + self.KernelMatrix[j, j] - 2 * self.KernelMatrix[i, j]
            if eta <= 0:
                print("eta<=0")
                return 0
            self.alphas[j] += self.labels[j] * (E_i - E_j) / eta
            self.alphas[j] = self.clipAlpha(self.alphas[j], L, H)
            self.__updateE_k(j)
            bi = self.b - E_i - self.labels[i] * (self.alphas[i] - alpha_i_old) * self.KernelMatrix[i, i] - \
                 self.labels[j] * (self.alphas[j] - alpha_j_old) * self.KernelMatrix[i, j]
            bj = self.b - E_j - self.labels[i] * (self.alphas[i] - alpha_j_old) * self.KernelMatrix[i, j] - \
                 self.labels[j] * (self.alphas[j] - alpha_j_old) * self.KernelMatrix[j, j]
            if 0 < self.alphas[i] < self.C:
                self.b = bi
            elif 0 < self.alphas[j] < self.C:
                self.b = bj
            else:
                self.b = (bi + bj) / 2
            return 1
        else:
            return 0

    def __SequentialMinimalOptimization(self):
        num = 0
        entireSet = True
        alphaPairsChanged = 0
        m = self.X.shape[0]
        while num < self.maxIter and (entireSet or alphaPairsChanged > 0):
            alphaPairsChanged = 0
            if entireSet:
                for i in range(m):
                    alphaPairsChanged += self.__innerL(i)
                    print("fullSet,iter: %d i:%d,pairs changed %d" % (num, i, alphaPairsChanged))
            else:
                nonBounds = np.nonzero((self.alphas.A > 0) * (self.alphas.A < self.C))[0]
                for i in nonBounds:
                    print("non-bound,iter:%d i:%d,pairs changed %d" % (num, i, alphaPairsChanged))
            num += 1
            if entireSet:
                entireSet = False
            elif alphaPairsChanged == 0:
                entireSet = True
            print("iteration number:%d" % num)
            print("entireSet:%s" % entireSet)
            print("alphaChanged:%d" % alphaPairsChanged)
        self.__calW()

    def getSupportVectors(self):
        supportVectors=self.X[self.alphas>0]
        return supportVectors

    def predict(self,X):
        y_pred=np.dot(X,self.w)
        y_pred[y_pred>0]=1
        y_pred[y_pred<=0]=-1
        return y_pred

测试代码如下：

from sklearn import datasets
from sklearn.model_selection import train_test_split

digits = datasets.load_digits()
X=digits.data
y=digits.target
X0=X[y==0]
X1=X[y==1]
y0=y[y==0]
y1=y[y==1]
X=np.concatenate((X0,X1))
y=np.concatenate((y0,y1))
X_train,X_test,y_train,y_test=train_test_split(X,y)


svm=SupportVectorMachine(maxIter=1e4)
svm.train(X,y)
y_pred=svm.predict(X)
print((y_pred==y).astype(np.int).mean())

参考文献：《机器学习实战》部分代码