最优化方法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 = C
        self.tagVal = np.sign
    def obj(self, mu):													#目标函数
        return 0.5 * (mu @ (self.Q @ mu)) - mu.sum()
    def ynormalize(self, y, trained):									#标签数据预处理
        if not trained:
            self.ymin = 0
            self.ymax = 1
        return (y - self.ymin) / (self.ymax - self.ymin)
    def fit(self, X, Y, mu = None):										#训练函数
        print("训练中...，稍候")
        m, n = X.shape
        self.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												#不等式约束函数1
        g2 = lambda x: self.C - x										#不等式约束函数2
        cons = [{'type': 'eq', 'fun': h},								#约束条件列表
                {'type': 'ineq', 'fun': g1},
                {'type': 'ineq', 'fun': g2}]
        res = minimize(self.obj, mu, constraints = cons)				#求解最优化问题
        self.mu0 = res.x
        self.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函数代码不变（见博文《最优化方法Python计算：无约束优化应用——线性回归分类器》）。
• 第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] = 1
        if X[i, j] == 'o':
            X[i, j] = -1
        if 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 23
 	     24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
 	     48 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$$

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