前言

• 在支持向量机中，理论逻辑很简单：最大化最小的几何间隔。但是实际编写代码过程中有一个小点需要注意。总是把二分类的类别分为0和1，这样就导致我的目标函数跟算法描述的就不一样，所以求解结果就不正确。
• 同时还有第二个要注意的就是凸优化包cvxpy中各种运算的表示方法，比如凸优化中常见的二次方程的表示，变量的默认形状等，要查看官方文档才能熟悉。变量的默认形状为列向量。

参考：

• 李航：统计学习方法
• CVXPY

分析

支持向量机算法中，我们的训练数据除了是两种类别以外，类别的编号也有要求，分别是1、-1，只有这样，我们才能求每个样本所对应的函数间隔$D=y_i(w{x}_i+b)$和几何间隔$D=y_i\frac{w{x}_i+b}{\Vert w{\Vert }_2}$，在这种类别标签的情况下，预测值$w{x}_i+b$与真实值$y_i$之间的乘积才有意义。
预测与真实相同，乘积才会是大于0的；预测与真实相反，乘积就是小于0的。只有这样后面的目标函数最大化几何间隔才有意义。
重新表述一下可分数据集上支持向量机的目标函数和约束条件：
$$\begin{array}{rlr}& \text{max}& \lambda \\ & \text{s.t.}& y_i\frac{w{x}_i+b}{\Vert w{\Vert }_2}\ge \lambda \end{array}$$

如果我们使用类别标签为0，1，那么当错误分类时，几何间隔为0，无法指导参数修改。所以必须要使用1，-1。

数据集线性可分情况下的支持向量机

此时有两种求法，一种是用原始算法，直接用cvxpy函数包求解原始的这个凸优化问题，并把问题变为下述形式：
$$\begin{array}{rlr}& \text{max}& \frac{t}{\Vert w{\Vert }_2}\\ & \text{s.t.}& y_i(w{x}_i+b)\ge t\end{array}$$

由于通过同比例放大w,b可以实现条件中左边的乘积大小的任意变换，所以我们修改t为$1$。上述凸优化问题就变为：
$$\begin{array}{rlr}& \text{min}& \Vert w{\Vert }_2\\ & \text{s.t.}& y_i(w{x}_i+b)\ge 1\end{array}$$

第二种方法就是将原始问题使用拉格朗日乘子法变换为对偶问题，将加入条件和拉格朗日乘子的拉格朗日函数进行求导，并将求导得到的关系式带入拉格朗日函数，这样就可以得到对偶问题。

原始问题凸优化包解法

import numpy as np
import cvxpy as cp
#生成符合要求的样本数据
np.random.seed(3)
train_x=np.random.randn(2,2)
weight_x=np.random.randn(2)
bias_x=np.random.randn(1)
train_y=np.where(train_x@weight_x+bias_x<0,-1,1).reshape(2,-1)
print(train_x,train_y)
#求解对偶凸优化问题
w=cp.Variable(2)
b=cp.Variable(1)
obj=cp.Minimize(1/2*cp.sum_squares(w))
cons=[train_y[0]*(train_x@w+b)[0]>=1,train_y[1]*(train_x@w+b)[1]>=1]
prob=cp.Problem(obj,cons)
result=prob.solve()
#输出拉格朗日乘子的和最优化的目标函数值以及权重w
print(f'w.value,result,b.value:{w.value,result,b.value}')

结果如下：

[[ 1.78862847  0.43650985][ 0.09649747 -1.8634927 ]] [[-1][ 1]]
w.value,result,b.value:(array([-0.41507783, -0.56418804]), 0.24529887505030906, array([-0.01130633]))

对偶问题凸优化包解法

本例子中我们使用了小规模的数据，只有两个样本，所以这两个样本肯定都是支持向量，也就是对应的拉格朗日乘子都不为0，对于大规模样本数据的情况，如果不在分界面上，那么对应的拉格朗日乘子为0，也就不是支持向量。拉格朗日乘子不为0的肯定就是支持向量。

