简介
Yuan-Hai Shao等人提出了一种改进的TWSVM,也称为孪生有界支持向量机(Twin bound Support vector Machine, TBSVM)。使用逐次超松驰(SOR)技术,针对优化问题,以提高训练过程的速度,还应用了结构风险最小化原则,在计算时间和分类精度方面,TBSVM比TWSVM更优秀。建议先读风险最小化(https://blog.csdn.net/LIUGXIN/article/details/112094477)
线性TBSVM
与TWSVM一样,TBSVM也是寻找两个非平行的决策超平面来判别数据。
线性TBSVM的原问题如下:
c1,c2,c3,c4是惩罚参数,目标函数第一项是最小化训练误差的额外正规化项。这里与TWSVM做个对比,显而易见,TBSVM多了一个正则化项。由于加入的这一项,使得结构风险最小化了。
在普通的SVM中,结构风险最小化是通过最大化两类之间的间隔来实现的,该间隔是通过两个支持超平面之间的欧氏距离来度量的。相应的,对于TWSVM两个类之间的间隔可以通过近端超平面wTx+b = 0和边界超平面wTx+b = -1之间的某种距离来测量。下面是一种可行性的距离(两直线之间的距离公式)。
为了求解原问题,我们得到其拉格朗日对偶式:
α,β为拉格朗日乘子向量,解得KKT条件为:
整合前两个式子得:
I为适当维度的单位矩阵,定义v1= [w1,b1]T,上式可以写为:
接下来我们得到了,这个问题的对偶问题:
这个式子乍一看与TWSVM的对偶式只是换了个c3,但是他们在含义上有着十分重要的差异。TWSVM中仅仅是一个不变的标量,在TBSVM中C3是一个权重因子,决定了正规化项和经验风险之间的权衡。因此,选择合适的c3,无论大小,都体现了结构风险最小化的原则。
同样的方式我们得到另一个对偶式:
解出两个超平面后,通过下式判断数据的分类。
非线性TBSVM
与TWSVM类似,得到如下超平面
CT= [A B]T,k为合适的核函数,对第一个平面构建原始优化问题:
c1,c3>0,均为标量,可得以下拉格朗日函数:
α = (α1,α2…αm1)T,β = (β1,β2…βm2)T,为拉格朗日乘子向量。得到KKT条件如下:
β>0=0 ,得到 0 =< α <= c,结合前两式,显然得到:
定义:
得到参数矩阵 z1 = [u1, b1]T:
得到对偶问题如下:
同理可得另一超平面的原始问题:
对偶式为:
参数矩阵为:
解出两平面参数矩阵的值,得到以下分类公式:
即将数据点分到距离更近的平面。
代码
import numpy as np
import scipy.optimize as optimize
import cvxopt
class tb_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
# 附加项参数
self.c3 = 1
self.c4 = 1
# 二次规划问题计算器 'cvxopt' or 'optimize'
self.solver = 'optimize'
# Kernel 设置
self.kernel_type = kernel_type
self.b = b # 双曲正切的常数乘子
self.c = c # 多项式,切线和线性样条的常数和
self.d = d # 多项式的权重
self.sigma = sigma # RBF 和 ERBF 的 sigma
# 计算核函数 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(self, p, q, C=None):
# 解决 SVM QP 优化问题:
# 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: 利用SciPy优化解决二次型问题
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] # 取得负样本数据
# 拼接AB为求k做准备
self.Ct = np.vstack((A, B))
m1 = A.shape[0]
m2 = B.shape[0]
e1 = np.ones((m1, 1))
e2 = np.ones((m2, 1))
I = np.eye(m1 + m2 + 1)
# 求AB对于核函数的值
K_ACt = self.kernel(A, self.Ct)
K_BCt = self.kernel(B, self.Ct)
# TBSVM 1
S = np.hstack((K_ACt, e1))
R = np.hstack((K_BCt, e2))
P1 = R @ np.linalg.pinv(S.T @ S + self.c3 * I) @ R.T
# TBSVM 2
L = np.hstack((K_ACt, e1))
J = np.hstack((K_BCt, e2))
P2 = L @ np.linalg.pinv(J.T @ J + self.c4 * I) @ 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 算法:
alpha = self.solveQP(P1, -e2, C=self.c3)
gamma = self.solveQP(P2, -e1, C=self.c4)
alpha = np.ravel(alpha['x'])
gamma = np.ravel(gamma['x'])
elif self.solver == 'optimize':
# Scipy optimize
b1 = optimize.Bounds(0, self.c3)
b2 = optimize.Bounds(0, self.c4)
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')
tb_svm = tb_svm(X_train, y_train)
tb_svm.train()
X_test, y_test = loadDataSet('testSet.txt')
y_hat = tb_svm.predict(X_test)
error = 0
for i in range(len(y_test)):
if y_test[i] != y_hat[i]:
error += 1
print(1-(error/len(y_test)))