论文复现：A Simple Decomposition Method for SupportVector Machines

 该版本代码未进行 working set 的选择，并使用了优化工具箱：

import numpy as np
from sklearn.metrics.pairwise import rbf_kernel
from sklearn import datasets
import matplotlib.pyplot as plt
import cvxpy as cp
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

from sklearn.svm import SVC

class BSVM():
    def __init__(self):
        self.X = self.y = None
        self.C = 1

    def fit(self, X, y):
        self.X = np.asarray(X, np.float)
        self.y = np.asarray(y, np.float)
        n, m = X.shape
        self.K = rbf_kernel(X=X,Y=X)
        self.Q = self.K * np.outer(y,y)
        self.H = self.Q + np.outer(y,y)
        self.a = np.zeros(n)


        G_1 = -np.eye(n)
        G_2 = np.eye(n)
        G = np.diag(np.vstack((G_1, G_2)))
        h_1 = np.zeros(n)
        h_2 = np.ones(n) * self.C
        h = np.hstack((h_1,h_2))


        self.alpha = cp.Variable(n)
        prob = cp.Problem(cp.Minimize((1/2) * cp.quad_form(self.alpha,self.H) - np.ones(n) @ self.alpha),
                          [G @ self.alpha <= h])

        prob.solve()
        print("alpha:::",self.alpha.value)

        for i in range(n):
            if self.alpha.value[i] <= 0:
                self.a[i] = 0
            elif self.alpha.value[i] >=1:
                self.a[i] = 1
            else:
                self.a[i] = self.alpha.value[i]
        print("a::",self.a)


    def predict(self,X):
        K_ = rbf_kernel(X=self.X,Y=X)
        # y_pred = self.alpha.value * self.y @ K_
        y_pred = self.a * self.y @ K_
        print("y_pred::",y_pred.shape)
        return np.sign(y_pred)

        #
        # print(y_pred)



if __name__ == '__main__':
    X,y = datasets.make_blobs(n_samples=500, n_features=2, centers=2, cluster_std=[4,4], random_state=11)
    print("原始标签：：",np.unique(y))
    y[y==0] = -1
    print("调整标签",np.unique(y))
    plt.scatter(X[:,0],X[:,1],c=y)
    plt.show()
    X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.3)


    model = BSVM()
    model.fit(X=X_train, y=y_train)
    y_pred = model.predict(X=X_test)

    acc = accuracy_score(y_true=y_test, y_pred=y_pred)

    model_svc = SVC(kernel='rbf')
    model_svc.fit(X=X_train, y=y_train)
    y_pred_1 = model_svc.predict(X=X_test)
    acc_1 = accuracy_score(y_true=y_test, y_pred=y_pred_1)


    print(acc)
    print(acc_1)

如下是使用scipy优化工具EQSQP 序列二次规划包，因此运行速度比较慢，本人认为使用有效集方法更加直接，速度可能会快一些。

import numpy as np
from sklearn.metrics.pairwise import rbf_kernel
from sklearn import datasets
import matplotlib.pyplot as plt
import cvxpy as cp
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

import scipy.optimize as so
from scipy.optimize import Bounds


from sklearn.svm import SVC
from time import time



