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

更多的2-分类问题$\{({\mathbf{x}}_i,y_i)\}$，${\mathbf{x}}_i\in {\text{R}}^n$， y i ∈ { − 1 , 1 } y_i\in\{-1,1\} yi​∈{−1,1}，$i=1,2,\cdots ,m$不是线性可分的，即不存在超平面将 { x i ∣ i = 1 , 2 , ⋯ , m } \{\boldsymbol{x}_i|i=1,2,\cdots,m\} {xi​∣i=1,2,⋯,m}按对应$y_i$的不同取值而分隔开（如下图(a)）。解决这类问题的方法是进行坐标变换，将样本特征数据 { x i } \{\boldsymbol{x}_i\} {xi​}通过合适的映射$\mathbf{\varphi }:{\text{R}}^n\to {\text{R}}^{n^{\prime }}$（$n^{\prime }\ge n$），在一个维度更高的空间内将问题转换为线性可分问题（如下图(b)）。

即求解有约束优化问题
$$\{\begin{array}{l}\min \frac{1}{2}{\mathbf{w}}^{\top }\mathbf{Hw}+{\mathbf{c}}^{\top }\varphi (\mathbf{z})\\ \text{s.t. }{y}_i(\varphi ({\mathbf{x}}_i{)}^{\top },1)\mathbf{w}\ge 1-\varphi (z_i)\\ \varphi (\mathbf{z})\ge \mathbf{o}\end{array}i=1,2,\cdots m(1)$$
的对偶问题
$$\{\begin{array}{l}\min \frac{1}{2}{\mathbf{\mu }}^{\top }{\mathbf{Q}}^{\prime }\mathbf{\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{Q}}^{\prime }=\left(\begin{array}{cccc}{y}_1{y}_1\mathbf{\varphi }({\mathbf{x}}_1{)}^{\top }\mathbf{\varphi }({\mathbf{x}}_1)& y_1{y}_2\mathbf{\varphi }({\mathbf{x}}_1{)}^{\top }\mathbf{\varphi }({\mathbf{x}}_2)& \cdots & y_1{y}_m\mathbf{\varphi }({\mathbf{x}}_1{)}^{\top }\mathbf{\varphi }({\mathbf{x}}_m)\\ y_2{y}_1\mathbf{\varphi }({\mathbf{x}}_2{)}^{\top }\mathbf{\varphi }({\mathbf{x}}_1)& y_2{y}_2\mathbf{\varphi }({\mathbf{x}}_2{)}^{\top }\mathbf{\varphi }({\mathbf{x}}_2)& \cdots & y_2{y}_m\mathbf{\varphi }({\mathbf{x}}_2{)}^{\top }\mathbf{\varphi }({\mathbf{x}}_m)\\ ⋮& ⋮& \ddots & ⋮\\ y_m{y}_1\mathbf{\varphi }({\mathbf{x}}_m{)}^{\top }\mathbf{\varphi }({\mathbf{x}}_1)& y_m{y}_2\mathbf{\varphi }({\mathbf{x}}_m{)}^{\top }\mathbf{\varphi }({\mathbf{x}}_2)& \cdots & y_m{y}_m\mathbf{\varphi }({\mathbf{x}}_m{)}^{\top }\mathbf{\varphi }({\mathbf{x}}_m)\end{array}\right)$$
直接寻求映射$\mathbf{\varphi }:{\text{R}}^n\to {\text{R}}^{n^{\prime }}$（$n^{\prime }\ge n$）并计算矩阵${\mathbf{Q}}^{\prime }$会使计算量爆炸式增长。克服这一困难的方法之一是所谓的“核函数技巧”：
将$\mathbf{\varphi }(\mathbf{x}{)}^{\top }\mathbf{\varphi }(\mathbf{y})$视为函数$\kappa :{\text{R}}^{n^{\prime }}\times {\text{R}}^{n^{\prime }}\to \text{R}$，则对 { x i ∣ x i ∈ R n , i = 1 , 2 , ⋯ , m } \{\boldsymbol{x}_i|\boldsymbol{x}_i\in\text{R}^n,i=1,2,\cdots,m\} {xi​∣xi​∈Rn,i=1,2,⋯,m}，$\kappa (\mathbf{x},\mathbf{y})$若矩阵 $${\mathbf{Q}}_{\kappa }=\left(\begin{array}{cccc}\kappa ({\mathbf{x}}_1,{\mathbf{x}}_1)& \kappa ({\mathbf{x}}_1,{\mathbf{x}}_2)& \cdots & \kappa ({\mathbf{x}}_1,{\mathbf{x}}_m)\\ \kappa ({\mathbf{x}}_2,{\mathbf{x}}_1)& \kappa ({\mathbf{x}}_2,{\mathbf{x}}_2)& \cdots & \kappa ({\mathbf{x}}_2,{\mathbf{x}}_m)\\ ⋮& ⋮& \ddots & ⋮\\ \kappa ({\mathbf{x}}_m,{\mathbf{x}}_1)& \kappa ({\mathbf{x}}_m,{\mathbf{x}}_2)& \cdots & \kappa ({\mathbf{x}}_m,{\mathbf{x}}_m)\end{array}\right)$$
半正定，称$\kappa (\mathbf{x},\mathbf{y})$为核函数。事实上向量间的内积$\kappa (\mathbf{x},\mathbf{y})={\mathbf{x}}^{\top }\mathbf{y}$就是一个核函数，内积的多项式$\kappa (\mathbf{x},\mathbf{y})=(\gamma {\mathbf{x}}^{\top }\mathbf{y}+1{)}^d$，$\gamma >0$，$d>0$也是一个核函数，径向基函数$\kappa (\mathbf{x},\mathbf{y})=e^{-\gamma \Vert \mathbf{x}-\mathbf{y}{\Vert }^2}$（$\gamma >0$）被称为万能核函数。对给定的核函数$\kappa (\mathbf{x},\mathbf{y})$，使得$\kappa (\mathbf{x},\mathbf{y})=\varphi (\mathbf{x}{)}^{\top }\varphi (\mathbf{y})$，二次规划(2)中的
$${\mathbf{Q}}^{\prime }=\left(\begin{array}{cccc}{y}_1{y}_1\kappa ({\mathbf{x}}_1,{\mathbf{x}}_1)& y_1{y}_2\kappa ({\mathbf{x}}_1,{\mathbf{x}}_2)& \cdots & y_1{y}_m\kappa ({\mathbf{x}}_1,{\mathbf{x}}_m)\\ y_2{y}_1\kappa ({\mathbf{x}}_2,{\mathbf{x}}_1)& y_2{y}_2\kappa ({\mathbf{x}}_2,{\mathbf{x}}_2)& \cdots & y_2{y}_m\kappa ({\mathbf{x}}_1,{\mathbf{x}}_m)\\ ⋮& ⋮& \ddots & ⋮\\ y_m{y}_1\kappa ({\mathbf{x}}_m,{\mathbf{x}}_1)& y_m{y}_2\kappa ({\mathbf{x}}_m,{\mathbf{x}}_2)& \cdots & y_m{y}_m\kappa ({\mathbf{x}}_m,{\mathbf{x}}_m)\end{array}\right),$$
则(2)为二次规划问题。设其最优解${\mathbf{\mu }}_0$，计算支持向量下标集 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|$，可算得参数
$$b_0=\frac{1}{m_s}\sum _{j\in s}\left(y_j-\sum _{i\in s}{\mu }_i{y}_i\kappa ({\mathbf{x}}_i,{\mathbf{x}}_j)\right)$$
进而算得决策函数
$$F(\mathbf{x})=\sum _{i\in s}{\mu }_{0_i}{y}_i\kappa (\mathbf{x},{\mathbf{x}}_i)+b_0$$
对新的样本特征数据$\mathbf{x}$，分类预测值为
$$\text{sign}(F(\mathbf{x})).$$
下列代码就以径向基函数为核函数实现先行不可分问题的支持向量机模型。

import numpy as np														#导入numpy
from scipy.optimize import minimize, LinearConstraint					#导入minimize等
rbf = lambda x, y, gamma: np.exp(-gamma * np.linalg.norm(x - y) ** 2)	#径向基函数
class SvmClassifier(Classification, SupervisedLearningModel):			#支持向量机分类器def __init__(self, C = 1e+4):										#构造函数self.kernel = rbfself.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)sigma2 = X.var() * nself.gamma = 1.0 / sigma2Qk = np.array([[self.kernel(self.X[i], self.X[j], self.gamma)for j in range(m)]for i in range(m)])self.Q = np.outer(self.Y, self.Y) * Qkh = lambda x: self.Y @ x									#等式约束函数g1 = lambda x: identity(x)									#不等式约束函数1g2 = lambda x: self.C - identity(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.b0 = (self.Y[self.support_]-(self.mu0[self.support_] * self.Y[self.support_])\@ Qk[np.ix_(self.support_, self.support_)]).mean()print("%d次迭代后完成训练。"%res.nit)def F(self, w= None, X = None):									#决策函数self.K = lambda x: np.array([[self.kernel(x[i], self.X[j], self.gamma)for j in self.support_]for i in range(x.shape[0])])return self.K(X) @ (self.mu0[self.support_] * self.Y[self.support_]) + self.b0

程序中，第3行定义径向基函数$rbf(\mathbf{x},\mathbf{y})=e^{-\gamma \Vert \mathbf{x}-\mathbf{y}{\Vert }^2}$。由于决策函数不再是线性的，故第4~48行定义的支持向量机分类器类SvmClassifier继承了SupervisedLearningModel（见博文《最优化方法Python计算：无约束优化应用——线性回归模型》）和Classification（见博文《最优化方法Python计算：无约束优化应用——线性回归分类器》）的属性与方法。类定义体中

• 第5~8行定义构造函数，C参数表示惩罚系数，缺省值为$10^4$。kernel表示核函数，第6行设置为rbf。第7行将标签值函数tagVal设置为Numpy的sign函数，以便计算分类预测值。
• 第9~10行定义目标函数obj，返回值为$\frac{1}{2}{\mathbf{\mu }}^{\top }{\mathbf{Q}}^{\prime }\mathbf{\mu }-{\mathbf{\mu }}^{\top }1$，即问题(2)的目标函数。
• 第11~15行重载标签数据归一化函数ynormalize。由于支持向量机模型中的标签数据不需要归一化，所以在训练时将self.ymin和self.ymax设置为0和1。第13行返回归一化后的标签数据。做了这样的调整，我们在进行预测操作时，式
  $$y=y\cdot (\max y-\min y)+\min y$$
  仍保持$y$的值不变，进而可保持SupervisedLearningModel的predict函数代码不变（见博文《最优化方法Python计算：无约束优化应用——线性回归模型》）。
• 第16~43行重载训练函数fit。对比程序4.1中定义的父类fit函数可见第17~25行的代码是保持一致的（仅将模型参数w改为mu）。
  - 第26~27行计算样本特征数据X的方差${\sigma }^2$，并将$\gamma$参数设置为$\frac{1}{n{\sigma }^2}$。其中，n为每个样本的特征个数。
  - 第28~30行计算核矩阵${\mathbf{Q}}_{\kappa }$。
  - 第31行计算问题(2)的目标函数系数矩阵${\mathbf{Q}}^{\prime }$。
  - 第32~34行分别定义问题(2)的等式约束条件$\mathbf{h}(\mathbf{\mu })={\mathbf{y}}^{\top }\mathbf{\mu }=\mathbf{o}$中的函数$\mathbf{h}(\mathbf{x})$；不等式约束函数${\mathbf{g}}_1(\mathbf{\mu })=\mathbf{\mu }\ge \mathbf{o}$中的函数${\mathbf{g}}_1(\mathbf{x})$；不等式约束函数${\mathbf{g}}_2(\mathbf{\mu })=\mathbf{c}-\mathbf{\mu }\ge \mathbf{o}$中的函数${\mathbf{g}}_2(\mathbf{x})$。注意，$\mathbf{h}(\mathbf{x})$、${\mathbf{g}}_1(\mathbf{x})$和${\mathbf{g}}_2(\mathbf{x})$均为向量函数。
  - 第35~37行定义问题(2)的约束条件列表cons。
  - 第38行调用minimize函数求解优化问题(2)，返回值赋予res。第39行将res的x属性赋予mu0。
  - 第40行按条件$0<{\mu }_{0_i}<C$计算支持向量下标集support_。
  - 第41~42行按式
  $$b_0=\frac{1}{m_s}\sum _{j\in s}\left(y_j-\sum _{i\in s}{\mu }_i{y}_i\kappa ({\mathbf{x}}_i,{\mathbf{x}}_j)\right)$$
  计算参数$b_0$。
• 第44~48行定义决策函数F。该函数的参数w和X分别表示权重向量和样本特征数据。由于支持向量机模型中不需要权重向量，所以该参数w没有实际意义。第45~47行定义辅助函数K(x)，第48行利用K(x)按式
  $$\sum _{i\in s}{\mu }_{0_i}{y}_i\kappa (\mathbf{x},{\mathbf{x}}_i)+b_0$$
  计算决策函数F(X)的值。
  例1 文件make_moons.csv（需要者可私信）包含了1000个3维空间的点，

X	Y	label
-1.129064193	-0.049046872	0
…	…	…
1.275997495	-0.424925687	1
…	…	…
-0.360415936	1.014813199	0

这些点投影到XY平面的图形如上图(a)所示。标签label取值0/1，用支持向量机训练分类器，需要将取值转换为-1/1。下列代码用该数据集中的66个训练一个SvmClassifier类对象，并用剩余数据对其进行测试。

import numpy as np										#导入numpy
data = np.loadtxt('make_moons.csv', delimiter=',',dtype = object, skiprows = 1)						#读取数据文件
X = np.array(data).astype(float)						#转换为数组
Y = np.array([-1.0 if y == 0 else 1.0 for y in X[:,2]])	#标签数据
X = np.delete(X, [2], axis = 1)
m = X.shape[0]											#样本数
print('共有%d个数据样本'%m)
a = np.arange(m)
np.random.seed(1052)									#随机种子
print('用%d个样本数据训练模型'%(m // 15))
train = np.random.choice(a, m // 15, replace = False)	#训练数据下标
test = np.setdiff1d(a, train)							#测试数据下标
moons = SvmClassifier()									#构造模型
moons.fit(X[train], Y[train])							#训练模型
print('支持向量：%s'%moons.support_)
acc=moons.score(X[test], Y[test])						#测试模型
print('用其余%d个样本测试模型，正确率为%.2f'%(m - m // 15, acc * 100)+ '%')

程序的第2~6行从文件中读取数据并转换为数组X，并从中拆分出标签数据Y。第7行读取样本数m。第12行随机选取m//15个样本作为训练数据集X[train]、Y[train]，其余样本作为测试数据集X[test]、Y[test]。第14行声明SvmClassifier类对象moons，作为支持向量机分类模型。注意，罚项系数C使用缺省值$10^4$。第15行调用moons的fit函数用X[train]、Y[train]对其进行训练。第16行输出支持向量下标集support_。第17行调用moons的score函数用X[test]、Y[test]对其进行测试。运行程序，输出

共有1000个数据样本
用66个样本数据训练模型
训练中...，稍候
100次迭代后完成训练。
支持向量：[ 2  6  8 16 25 26 27 29 32 34 37 46 53 55 59 62]
用其余934个样本测试模型，正确率为100.00%

训练所得模型测试效果不错！
例2下列代码用井字棋数据文件tic-tac-toe.csv（数据集结构见博文《最优化方法Python计算：有约束优化应用——近似线性可分问题支持向量机》）的数据训练并测试由SvmClassifier实现的支持向量机分类模型，罚项系数C设置为$10^4$。

import numpy as np														#导入numpy
data = np.loadtxt('tic-tac-toe.csv', delimiter = ',', dtype = object)	#读取数据文件
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=100
np.random.seed(1264)													#随机种子
print('随机抽取%d个样本作为训练数据。'%(m1))
train=np.random.choice(a,m1,replace=False)								#训练数据下标
test = np.setdiff1d(a,train)											#测试数据下标
tictactoe = SvmClassifier(C = 100.0)									#创建模型
tictactoe.fit(X[train],Y[train])										#训练模型
print('支持向量：%s'%tictactoe.support_)
acc=tictactoe.score(X[test], Y[test]) * 100								#测试模型
print('对其余%d个样本数据测试，正确率：%.2f'%(m-m1,acc)+'%')

程序的第2~18行从文件中读取数据并转换为数组X，并从中拆分出标签数据Y。第19~24行在数据集中随机选取100个作为训练用数据X[train]，Y[train]，其余的作为测试用数据X[test]，Y[test]。第25行声明SvmClassifier类对象tictactoc，创建核函数为rbf，正则化参数为$10^2$的支持向量机模型。第26行训练该模型，第28行测试模型。运行程序，输出

共有958个数据样本
随机抽取100个样本作为训练数据。
训练中...，稍候
100次迭代后完成训练。
支持向量：[ 1  2  3  5  7  9 10 11 14 15 16 17 18 19 20 22 23 24 25 26 29 31 33 3435 37 40 43 46 47 48 49 50 51 52 55 58 60 61 62 67 69 70 72 76 78 82 8386 89 90 93 94 95 96 98 99]
对其余858个样本数据测试，正确率：99.30%

与博文《最优化方法Python计算：有约束优化应用——近似线性可分问题支持向量机》中同一问题的计算结果比较，可见对此问题SvmClassifier类模型优于ALineSvc类模型。
写博不易，敬请支持：
如果阅读本文于您有所获，敬请点赞、评论、收藏，谢谢大家的支持！