import numpy as np
import cvxpy as cp
#生成符合要求的样本数据
np.random.seed(3)
train_x=np.random.randn(2,2)
weight_x=np.random.randn(2)
bias_x=np.random.randn(1)
train_y=np.where(train_x@weight_x+bias_x<0,-1,1).reshape(2,-1)
print(train_x,train_y)
#求解对偶凸优化问题
alpha=cp.Variable(2)
obj=cp.Minimize(1/2*cp.quad_form(alpha,(train_x@train_x.T)*(train_y@train_y.T))-cp.sum(alpha))
cons=[alpha>=0,train_y.T@alpha>=0,train_y.T@alpha<=0]
prob=cp.Problem(obj,cons)
result=prob.solve()
#输出拉格朗日乘子的和最优化的目标函数值以及权重w
print(f'alpha.value,result,w:{alpha.value,result,np.array(alpha.value.reshape(2,-1)*train_y*train_x).sum(axis=0)}')
#检验支持向量机求出的分离面是否与这两个样本之间的连线垂直
w=np.array(alpha.value.reshape(2,-1)*train_y*train_x).sum(axis=0)
np.array(train_x[0,:]-train_x[1,:]).reshape(1,-1)@np.array([-w[1]/w[0],1]).reshape(2,-1)

结果如下：

[[ 1.78862847  0.43650985][ 0.09649747 -1.8634927 ]] [[-1][ 1]]
alpha.value,result,w:(array([0.24529888, 0.24529888]), -0.24529887505030898, array([-0.41507783, -0.56418804]))
array([[-4.4408921e-16]])

至于参数b，我们可以通过支持向量所对应的等式求出。
我们可以观察一下原始问题与对偶问题的解答是否一致。

数据集线性不可分情况下的线性支持向量机与软间隔最大化

软间隔顾名思义就是给原来的间隔留下一点宽容量，给那些不容易分正确的留一点余地。同时对于这些余地进行惩罚所得到的分割面，以上面线性可分的数据做演示。

import numpy as np
import cvxpy as cp
#生成符合要求的样本数据
np.random.seed(3)
train_x=np.random.randn(2,2)
weight_x=np.random.randn(2)
bias_x=np.random.randn(1)
train_y=np.where(train_x@weight_x+bias_x<0,-1,1).reshape(2,-1)
print(train_x,train_y)
#求解对偶凸优化问题
w=cp.Variable(2)
b=cp.Variable(1)
softgap=cp.Variable(2)
obj=cp.Minimize(1/2*cp.sum_squares(w)+100*cp.sum(softgap))
cons=[train_y[0]*(train_x@w+b)[0]>=1-softgap[0],train_y[1]*(train_x@w+b)[1]>=1-softgap[1],softgap>=0]
prob=cp.Problem(obj,cons)
result=prob.solve()
#输出拉格朗日乘子的和最优化的目标函数值以及权重w
print(f'w.value,result,b.value,softgap.value:{w.value,result,b.value,softgap.value}')

结果如下：

[[ 1.78862847  0.43650985][ 0.09649747 -1.8634927 ]] [[-1][ 1]]
w.value,result,b.value,softgap.value:(array([-0.41507783, -0.56418804]), 0.24529887505030906, array([-0.01130633]), array([-2.90746355e-22,  1.53881391e-22]))

可以看出在线性可分的的情况下软间隔不起作用。
那么我们制造一些线性不可分的数据，来测试一下。

import numpy as np
import cvxpy as cp
#生成符合要求的样本数据
np.random.seed(3)
train_x=np.array([[0,0],[1,0],[2,0]])
train_y=np.array([-1,1,-1]).reshape(3,-1)
print(train_x,train_y)
#求解对偶凸优化问题
w=cp.Variable(2)
b=cp.Variable(1)
softgap=cp.Variable(3)
obj=cp.Minimize(1/2*cp.sum_squares(w)+0.1*cp.sum(softgap))
cons=[train_y[0]*(train_x@w+b)[0]>=1-softgap[0],train_y[1]*(train_x@w+b)[1]>=1-softgap[1],train_y[2]*(train_x@w+b)[2]>=1-softgap[2],softgap>=0]
prob=cp.Problem(obj,cons)
result=prob.solve()
#输出拉格朗日乘子的和最优化的目标函数值以及权重w
print(f'w.value,result,b.value,softgap.value:{w.value,result,b.value,softgap.value}')

结果如下：

[[0 0][1 0][2 0]] [[-1][ 1][-1]]
w.value,result,b.value,softgap.value:(array([9.07653476e-18, 0.00000000e+00]), 0.2, array([-1.]), array([ 8.59013373e-23,  2.00000000e+00, -8.59013423e-23]))