孪生双子支持向量机（TWSVM）

线性可分支持向量机

孪生双子支持向量机与广义特征支持向量类似都是求取两个非平行的超平面来分离数据。（广义特征支持向量机https://blog.csdn.net/LIUGXIN/article/details/112132081）但是他们在本质上是不同的。TWSVM对中的两个二次规划问题都有一个典型的SVM公式。
通过求解以下公式得到TWSVM分类器二次规划问题：

这里的c1、c2 为常数且 c1、c2 > 0，e1、e2为合适维度的单位向量。
该算法找到两个超平面，每一类对应一个超平面，并根据该超平面对给定点最接近的点进行分类。上式目标函数的第一项为超平面到每一类数据距离的最小平方和。因此，最小化目标函数意味着，使一个超平面尽量离其中的一类数据近（class 1），约束要求超平面距离另一类数据(class -1)的距离至少大于1。当超平面距离小于最小距离1时，使用一组误差变量来测量误差。目标函数的第二项目标函数最小化误差变量的总和，试图最小化由于数据属于-1类而错分类的点。

TWSVM构建了两个更小的QPP问题，所以计算速度要比SVM更快，大约为：

针对TWSVM1构造拉格朗日函数：

α = (α1，α2…αm1)T，β = (β1，β2…βm2)T，解得TWSVM1的KKT条件为：

因为β>0,则 0< α <c1,前两个式子结合得到

使：

参数向量u = [w(1),b(1)]T，上式可以写为：

引入正则化项得到：

为了方便以下未添加正则化项，根据KKT条件，得到TWSVM1的对偶式如下：

TWSVM2与TWSVM1相同，对偶式如下：

P=[A e1]，Q=[B e2]，参数向量v = [w(2),b(2)]T，得到：

HT H和QTQ矩阵的大小是(n+1)x(n+1)，一般来说，n应该远远小于样本数量。
如何判断新的数据属于哪一类呢？
解上面矩阵，得到两个超平面

判断新的样本点距离那个超平面更接近，就将新的超平面分为相应的类。

| |代表的是点到平面的垂直距离。

非线性核分类

对于孪生双子支持向量机不是从线性扩展而来，而是直接从核出发构造优化问题。
补充：对于这个怎么得来的可能有些人不好理解，这儿可以看做是先对非线性数据进行核处理，也就是将非线性数据升到一个高维度中，在这个维度中，数据是线性可分的，然后在对此时的数据进行分类。

K为一个适合的核函数，可以选择线性核，高斯核等。
使用线性核：

引出优化问题：

得到KKT条件：

进一步得到：

使：

参数向量：

进一步简写为：

得到其对偶式：

相应的我们得到TWSVM2的优化问题如下：

同理可得TWSVM2的对偶式：

得到两个超平面后，对新样本点的预测方式与线性的一样。

代码

代码参考github https://github.com/vivamoto/classifier/tree/master/svm

