最优化方法Python计算：有约束优化应用——近似线性可分问题支持向量机

二分问题的数据集$\{({\mathbf{x}}_i,y_i)\}$，$i=1,2,\cdots ,m$中，特征数据 { x i } \{\boldsymbol{x}_i\} {xi​}未必能被一块超平面按其标签值 y i ∈ { − 1 , 1 } y_i\in\{-1,1\} yi​∈{−1,1}分隔，即二分问题未必是线性可分的。如果存在超平面$\pi$：$({\mathbf{x}}^{\top },1){\mathbf{w}}_0=0$，使得大多数样本点的特征数据被$\pi$按标签数据取值分隔，满足约束$y_i({\mathbf{x}}_i,1){\mathbf{w}}_0\ge 1$，但有少数样本点特征数据不满足此约束而交叉纠结，如下图所示。这样的问题，称为近似线性可分问题。图中的超平面$({\mathbf{x}}^{\top },1){\mathbf{w}}_0=\pm 1$称为软边界。

解决近似线性可分问题的支持向量机模型引入松弛变量$\mathbf{z}=\left(\begin{array}{c}{z}_1\\ z_2\\ ⋮\\ z_m\end{array}\right)$，构造优化问题
$$\{\begin{array}{l}\min \frac{1}{2}{\mathbf{w}}^{\top }\mathbf{Hw}+{\mathbf{c}}^{\top }\mathbf{z}\\ \text{s.t. }{y}_i({\mathbf{x}}_i,1)\mathbf{w}\ge 1-z_i\\ z_i\ge 0\end{array}i=1,2,\cdots m.(1)$$
其中，$\mathbf{H}=\left(\begin{array}{ccccc}1& 0& \cdots & 0& 0\\ 0& 1& \cdots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & 1& 0\\ 0& 0& \cdots & 0& 0\end{array}\right)$，$\mathbf{c}=\left(\begin{array}{c}C\\ C\\ ⋮\\ C\end{array}\right)$，正则化参数$C>0$为松弛变量的罚因子。其值越大，模型对误分类的惩罚越强；反之，模型将侧重于分类间隔。当C足够大时，上述问题退化为问题。
$$\{\begin{array}{l}\min \frac{1}{2}{\mathbf{w}}^{\top }\mathbf{Hw}\\ \text{s.t.}{y}_i({\mathbf{x}}_i^{\top },1)\mathbf{w}\ge 1i=1,2,\cdots ,m.\end{array}$$
根据强对偶定理，可以证明二次规划问题(1)的最优解${\mathbf{w}}_0$与其对偶问题
$$\{\begin{array}{l}\min \frac{1}{2}{\mathbf{\mu }}^{\top }\mathbf{Q\mu }-{\mathbf{\mu }}^{\top }1\\ \text{s.t. }{\mathbf{\mu }}^{\top }\mathbf{y}=0\\ \mathbf{o}\le \mathbf{\mu }\le \mathbf{c}\end{array}.(2)$$
的最优解${\mathbf{\mu }}_0$等价。其中，$\mathbf{Q}={\mathbf{X}}^{\top }\mathbf{X}$，$\mathbf{X}=(y_1{\mathbf{x}}_1,y_2{\mathbf{x}}_2,\cdots ,y_m{\mathbf{x}}_m)$，$\mathbf{\mu }\in {\text{R}}^m$。等价的意义为
$${\mathbf{w}}_0=\left(\begin{array}{c}\sum _{i\in s}{\mu }_{0_i}{y}_i{\mathbf{x}}_i\\ \frac{1}{m_s}\sum _{j\in s}\left(y_j-\sum _{i\in s}{\mu }_{0_i}{y}_i{\mathbf{x}}_j^{\top }{\mathbf{x}}_i\right)\end{array}\right).$$
其中，
s = { i ∣ 0 ≤ i ≤ m , 0 < μ 0 i < C } s=\{i|0\leq i\leq m, 0<\mu_{0_i}<C\} s={i∣0≤i≤m,0<μ0i​​<C}
为支持向量下标集，$m_s$为$s$支持向量个数。
下列代码实现通过求解二次规划问题(2)的近似线性可分问题支持向量机模型。

import numpy as np
from scipy.optimize import minimize
class ALineSvc(Classification, LineModel):								#近似线性可分支持向量机分类器def __init__(self, C = 1e+4):										#构造函数self.C = Cself.tagVal = np.signdef obj(self, mu):													#目标函数return 0.5 * (mu @ (self.Q @ mu)) - mu.sum()def ynormalize(self, y, trained):									#标签数据预处理if not trained:self.ymin = 0self.ymax = 1return (y - self.ymin) / (self.ymax - self.ymin)def fit(self, X, Y, mu = None):										#训练函数print("训练中...，稍候")m, n = X.shapeself.scalar = (len(X.shape) == 1)self.X, self.Y = self.pretreat(X, Y)if not isinstance(mu, np.ndarray):if mu == None:mu = np.random.random(m)else:mu = np.array([mu] * m)Qk = np.array([[self.X[i] @ self.X[j]							#内积矩阵for j in range(m)]for i in range(m)])self.Q=np.outer(self.Y, self.Y) * Qk							#系数矩阵h = lambda x: self.Y @ x										#等式约束条件函数g1 = lambda x: x												#不等式约束函数1g2 = lambda x: self.C - x										#不等式约束函数2cons = [{'type': 'eq', 'fun': h},								#约束条件列表{'type': 'ineq', 'fun': g1},{'type': 'ineq', 'fun': g2}]res = minimize(self.obj, mu, constraints = cons)				#求解最优化问题self.mu0 = res.xself.support_ = np.where((self.mu0 > 1e-5) &					#支持向量下标(self.mu0 < self.C))[0]self.w0 = (self.mu0[self.support_] * self.Y[self.support_]) @\self.X[self.support_,0:n]									#模型参数前n项b0 = (self.Y[self.support_]-(self.mu0[self.support_] * self.Y[self.support_])\@ Qk[np.ix_(self.support_, self.support_)]).mean()	#模型参数最末项self.w0 = np.append(self.w0, b0)								#整合模型参数self.coef_, self.intercept_ = self.coef_inte()					#计算超平面系数和截距print("%d次迭代后完成训练。"%res.nit)

程序中第3~44行定义的近似线性可分支持向量机分类器类ALineSvc继承了Classification（见博文《最优化方法Python计算：无约束优化应用——线性回归分类器》）和LineModel（见博文《最优化方法Python计算：无约束优化应用——线性回归模型》）的属性与方法。类定义体中

• 第4~6行定义构造函数，设置松弛变量上限C的默认值为$10^4$。第6行将标签值函数tagVal设置为Numpy的sign函数，以适于计算分类预测值。
• 第7~8行定义问题(2)的目标函数obj，返回值为$\frac{1}{2}{\mathbf{\mu }}^{\top }\mathbf{Q}\mathbf{\mu }-{\mathbf{\mu }}^{\top }1$。
• 第9~11行重载标签数据归一化函数ynormalize。由于近似线性可分支持向量机模型中的标签数据不需要归一化，所以在训练时将self.ymin和self.ymax设置为0和1。第12行返回归一化后的标签数据。做了这样的调整，我们在进行预测操作时，按式
  $$y=y\cdot (\max y-\min y)+\min y$$
  仍保持$y$的值不变，进而可保持LineModel的predict函数代码不变（见程序4.1及其说明）。
• 第14~44行重载训练函数fit。对比程序4.1中定义的父类fit函数可见第15~23行的代码是保持一致的。
  - 第24~26行计算内积矩阵$${\mathbf{Q}}_k=\left(\begin{array}{cccc}{\mathbf{x}}_1^{\top }{\mathbf{x}}_1& {\mathbf{x}}_1^{\top }{\mathbf{x}}_2& \cdots & {\mathbf{x}}_1^{\top }{\mathbf{x}}_m\\ {\mathbf{x}}_2^{\top }{\mathbf{x}}_1& {\mathbf{x}}_2^{\top }{\mathbf{x}}_2& \cdots & {\mathbf{x}}_2^{\top }{\mathbf{x}}_m\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{x}}_m^{\top }{\mathbf{x}}_1& {\mathbf{x}}_m^{\top }{\mathbf{x}}_2& \cdots & {\mathbf{x}}_m^{\top }{\mathbf{x}}_m\end{array}\right)$$
  第27行计算目标函数系数矩阵$$\mathbf{Q}=\left(\begin{array}{cccc}{y}_1{y}_1{\mathbf{x}}_1^{\top }{\mathbf{x}}_1& y_1{y}_2{\mathbf{x}}_1^{\top }{\mathbf{x}}_2& \cdots & y_1{y}_m{\mathbf{x}}_1^{\top }{\mathbf{x}}_m\\ y_2{y}_1{\mathbf{x}}_2^{\top }{\mathbf{x}}_1& y_2{y}_2{\mathbf{x}}_2^{\top }{\mathbf{x}}_2& \cdots & y_2{y}_m{\mathbf{x}}_2^{\top }{\mathbf{x}}_m\\ ⋮& ⋮& \ddots & ⋮\\ y_m{y}_1{\mathbf{x}}_m^{\top }{\mathbf{x}}_1& y_m{y}_2{\mathbf{x}}_m^{\top }{\mathbf{x}}_2& \cdots & y_m{y}_m{\mathbf{x}}_m^{\top }{\mathbf{x}}_m\end{array}\right).$$
  - 第28~30行定义等式约束条件函数h和不等式约束条件函数g1和g2。第31~33行构造约束条件列表cons。第34行调用minimize函数求解优化问题(2)，返回值赋予res。第35行将res的x属性赋予mu0，为问题(2)的最优解${\mathbf{\mu }}_0$。
  - 第36~37行按条件
  $$0<{\mu }_{0_i}<C,i=1,2,\cdots ,m$$
  查找支持向量的下标集$s$并赋予属性support_。
  - 第38~39行计算原问题(1)的最优解${\mathbf{w}}_0$的前$n$个元素$\sum _{i\in s}{\mu }_{0_i}{y}_i{\mathbf{x}}_i$。第40~41行计算原问题(1)的最优解${\mathbf{w}}_0$的最后一个元素$b_0$。第43行连接两部分构成${\mathbf{w}}_0$赋予属性w0。
  - 第44行调用父类的coef_inte函数（见博文《最优化方法Python计算：无约束优化应用——线性回归模型》）计算决策面的系数coef_与截距intercept_。
  综合案例
  井字棋是很多人儿时喜爱的游戏。方形棋盘被井字形的4条直线分成9个空格。两个玩家各执标有一符：或为“$\times$”，或为“$\circ$”的棋子，轮流在棋盘中放置一子。一方所执符号棋子占据棋盘中一直线（水平或竖直或对角）上的三个空格（见题图），即胜出。文件tic-tac-toe.csv（来自UC Irvine Machine Learning Repository）罗列了958例棋盘格局及其胜负判断

x	x	x	x	o	o	x	o	o	positive
x	x	x	x	o	o	o	x	o	positive
$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
b	b	b	b	o	o	x	x	x	positive
x	x	o	x	x	o	o	b	o	negative
$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
o	x	o	o	x	x	x	o	x	negative
o	o	x	x	x	o	o	x	x	negative

文件中，每一行表示一种棋盘格局及对该格局胜负方的判定。前9个数据表示按行优先排列的棋盘格局，第10个数据“positive”表示执“$\times$”者胜，“negative”表示执“$\circ$”者胜。例如，第一行中前9个数据“x x x x o o x o o”表示棋盘格局

由于执“$\times$”者占满了第1行和第1列，故判断执“$\times$”者胜出，即第10个数据为“positive”。注意，格局数据中“b”表示对应格子未置入任何棋子。例如，数据“b b b b o o x x x”表示棋盘格局

下列代码试图用文件中的70例数据训练一个ALineSvc支持向量机分类模型，用其余数据对其进行测试。

import numpy as np
data = np.loadtxt('tic-tac-toe.csv',		#读取数据文件delimiter = ',', dtype = str)
X = np.array(data)							#转换为数组
Y = X[:, 9]									#读取标签数据
X = np.delete(X, [9], axis = 1)				#去掉标签列
m, n=X.shape
print('共有%d个数据样本'%m)
for i in range(m):for j in range(n):if X[i, j] == 'x':X[i, j] = 1if X[i, j] == 'o':X[i, j] = -1if X[i, j] == 'b':X[i, j] = 0
X = X.astype(float)
Y = np.array([1 if y == 'positive' else		#类别数值化-1 for y in Y]).astype(int)
a = np.arange(m)							#数据项下标
m1=70										#训练集样本容量
np.random.seed(3489)
print('随机抽取%d个样本作为训练数据。'%(m1))
train = np.random.choice(a,m1,replace=False)
test = np.setdiff1d(a,train)				#检测数据下标
tictactoe = ALineSvc(C = 10)				#创建模型
tictactoe.fit(X[train], Y[train])			#训练模型
print('系数：', tictactoe.coef_)
print('截距：%.4f'%tictactoe.intercept_)
support = tictactoe.support_				#支持向量下标
print('支持向量：%s'%support)
acc=tictactoe.score(X[test], Y[test]) * 100
print('对其余%d个样本数据测试，正确率：%.2f'%(m - m1,acc) + '%')

程序的第2~6行从文件中读取数据并转换为数组X，并从中拆分出标签数据Y。由于表示格局的特征数据和表示胜负的标签数据都是字符型，应转换为数值类型。对特征数据，作映射：
$$\text{x}\to 1,\text{o}\to -1,\text{b}\to 0$$
第9~16行的嵌套for循环完成9个格局数据按上述对应关系的类型转换。第18~19行将标签数据按映射
$$\text{positive}\to 1,\text{negative}\to -1$$
进行类型转换。第20~25行从958个棋盘格局数据中随机选取70（m1）个作为训练数据集X[train]、Y[train]，其余888个作为训练数据集X[test]、Y[test]。第26行声明ALineSvc类对象tictactoe作为支持向量机分类模型，设置罚项系数C为10.0。第27行调用tictactoe的fit函数用X[train]、Y[train]对其进行训练。第32行调用tictactoe的score函数用X[test]、Y[test]对其进行测试。运行程序，输出

共有958个数据样本
随机抽取100个样本作为训练数据。
训练中...，稍候
15次迭代后完成训练。
系数： [2. 2. 2. 2. 2. 2. 2. 2. 2.]
截距：-19.0286
支持向量：[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 2324 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 4748 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69]
对其余888个样本数据测试，正确率：98.31%

意为用70个训练数据，经过15次迭代完成训练。对其余888个数据对训练所得模型测试，正确率为98.31%。训练给出了决策面
$$\pi :2x_1+2x_2+2x_3+2x_4+2x_5+2x_6+2x_7+2x_8+2x_9-19.03=0$$

写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！