class BSVM():
    def __init__(self):
        self.X = self.y = None
        self.C = 1

    def fit(self, X, y):
        self.X = np.asarray(X, np.float)
        self.y = np.asarray(y, np.float)
        self.n, self.m = X.shape
        n, m = X.shape
        self.K = rbf_kernel(X=X,Y=X)
        self.Q = self.K * np.outer(y,y)
        self.H = self.Q + np.outer(y,y)
        self.a = np.zeros(n)

        #
        # G_1 = -np.eye(n)
        # G_2 = np.eye(n)
        # G = np.diag(np.vstack((G_1, G_2)))
        # h_1 = np.zeros(n)
        # h_2 = np.ones(n) * self.C
        # h = np.hstack((h_1,h_2))

        bounds = Bounds(np.zeros(self.n), self.C * np.ones(self.n))
        alpha_0 = np.zeros(self.n)

        res = so.minimize(fun=self.func, x0=alpha_0, method="trust-constr", jac=self.Jac, hess=self.Hessian, options={'verbose':1}, bounds=bounds)
        self.a = res.x


        print("self.a:::",self.a)



    def func(self,alpha):
        '''The dual function'''
        return (1/2) * alpha.T @ self.H @ alpha - np.ones(self.n) @ alpha

    def Jac(self, alpha):
        '''The gradient '''
        return self.H @ alpha - np.ones(self.n)

    def Hessian(self, alpha):

        return self.H




    def predict(self,X):
        K_ = rbf_kernel(X=self.X,Y=X)
        # y_pred = self.alpha.value * self.y @ K_
        y_pred = self.a * self.y @ K_
        print("y_pred::",y_pred.shape)
        return np.sign(y_pred)

        #
        # print(y_pred)



if __name__ == '__main__':
    X,y = datasets.make_blobs(n_samples=500, n_features=2, centers=2, cluster_std=[4,4], random_state=13)
    print("原始标签：：",np.unique(y))
    X = X + 20.
    y[y==0] = -1
    print("调整标签",np.unique(y))
    plt.scatter(X[:,0],X[:,1],c=y)
    plt.show()
    X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.2)

    # ======================
    s_t = time()
    model = BSVM()
    model.fit(X=X_train, y=y_train)
    y_pred = model.predict(X=X_test)
    acc = accuracy_score(y_true=y_test, y_pred=y_pred)
    e_t = time()
    print("time:", e_t - s_t)


    # ======================
    s_t = time()
    model_svc = SVC(kernel='rbf',gamma='scale')
    model_svc.fit(X=X_train, y=y_train)
    y_pred_1 = model_svc.predict(X=X_test)
    acc_1 = accuracy_score(y_true=y_test, y_pred=y_pred_1)
    e_t = time()
    print("time:",e_t - s_t)


    print(acc)
    print(acc_1)

完整版

import cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics.pairwise import rbf_kernel
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import accuracy_score
from sklearn import datasets
from copy import deepcopy

from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
from time import time
from sklearn.svm import SVC