import numpy as np
import scipy.optimize as optimize
import cvxopt
class tw_svm():
    def __init__(self, X_train, y_train, c1=1, c2=1, kernel_type='rbf', b=1, c=1, d=2, sigma=3):
        self.X_train = X_train
        self.y_train = y_train

        # Kernel 设置
        self.kernel_type = kernel_type
        self.b = b  # 双曲正切的常数乘子
        self.c = c  # 多项式，切线和线性样条的常数和
        self.d = d  # 多项式的权重
        self.sigma = sigma  # RBF 和 ERBF 的 sigma
        # Twin SVM settings
        self.c1 = c1  # TW-SVM 软间隔 1
        self.c2 = c2  # TW-SVM 软间隔 2
        # SVM settings
        self.C = 10                 # soft margin
        self.solver = 'optimize'   # Quadratic problem solver. 'cvxopt' or 'optimize'
        # #
        # self.alpha = None           # 拉格朗如乘子 (SVM and LSSVM)
        # self.SV = None              # 支持向量
        # self.bias = 1               # 偏置
    # 计算核函数 k(x1, x2)
    def kernel(self, x1, x2):

        if x1.ndim == 1:
            x1 = np.array([x1]).T

        if x2.ndim == 1:
            x2 = np.array([x2]).T

        m1 = x1.shape[0]
        m2 = x2.shape[0]
        k = np.zeros((m1, m2))

        t = self.kernel_type  # 选择使用的核函数
        b = self.b
        c = self.c
        d = self.d
        sigma = self.sigma
        for i in range(m1):
            for j in range(m2):
                # 线性核函数
                if t == 'linear':
                    k[i, j] = x1[i] @ x2[j]
                # 多项式核函数
                elif t == 'poly':
                    k[i, j] = (x1[i] @ x2[j] + c) ** d
                # 高斯核函数
                elif t == 'rbf':
                    k[i, j] = np.exp(-(x1[i] - x2[j]) @ (x1[i] - x2[j]) / (2 * sigma ** 2))
                # 指数径向基函数
                elif t == 'erbf':
                    k[i, j] = np.exp(-np.abs(x1[i] - x2[j]) / (2 * sigma ** 2))
                # H双曲正切核函数
                elif t == 'tanh':
                    k[i, j] = np.tanh(b * (x1[i] @ x2[j]) + c)
                # 线性样条核函数
                elif t == 'lspline':
                    k[i, j] = c + x1[i] * x2[j] + x1[i] * x2[j] * min(x1[i], x2[j]) + 1 / 2 * (x1[i] + x2[j]) * min(
                        x1[i], x2[j]) ** d
        return k


    # TWSVM: Solve quadratic problem
    def solveQP(p, q, x0, C=None):
        # "Quadratic Programming with Python and CVXOPT"
        # Solve SVM QP optimization problem:
        # min(x) 1/2 * x.T P x + q.T x
        # s.t. Gx<=h
        #      Ax=b
        # Input parameters of CVXOPT library
        N = q.shape[0]
        # construct P, q, G, h matrices for CVXOPT
        p = cvxopt.matrix(p, tc='d')
        q = cvxopt.matrix(q, tc='d')
        if C is None or C == 0:  # hard-margin SVM
            g = cvxopt.matrix(np.diag(np.ones(N) * -1), tc='d')
            h = cvxopt.matrix(np.zeros((N, 1)), tc='d')
        else:  # soft-margin SVM
            g = cvxopt.matrix(np.vstack((np.diag(np.ones(N) * -1), np.eye(N))), tc='d')
            h = cvxopt.matrix(np.hstack((np.zeros(N), np.ones(N) * C)), tc='d')
        # solve QP problem
        x = cvxopt.solvers.qp(p, q, g, h)

        return x

    # TWSVM: Solve quadratic problem using SciPy optimize
    def minimize(self, alpha, M):
        g = -np.sum(alpha) + 0.5 * alpha.T @ M @ alpha
        return g


    def train(self):
        X_train = self.X_train
        y_train = self.y_train
        if y_train.ndim == 1:
            y_train = np.array([y_train]).T
        A = X_train[y_train[:, 0] == 1]   # 取得正样本数据
        B = X_train[y_train[:, 0] == -1]  # 取得负样本数据
        self.Ct = np.vstack((A, B))       # 将样本拼接，用来计算核函数

        m1 = A.shape[0]
        m2 = B.shape[0]

        e1 = np.ones((m1, 1))
        e2 = np.ones((m2, 1))

        K_ACt = self.kernel(A, self.Ct)
        K_BCt = self.kernel(B, self.Ct)

        # TWSVM 1
        S = np.hstack((K_ACt, e1))
        R = np.hstack((K_BCt, e2))
        P1 = R @ np.linalg.pinv(S.T @ S) @ R.T

        # TWSVM 2
        L = np.hstack((K_ACt, e1))
        J = np.hstack((K_BCt, e2))
        P2 = L @ np.linalg.pinv(J.T @ J) @ L.T

        # 初始化 alpha 和 gamma
        alpha0 = np.zeros((np.size(R, 0), 1))
        gamma0 = np.zeros((np.size(L, 0), 1))

        # 拉格朗日乘数
        if self.solver == 'cvxopt':
            # CVXOPT algorithm:
            alpha = self.solveQP(P1, -e2, C=self.c1)
            gamma = self.solveQP(P2, -e1, C=self.c2)

            alpha = np.ravel(alpha['x'])
            gamma = np.ravel(gamma['x'])
        elif self.solver == 'optimize':
            # Scipy optimize
            b1 = optimize.Bounds(0, self.c1)
            b2 = optimize.Bounds(0, self.c2)
            alpha = optimize.minimize(self.minimize, x0=alpha0, args=P1, method='L-BFGS-B', bounds=b1).x
            gamma = optimize.minimize(self.minimize, x0=gamma0, args=P2, method='L-BFGS-B', bounds=b2).x

        if alpha.ndim == 1:
            alpha = np.array([alpha]).T

        if gamma.ndim == 1:
            gamma = np.array([gamma]).T

        # 解出参数
        epsilon = 1e-16
        I = np.eye(len(S.T @ S))
        self.z1 = -np.linalg.pinv(S.T @ S + epsilon * I) @ R.T @ alpha

        I = np.eye(len(J.T @ J))
        self.z2 = np.linalg.pinv(J.T @ J + epsilon * I) @ L.T @ gamma

        return

    def predict(self, X_test):
        # 输入参数:
        # z1, z2
        # Ct
        u1 = self.z1[:-1]   # TWSVM 1 的权重
        b1 = self.z1[-1:]   # TWSVM 1 的偏置
        u2 = self.z2[:-1]   # TWSVM 2 的权重
        b2 = self.z2[-1:]   # TWSVM 2 的偏置

        K_XCt = self.kernel(X_test, self.Ct)

        # 得到两个非平行的超平面
        surface1 = K_XCt @ u1 + b1
        surface2 = K_XCt @ u2 + b2

        # 计算点到两个超平面的距离
        dist1 = abs(surface1)  # class 1
        dist2 = abs(surface2)  # class -1

        # 根据数据点距离两个超平面的远近来判断所属的类别
        y_hat = np.argmax(np.hstack((dist1, dist2)), axis=1)
        if y_hat.ndim == 1:
            y_hat = np.array([y_hat]).T
        y_hat = np.where(y_hat == 0, -1, y_hat)
        return y_hat


def loadDataSet(fileName):
    dataMat = []
    labelMat = []
    fr = open(fileName)
    for line in fr.readlines():  # 逐行读取，滤除空格等
        lineArr = line.strip().split('\t')
        dataMat.append([float(lineArr[0]), float(lineArr[1])])  # 添加数据
        labelMat.append(float(lineArr[2]))  # 添加标签
    return np.array(dataMat), np.array(labelMat)


if __name__ == '__main__':
    X_train, y_train = loadDataSet('testSet1.txt')
    tw_svm = tw_svm(X_train, y_train)
    tw_svm.train()
    y_hat = tw_svm.predict(X_train)
    error = 0
    for i in range(len(y_train)):
        if y_train[i] != y_hat[i]:
            error += 1
    print(1-(error/len(y_train)))

数据https://download.csdn.net/download/LIUGXIN/14027036