Smooth Support Vector Machine - Python实现

• 算法特征:
  ①. 所有点尽可能正确分开; ②. 极大化margin; ③. 极小化非线性可分之误差.
• 算法推导:
  Part Ⅰ
  线性可分之含义:
  包含同类型所有数据点的最小凸集合彼此不存在交集.
  引入光滑化手段:
  plus function: \begin{equation*}(x)_{+} = \max \{ x, 0 \}\end{equation*}
  p-function:  \begin{equation*}p(x, \beta) =  \frac{1}{\beta}\ln(\mathrm{e}^{\beta x} + 1)\end{equation*}
  可采用p-function近似plus function. 相关函数图像如下:
  下面给出p-function的1阶及2阶导数:
  1st order (sigmoid function): \begin{equation*} s(x, \beta) = \frac{1}{\mathrm{e}^{-\beta x} + 1} \end{equation*}
  2nd order (delta function): \begin{equation*} \delta (x, \beta) = \frac{\beta \mathrm{e}^{\beta x}}{(\mathrm{e}^{\beta x} + 1)^2} \end{equation*}
  相关函数图像如下:
  Part Ⅱ
  SVM决策超平面方程如下:
  \begin{equation}\label{eq_1}
  h(x) = W^{\mathrm{T}}x + b
  \end{equation}
  其中, $W$是超平面法向量, $b$是超平面偏置参数.
  本文拟采用$L_2$范数衡量非线性可分之误差, 则SVM原始优化问题如下:
  \begin{equation}
  \begin{split}
  \min &\quad\frac{1}{2}\left\|W\right\|_2^2 + \frac{c}{2}\left\|\xi\right\|_2^2 \\
  \mathrm{s.t.} &\quad D(A^{\mathrm{T}}W + \mathbf{1}b) \geq \mathbf{1} - \xi
  \end{split}\label{eq_2}
  \end{equation}
  其中, $A = [x^{(1)}, x^{(2)}, \cdots , x^{(n)}]$, $D = diag([\bar{y}^{(1)}, \bar{y}^{(2)}, \cdots, \bar{y}^{(n)}])$, $\xi$是由所有样本非线性可分之误差组成的列矢量.
  根据上述优化问题(\ref{eq_2}), 误差$\xi$可采用plus function等效表示:
  \begin{equation}\label{eq_3}
  \xi = (\mathbf{1} - D(A^{\mathrm{T}}W + \mathbf{1}b))_{+}
  \end{equation}
  由此可将该有约束最优化问题转化为如下无约束最优化问题:
  \begin{equation}\label{eq_4}
  \min\quad \frac{1}{2}\left\|W\right\|_2^2 + \frac{c}{2}\left\| \left(\mathbf{1} - D(A^{\mathrm{T}}W + \mathbf{1}b)\right)_{+} \right\|_2^2
  \end{equation}
  同样, 根据优化问题(\ref{eq_2})KKT条件中稳定性条件有:
  \begin{equation}\label{eq_5}
  W = AD^{\mathrm{T}}\alpha
  \end{equation}
  其中, $\alpha$为对偶变量. 代入式(\ref{eq_4})变数变换可得:
  \begin{equation}\label{eq_6}
  \min\quad \frac{1}{2}\alpha^{\mathrm{T}}D^{\mathrm{T}}A^{\mathrm{T}}AD\alpha + \frac{c}{2}\left\| \left(\mathbf{1} - D(A^{\mathrm{T}}AD^{\mathrm{T}}\alpha + \mathbf{1}b)\right)_{+} \right\|_2^2
  \end{equation}
  由于权重参数$c$的存在, 上述最优化问题等效如下多目标最优化问题:
  \begin{equation}
  \left\{
  \begin{split}
  &\min\quad \frac{1}{2}\alpha^{\mathrm{T}}D^{\mathrm{T}}A^{\mathrm{T}}AD\alpha \\
  &\min\quad \frac{1}{2}\left\| \left(\mathbf{1} - D(A^{\mathrm{T}}AD^{\mathrm{T}}\alpha + \mathbf{1}b)\right)_{+} \right\|_2^2
  \end{split}
  \right.
  \label{eq_7}
  \end{equation}
  由于$D^TA^TAD\succeq 0$, 有如下等效:
  \begin{equation}\label{eq_8}
  \min\quad \frac{1}{2}\alpha^{\mathrm{T}}D^{\mathrm{T}}A^{\mathrm{T}}AD\alpha \quad\Rightarrow\quad \min\quad \frac{1}{2}\alpha^{\mathrm{T}}\alpha
  \end{equation}
  根据Part Ⅰ中光滑化手段, 有如下等效:
  \begin{equation}\label{eq_9}
  \min\quad \frac{1}{2}\left\| \left(\mathbf{1} - D(A^{\mathrm{T}}AD^{\mathrm{T}}\alpha + \mathbf{1}b)\right)_{+} \right\|_2^2 \quad\Rightarrow\quad \min\quad \frac{1}{2}\left\| p(\mathbf{1} - DA^{\mathrm{T}}AD\alpha - D\mathbf{1}b, \beta)\right\|_2^2, \quad\text{where $\beta\to\infty$}
  \end{equation}
  结合式(\ref{eq_8})、式(\ref{eq_9}), 优化问题(\ref{eq_6})等效如下:
  \begin{equation}\label{eq_10}
  \min\quad \frac{1}{2}\alpha^{\mathrm{T}}\alpha + \frac{c}{2}\left\| p(\mathbf{1} - DA^{\mathrm{T}}AD\alpha - D\mathbf{1}b, \beta)\right\|_2^2, \quad\text{where $\beta\to\infty$}
  \end{equation}
  为增强模型的分类能力, 引入kernel trick, 得最终优化问题:
  \begin{equation}\label{eq_11}
  \min\quad \frac{1}{2}\alpha^{\mathrm{T}}\alpha + \frac{c}{2}\left\| p(\mathbf{1} - DK(A^{\mathrm{T}}, A)D\alpha - D\mathbf{1}b, \beta)\right\|_2^2, \quad\text{where $\beta\to\infty$}
  \end{equation}
  其中, $K$代表kernel矩阵, 内部元素$K_{ij}=k(x^{(i)}, x^{(j)})$由kernel函数决定. 常见的kernel函数如下:
  ①. 多项式核: $k(x^{(i)}, x^{(j)}) = ({x^{(i)}}\cdot x^{(j)})^n$
  ②. 高斯核: $k(x^{(i)}, x^{(j)}) = \mathrm{e}^{-\mu\| x^{(i)} - x^{(j)}\|_2^2}$
  结合式(\ref{eq_1})、式(\ref{eq_5})及kernel trick可得最终决策超平面方程:
  \begin{equation}
  h(x) = \alpha^{\mathrm{T}}DK(A^{\mathrm{T}}, x) + b
  \end{equation}
  
  Part Ⅲ
  以下给出优化相关协助信息. 拟设式(\ref{eq_11})之目标函数符号为$J$, 并令:
  \begin{equation}
  J = J_1 + J_2\quad\quad\text{其中,}
  \left\{
  \begin{split}
  &J_1 = \frac{1}{2}\alpha^{\mathrm{T}}\alpha   \\
  &J_2 = \frac{c}{2}\left\| p(\mathbf{1} - DK(A^{\mathrm{T}}, A)D\alpha - D\mathbf{1}b, \beta)\right\|_2^2
  \end{split}
  \right.
  \end{equation}
  $J_1$相关:
  \begin{gather*}
  \frac{\partial J_1}{\partial \alpha} = \alpha \quad \frac{\partial J_1}{\partial b} = 0  \\
  \frac{\partial^2J_1}{\partial\alpha^2} = I \quad \frac{\partial^2J_1}{\partial\alpha\partial b} = 0 \quad \frac{\partial^2J_1}{\partial b\partial\alpha} = 0 \quad \frac{\partial^2J_1}{\partial b^2} = 0
  \end{gather*}
  \begin{gather*}
  \nabla J_1 = \begin{bmatrix} \frac{\partial J_1}{\partial \alpha} \\ \frac{\partial J_1}{\partial b}\end{bmatrix} \\ \\
  \nabla^2 J_1 = \begin{bmatrix}\frac{\partial^2J_1}{\partial\alpha^2} & \frac{\partial^2J_1}{\partial\alpha\partial b} \\ \frac{\partial^2J_1}{\partial b\partial\alpha} & \frac{\partial^2J_1}{\partial b^2}\end{bmatrix}
  \end{gather*}
  $J_2$相关:
  \begin{gather*}
  z_i = 1 - \sum_{j}K_{ij}\bar{y}^{(i)}\bar{y}^{(j)}\alpha_j - \bar{y}^{(i)}b   \\
  \frac{\partial J_2}{\partial \alpha_k} = -c\sum_{i}\bar{y}^{(i)}\bar{y}^{(k)}K_{ik}p(z_i, \beta)s(z_i, \beta)   \\
  \frac{\partial J_2}{\partial b} = -c\sum_{i}\bar{y}^{(i)}p(z_i, \beta)s(z_i, \beta)   \\
  \frac{\partial^2 J_2}{\partial \alpha_k\partial \alpha_l} = c\sum_{i}\bar{y}^{(k)}\bar{y}^{(l)}K_{ik}K_{il}[s^2(z_i, \beta) + p(z_i, \beta)\delta(z_i, \beta)]   \\
  \frac{\partial^2 J_2}{\partial \alpha_k\partial b} = c\sum_{i}\bar{y}^{(k)}K_{ik}[s^2(z_i, \beta) + p(z_i, \beta)\delta(z_i, \beta)]   \\
  \frac{\partial^2 J_2}{\partial b\partial\alpha_l} = c\sum_{i}\bar{y}^{(l)}K_{il}[s^2(z_i, \beta) + p(z_i, \beta)\delta(z_i, \beta)]   \\
  \frac{\partial^2 J_2}{\partial b^2} = c\sum_{i}[s^2(z_i, \beta) + p(z_i, \beta)\delta(z_i, \beta)]
  \end{gather*}
  \begin{gather*}
  \nabla J_2 = \begin{bmatrix} \frac{\partial J_2}{\partial\alpha_k} \\ \frac{\partial J_2}{\partial b} \end{bmatrix}   \\ \\
  \nabla^2 J_2 = \begin{bmatrix} \frac{\partial^2 J_2}{\partial\alpha_k\partial\alpha_l} &  \frac{\partial^2 J_2}{\partial\alpha_k\partial b} \\ \frac{\partial^2 J_2}{\partial b\partial\alpha_l} & \frac{\partial^2 J_2}{\partial b^2}  \end{bmatrix}
  \end{gather*}
  目标函数$J$之gradient及Hessian如下:
  \begin{gather}
  \nabla J = \nabla J_1 + \nabla J_2   \\
  H = \nabla^2 J = \nabla^2 J_1 + \nabla^2 J_2
  \end{gather}
  Part Ⅳ
  现以二分类为例进行算法实施:
  \begin{equation*}
  \begin{cases}
  h(x) \geq 0 & y = 1 & \text{正例} \\
  h(x) < 0 & y = -1 & \text{负例}
  \end{cases}
  \end{equation*}
• 代码实现:
  

  # Smooth Support Vector Machine之实现
  
  import numpy
  from matplotlib import pyplot as plt
  
  
  def spiral_point(val, center=(0, 0)):
      rn = 0.4 * (105 - val) / 104
      an = numpy.pi * (val - 1) / 25
  
      x0 = center[0] + rn * numpy.sin(an)
      y0 = center[1] + rn * numpy.cos(an)
      z0 = -1
      x1 = center[0] - rn * numpy.sin(an)
      y1 = center[1] - rn * numpy.cos(an)
      z1 = 1
  
      return (x0, y0, z0), (x1, y1, z1)
  
  
  def spiral_data(valList):
      dataList = list(spiral_point(val) for val in valList)
      data0 = numpy.array(list(item[0] for item in dataList))
      data1 = numpy.array(list(item[1] for item in dataList))
      return data0, data1
  
  
  # 生成训练数据集
  trainingValList = numpy.arange(1, 101, 1)
  trainingData0, trainingData1 = spiral_data(trainingValList)
  trainingSet = numpy.vstack((trainingData0, trainingData1))
  # 生成测试数据集
  testValList = numpy.arange(1.5, 101.5, 1)
  testData0, testData1 = spiral_data(testValList)
  testSet = numpy.vstack((testData0, testData1))
  
  
  class SSVM(object):
  
      def __init__(self, trainingSet, c=1, mu=1, beta=100):
          self.__trainingSet = trainingSet                     # 训练集数据
          self.__c = c                                         # 误差项权重
          self.__mu = mu                                       # gaussian kernel参数
          self.__beta = beta                                   # 光滑化参数
  
          self.__A, self.__D = self.__get_AD()
  
  
      def get_cls(self, x, alpha, b):
          A, D = self.__A, self.__D
          mu = self.__mu
  
          x = numpy.array(x).reshape((-1, 1))
          KAx = self.__get_KAx(A, x, mu)
          clsVal = self.__calc_hVal(KAx, D, alpha, b)
          if clsVal >= 0:
              return 1
          else:
              return -1
  
  
      def get_accuracy(self, dataSet, alpha, b):
          '''
          正确率计算
          '''
          rightCnt = 0
          for row in dataSet:
              clsVal = self.get_cls(row[:2], alpha, b)
              if clsVal == row[2]:
                  rightCnt += 1
          accuracy = rightCnt / dataSet.shape[0]
          return accuracy
  
  
      def optimize(self, maxIter=100, epsilon=1.e-9):
          '''
          maxIter: 最大迭代次数
          epsilon: 收敛判据, 梯度趋于0则收敛
          '''
          A, D = self.__A, self.__D
          c = self.__c
          mu = self.__mu
          beta = self.__beta
  
          alpha, b = self.__init_alpha_b((A.shape[1], 1))
          KAA = self.__get_KAA(A, mu)
  
          JVal = self.__calc_JVal(KAA, D, c, beta, alpha, b)
          grad = self.__calc_grad(KAA, D, c, beta, alpha, b)
          Hess = self.__calc_Hess(KAA, D, c, beta, alpha, b)
  
          for i in range(maxIter):
              print("iterCnt: {:3d},   JVal: {}".format(i, JVal))
              if self.__converged1(grad, epsilon):
                  return alpha, b, True
  
              dCurr = -numpy.matmul(numpy.linalg.inv(Hess), grad)
              ALPHA = self.__calc_ALPHA_by_ArmijoRule(alpha, b, JVal, grad, dCurr, KAA, D, c, beta)
  
              delta = ALPHA * dCurr
              alphaNew = alpha + delta[:-1, :]
              bNew = b + delta[-1, -1]
              JValNew = self.__calc_JVal(KAA, D, c, beta, alphaNew, bNew)
              if self.__converged2(delta, JValNew - JVal, epsilon ** 2):
                  return alphaNew, bNew, True
  
              alpha, b, JVal = alphaNew, bNew, JValNew
              grad = self.__calc_grad(KAA, D, c, beta, alpha, b)
              Hess = self.__calc_Hess(KAA, D, c, beta, alpha, b)
          else:
              if self.__converged1(grad, epsilon):
                  return alpha, b, True
          return alpha, b, False
  
  
      def __converged2(self, delta, JValDelta, epsilon):
          val1 = numpy.linalg.norm(delta, ord=numpy.inf)
          val2 = numpy.abs(JValDelta)
          if val1 <= epsilon or val2 <= epsilon:
              return True
          return False
  
  
      def __calc_ALPHA_by_ArmijoRule(self, alphaCurr, bCurr, JCurr, gCurr, dCurr, KAA, D, c, beta, C=1.e-4, v=0.5):
          i = 0
          ALPHA = v ** i
          delta = ALPHA * dCurr
          alphaNext = alphaCurr + delta[:-1, :]
          bNext = bCurr + delta[-1, -1]
          JNext = self.__calc_JVal(KAA, D, c, beta, alphaNext, bNext)
          while True:
              if JNext <= JCurr + C * ALPHA * numpy.matmul(dCurr.T, gCurr)[0, 0]: break
              i += 1
              ALPHA = v ** i
              delta = ALPHA * dCurr
              alphaNext = alphaCurr + delta[:-1, :]
              bNext = bCurr + delta[-1, -1]
              JNext = self.__calc_JVal(KAA, D, c, beta, alphaNext, bNext)
          return ALPHA
  
  
      def __converged1(self, grad, epsilon):
          if numpy.linalg.norm(grad, ord=numpy.inf) <= epsilon:
              return True
          return False
  
  
      def __p(self, x, beta):
          term = x * beta
          if term > 10:
              val = x + numpy.math.log(1 + numpy.math.exp(-term)) / beta
          else:
              val = numpy.math.log(numpy.math.exp(term) + 1) / beta
          return val
  
  
      def __s(self, x, beta):
          term = x * beta
          if term > 10:
              val = 1 / (numpy.math.exp(-beta * x) + 1)
          else:
              term1 = numpy.math.exp(term)
              val = term1 / (1 + term1)
          return val
  
  
      def __d(self, x, beta):
          term = x * beta
          if term > 10:
              term1 = numpy.math.exp(-term)
              val = beta * term1 / (1 + term1) ** 2
          else:
              term1 = numpy.math.exp(term)
              val = beta * term1 / (term1 + 1) ** 2
          return val
  
  
      def __calc_Hess(self, KAA, D, c, beta, alpha, b):
          Hess_J1 = self.__calc_Hess_J1(alpha)
          Hess_J2 = self.__calc_Hess_J2(KAA, D, c, beta, alpha, b)
          Hess = Hess_J1 + Hess_J2
          return Hess
  
  
      def __calc_Hess_J2(self, KAA, D, c, beta, alpha, b):
          Hess_J2 = numpy.zeros((KAA.shape[0] + 1, KAA.shape[0] + 1))
          Y = numpy.matmul(D, numpy.ones((D.shape[0], 1)))
          YY = numpy.matmul(Y, Y.T)
          KAAYY = KAA * YY
  
          z = 1 - numpy.matmul(KAAYY, alpha) - Y * b
          p = numpy.array(list(self.__p(z[i, 0], beta) for i in range(z.shape[0]))).reshape((-1, 1))
          s = numpy.array(list(self.__s(z[i, 0], beta) for i in range(z.shape[0]))).reshape((-1, 1))
          d = numpy.array(list(self.__d(z[i, 0], beta) for i in range(z.shape[0]))).reshape((-1, 1))
          term = s * s + p * d
  
          for k in range(Hess_J2.shape[0] - 1):
              for l in range(k + 1):
                  val = c * Y[k, 0] * Y[l, 0] * numpy.sum(KAA[:, k:k+1] * KAA[:, l:l+1] * term)
                  Hess_J2[k, l] = Hess_J2[l, k] = val
              val = c * Y[k, 0] * numpy.sum(KAA[:, k:k+1] * term)
              Hess_J2[k, -1] = Hess_J2[-1, k] = val
          val = c * numpy.sum(term)
          Hess_J2[-1, -1] = val
          return Hess_J2
  
  
      def __calc_Hess_J1(self, alpha):
          I = numpy.identity(alpha.shape[0])
          term = numpy.hstack((I, numpy.zeros((I.shape[0], 1))))
          Hess_J1 = numpy.vstack((term, numpy.zeros((1, term.shape[1]))))
          return Hess_J1
  
  
      def __calc_grad(self, KAA, D, c, beta, alpha, b):
          grad_J1 = self.__calc_grad_J1(alpha)
          grad_J2 = self.__calc_grad_J2(KAA, D, c, beta, alpha, b)
          grad = grad_J1 + grad_J2
          return grad
  
  
      def __calc_grad_J2(self, KAA, D, c, beta, alpha, b):
          grad_J2 = numpy.zeros((KAA.shape[0] + 1, 1))
          Y = numpy.matmul(D, numpy.ones((D.shape[0], 1)))
          YY = numpy.matmul(Y, Y.T)
          KAAYY = KAA * YY
  
          z = 1 - numpy.matmul(KAAYY, alpha) - Y * b
          p = numpy.array(list(self.__p(z[i, 0], beta) for i in range(z.shape[0]))).reshape((-1, 1))
          s = numpy.array(list(self.__s(z[i, 0], beta) for i in range(z.shape[0]))).reshape((-1, 1))
          term = p * s
  
          for k in range(grad_J2.shape[0] - 1):
              val = -c * Y[k, 0] * numpy.sum(Y * KAA[:, k:k+1] * term)
              grad_J2[k, 0] = val
          grad_J2[-1, 0] = -c * numpy.sum(Y * term)
          return grad_J2
  
  
      def __calc_grad_J1(self, alpha):
          grad_J1 = numpy.vstack((alpha, [[0]]))
          return grad_J1
  
  
      def __calc_JVal(self, KAA, D, c, beta, alpha, b):
          J1 = self.__calc_J1(alpha)
          J2 = self.__calc_J2(KAA, D, c, beta, alpha, b)
          JVal = J1 + J2
          return JVal
  
  
      def __calc_J2(self, KAA, D, c, beta, alpha, b):
          tmpOne = numpy.ones((KAA.shape[0], 1))
          x = tmpOne - numpy.matmul(numpy.matmul(numpy.matmul(D, KAA), D), alpha) - numpy.matmul(D, tmpOne) * b
          p = numpy.array(list(self.__p(x[i, 0], beta) for i in range(x.shape[0])))
          J2 = numpy.sum(p * p) * c / 2
          return J2
  
  
      def __calc_J1(self, alpha):
          J1 = numpy.sum(alpha * alpha) / 2
          return J1
  
  
      def __get_KAA(self, A, mu):
          KAA = numpy.zeros((A.shape[1], A.shape[1]))
          for rowIdx in range(KAA.shape[0]):
              for colIdx in range(rowIdx + 1):
                  x1 = A[:, rowIdx:rowIdx+1]
                  x2 = A[:, colIdx:colIdx+1]
                  val = self.__calc_gaussian(x1, x2, mu)
                  KAA[rowIdx, colIdx] = KAA[colIdx, rowIdx] = val
          return KAA
  
  
      def __init_alpha_b(self, shape):
          '''
          alpha, b之初始化
          '''
          alpha, b = numpy.zeros(shape), 0
          return alpha, b
  
  
      def __get_KAx(self, A, x, mu):
          KAx = numpy.zeros((A.shape[1], 1))
          for rowIdx in range(KAx.shape[0]):
              x1 = A[:, rowIdx:rowIdx+1]
              val = self.__calc_gaussian(x1, x, mu)
              KAx[rowIdx, 0] = val
          return KAx
  
  
      def __calc_hVal(self, KAx, D, alpha, b):
          hVal = numpy.matmul(numpy.matmul(alpha.T, D), KAx)[0, 0] + b
          return hVal
  
  
      def __calc_gaussian(self, x1, x2, mu):
          val = numpy.math.exp(-mu * numpy.linalg.norm(x1 - x2) ** 2)
          # val = numpy.sum(x1 * x2)
          return val
  
  
      def __get_AD(self):
          A = self.__trainingSet[:, :2].T
          D = numpy.diag(self.__trainingSet[:, 2])
          return A, D
  
  
  class SpiralPlot(object):
  
      def spiral_data_plot(self, trainingData0, trainingData1, testData0, testData1):
          fig = plt.figure(figsize=(5, 5))
          ax1 = plt.subplot()
          ax1.scatter(trainingData1[:, 0], trainingData1[:, 1], c="red", marker="o", s=10, label="training data - Positive")
          ax1.scatter(trainingData0[:, 0], trainingData0[:, 1], c="blue", marker="o", s=10, label="training data - Negative")
          ax1.scatter(testData1[:, 0], testData1[:, 1], c="red", marker="x", s=10, label="test data - Positive")
          ax1.scatter(testData0[:, 0], testData0[:, 1], c="blue", marker="x", s=10, label="test data - Negative")
          ax1.set(xlim=(-0.5, 0.5), ylim=(-0.5, 0.5), xlabel="$x_1$", ylabel="$x_2$")
          plt.legend(fontsize="x-small")
          fig.tight_layout()
          fig.savefig("spiral.png", dpi=100)
          plt.close()
  
  
      def spiral_pred_plot(self, trainingData0, trainingData1, testData0, testData1, ssvmObj, alpha, b):
          x = numpy.linspace(-0.5, 0.5, 500)
          y = numpy.linspace(-0.5, 0.5, 500)
          x, y = numpy.meshgrid(x, y)
          z = numpy.zeros(shape=x.shape)
          for rowIdx in range(x.shape[0]):
              for colIdx in range(x.shape[1]):
                  z[rowIdx, colIdx] = ssvmObj.get_cls((x[rowIdx, colIdx], y[rowIdx, colIdx]), alpha, b)
          cls2color = {-1: "blue", 0: "white", 1: "red"}
  
          fig = plt.figure(figsize=(5, 5))
          ax1 = plt.subplot()
          ax1.contourf(x, y, z, levels=[-1.5, -0.5, 0.5, 1.5], colors=["blue", "white", "red"], alpha=0.3)
          ax1.scatter(trainingData1[:, 0], trainingData1[:, 1], c="red", marker="o", s=10, label="training data - Positive")
          ax1.scatter(trainingData0[:, 0], trainingData0[:, 1], c="blue", marker="o", s=10, label="training data - Negative")
          ax1.scatter(testData1[:, 0], testData1[:, 1], c="red", marker="x", s=10, label="test data - Positive")
          ax1.scatter(testData0[:, 0], testData0[:, 1], c="blue", marker="x", s=10, label="test data - Negative")
          ax1.set(xlim=(-0.5, 0.5), ylim=(-0.5, 0.5), xlabel="$x_1$", ylabel="$x_2$")
          plt.legend(loc="upper left", fontsize="x-small")
          fig.tight_layout()
          fig.savefig("pred.png", dpi=100)
          plt.close()
  
  
  
  if __name__ == "__main__":
      ssvmObj = SSVM(trainingSet, c=0.1, mu=250, beta=100)
      alpha, b, tab = ssvmObj.optimize()
      accuracy1 = ssvmObj.get_accuracy(trainingSet, alpha, b)
      print("Accuracy on trainingSet is {}%".format(accuracy1 * 100))
      accuracy2 = ssvmObj.get_accuracy(testSet, alpha, b)
      print("Accuracy on testSet is {}%".format(accuracy2 * 100))
  
      spObj = SpiralPlot()
      spObj.spiral_data_plot(trainingData0, trainingData1, testData0, testData1)
      spObj.spiral_pred_plot(trainingData0, trainingData1, testData0, testData1, ssvmObj, alpha, b)

  笔者所用训练集、测试集数据分布如下:
  很显然, 此数据集非线性可分.
• 结果展示:
  此模型在训练集、测试集上的准确率均达到100%.
• 使用建议:
  ①. 通过核函数映射至高维后, 支持向量可能存在很多个;
  ②. SSVM本质上求解的仍然是原问题, 此时$\alpha$仅作为变数变换引入, 无需满足KKT条件中对偶可行性;
  ③. 数据集是否需要归一化, 应按实际情况予以考虑.
• 参考文档:
  Yuh-Jye Lee and O. L. Mangasarian. "SSVM: A Smooth Support Vector Machine for Classification", Computational Optimization and Applications, 20, (2001) 5-22.