class Bsvm():
    def __init__(self):
        self.X = self.y = None

    def fit(self, X, y, C=1, q=10, tol=1E-3, epoch=50):
        self.X, self.y = np.asarray(X, np.float), np.asarray(y, np.float)
        self.tol = tol
        self.epoch = epoch
        self.C = C
        self.q = q
        self.n, self.m = X.shape
        self.Kernel = rbf_kernel(X=X, Y=X)
        self.alpha = np.zeros(self.n)
        self.idx_list = [_ for _ in range(self.n)]
        self.ws = None   # ====working set
        self.yy = np.outer(y,y)

        self.Q = self.Kernel * self.yy
        self.H = self.Q + self.yy


        for i in range(epoch):
            self.get_work_set()
            self.optimization()


    def get_work_set(self):
        work_set = []
        grad = self.H @ self.alpha - 1
        v_flag = np.zeros(self.n)
        v_value = np.zeros(self.n)
        for i in range(self.n):
            if self.alpha[i] == 0:
                v_value[i] = min(grad[i],0)
                if grad[i] >= 0:
                    v_flag[i] = 1
            elif self.alpha[i] == self.C:
                v_value[i] = -max(grad[i],0)
                if grad[i] <= 0:
                    v_flag[i] = 1
            elif 0 < self.alpha[i] < self.C:
                v_value[i] = - abs(grad[i])
                if grad[i] == 0:
                    v_flag[i] = 1

        r = sum(v_flag)
        ord_ids_list = np.argsort(v_value)
        if r > 0:
            first_num = min(int(self.q/2),r)
            second_num = self.q - first_num
            r_idx_list = np.where(v_flag==1)[0]
            for i in range(int(first_num)):
                work_set.append(r_idx_list[i])

            count = 0
            for i in range(self.n):
                if count >= second_num:
                    break
                if ord_ids_list[i] not in work_set:
                    work_set.append(ord_ids_list[i])
                    count += 1

        elif r == 0:
            half_q = int(self.q/2)
            positive_num = 0
            negative_num = 0
            for i in range(self.n):
                if positive_num >= half_q and negative_num >= half_q:
                    break
                if self.y[ord_ids_list[i]] == 1 and positive_num < half_q:
                    work_set.append(ord_ids_list[i])
                    positive_num += 1
                elif self.y[ord_ids_list[i]] == -1 and negative_num < half_q:
                    work_set.append(ord_ids_list[i])
                    negative_num += 1

        self.ws = work_set

    def optimization(self):
        self.ns = deepcopy(self.idx_list)
        for idx in self.ws:
            self.ns.remove(idx)

        Q_BB = self.Q[:,self.ws][self.ws,:]
        # print("Q_BB::",Q_BB)
        # print("----")
        # print(self.y[self.ws])
        # print(np.outer(self.y[self.ws],self.y[self.ws]))
        M_BB = Q_BB + np.outer(self.y[self.ws],self.y[self.ws])

        Q_BN = self.Q[:,self.ns][self.ws,:]
        print("Q_BN:",Q_BN.shape)
        print(np.outer(self.y[self.ws],self.y[self.ns]).shape)

        print(1-(Q_BN + np.outer(self.y[self.ws],self.y[self.ns])) @ self.alpha[self.ns])
        print("=======")
        M_BN = 1-(Q_BN + np.outer(self.y[self.ws],self.y[self.ns])) @ self.alpha[self.ns]


        alpha_ws = cp.Variable(len(self.ws))
        G1 = -np.eye(self.q)
        G2 = np.eye(self.q)

        G = np.diag(np.vstack((G1,G2)))
        print("G:",G.shape)
        print("alhpa_ws:",alpha_ws.shape)

        h1 = np.zeros(len(self.ws))
        h2 = np.ones(len(self.ws)) * self.C
        h = np.hstack((h1,h2))

        print("h:",h.shape)
        problem = cp.Problem(cp.Minimize((1/2) * cp.quad_form(x=alpha_ws,P=M_BB) - M_BN @ alpha_ws ),[G @ alpha_ws <= h])
        problem.solve()

        for i in range(len(self.ws)):
            if alpha_ws.value[i] <= self.tol:
                self.alpha[self.ws[i]] = 0

            elif alpha_ws.value[i] >= self.C - self.tol:
                self.alpha[self.ws[i]] = self.C

            else:
                self.alpha[self.ws[i]] = alpha_ws.value[i]

        #
        # print(alpha_ws.value.shape)
        # self.alpha[self.ws] = alpha_ws.value

    def predict(self, X):
        K_ = rbf_kernel(X=self.X, Y=X)
        # y_pred = self.alpha.value * self.y @ K_
        y_pred = self.alpha * self.y @ K_
        print("y_pred::", y_pred.shape)
        return np.sign(y_pred)





if __name__ == '__main__':
    X,y = datasets.make_blobs(n_samples=500, n_features=2, centers=2, cluster_std=[4,4], random_state=11)
    print("原始标签：：",np.unique(y))
    y[y==0] = -1
    print("调整标签",np.unique(y))
    plt.scatter(X[:,0],X[:,1],c=y)
    plt.show()
    X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.3)

    s_t = time()
    model = Bsvm()
    model.fit(X=X_train, y=y_train)
    y_pred = model.predict(X=X_test)
    acc = accuracy_score(y_true=y_test, y_pred=y_pred)
    e_t = time()

    print("time::",e_t - s_t)

    model_svc = SVC(kernel='rbf')
    model_svc.fit(X=X_train, y=y_train)
    y_pred_1 = model_svc.predict(X=X_test)
    acc_1 = accuracy_score(y_true=y_test, y_pred=y_pred_1)


    print(acc)
    print(acc_1)