Smooth Support Vector Clustering

• 算法特征:
  ①所有点尽可能落在球内; ②极小化球半径; ③极小化包络误差.
• 算法推导:
  Part Ⅰ: 模型训练阶段 $Rightarrow$ 形成聚类轮廓
  SVC轮廓线方程如下:
  egin{equation}
  h(x) = |x - a|_2 = sqrt{x^{mathrm{T}}x - 2a^{mathrm{T}}x + a^{mathrm{T}}a}
  label{eq_1}
  end{equation}
  此即样本点$x$距离球心$a$的$L_2$距离.
  本文拟采用$L_2$范数衡量包络误差, 则SVC原始优化问题如下:
  egin{equation}
  egin{split}
  minquad &frac{1}{2}R^2 + frac{c}{2}sum_{i}xi_i^2 \
  mathrm{s.t.}quad &|x^{(i)} - a|_2 leq R + xi_i quadquad forall i = 1, 2, cdots, n
  end{split}
  label{eq_2}
  end{equation}
  其中, $R$是球半径, $a$是球心, $xi_i$是第$i$个样本的包络误差, $c$是权重参数.
  根据上述优化问题( ef{eq_2}), 包络误差$xi_i$可采用plus function等效表示:
  egin{equation}
  xi_i = (|x^{(i)} - a|_2 - R)_+
  label{eq_3}
  end{equation}
  由此可将该有约束最优化问题转化为如下无约束最优化问题:
  egin{equation}
  minquad frac{1}{2}R^2 + frac{c}{2}sum_{i}(|x^{(i)} - a|_2 - R)_+^2
  label{eq_4}
  end{equation}
  同样, 根据优化问题( ef{eq_2})KKT条件, 进行如下变数变换:
  egin{equation}
  a = Aalpha
  label{eq_5}
  end{equation}
  其中, $A = [x^{(1)}, x^{(2)}, cdots, x^{(n)}]$, $alpha = [alpha_1, alpha_2, cdots, alpha_n]^{mathrm{T}}$. $alpha$由原变量及对偶变量组成, 此处无需关心其具体形式. 将上式代入式( ef{eq_4})可得:
  egin{equation}
  minquad frac{1}{2}R^2 + frac{c}{2}sum_{i}(|x^{(i)} - Aalpha|_2 - R)_+^2
  label{eq_6}
  end{equation}
  由于权重参数$c$的存在, 上述最优化问题等效如下多目标最优化问题:
  egin{equation}
  egin{cases}
  minquad frac{1}{2}R^2   \
  minquad frac{1}{2}sumlimits_{i}(|x^{(i)} - Aalpha|_2 - R)_+^2
  end{cases}
  label{eq_7}
  end{equation}
  根据光滑化手段, 有如下等效:
  egin{equation}
  minquad frac{1}{2}sum_{i}(|x^{(i)} - Aalpha|_2 - R)_+^2 quadRightarrowquad minquad frac{1}{2}sum_i p^2(|x^{(i)} - Aalpha|_2 - R, eta), quad ext{where }eta ightarrowinfty
  label{eq_8}
  end{equation}
  光滑化手段详见:
  https://www.cnblogs.com/xxhbdk/p/12275567.html
  https://www.cnblogs.com/xxhbdk/p/12425701.html
  结合式( ef{eq_7})、式( ef{eq_8}), 优化问题( ef{eq_6})等效如下:
  egin{equation}
  minquad frac{1}{2}R^2 + frac{c}{2}sum_i p^2(sqrt{x^{(i)mathrm{T}}x^{(i)} - 2x^{(i)mathrm{T}}Aalpha + alpha^{mathrm{T}}A^{mathrm{T}}Aalpha} - R, eta), quad ext{where }eta ightarrowinfty
  label{eq_9}
  end{equation}
  为增强模型的聚类能力, 引入kernel trick, 得最终优化问题:
  egin{equation}
  minquad frac{1}{2}R^2 + frac{c}{2}sum_i p^2(sqrt{K(x^{(i)mathrm{T}}, x^{(i)}) - 2K(x^{(i)mathrm{T}}, A)alpha + alpha^{mathrm{T}}K(A^{mathrm{T}}, A)alpha} - R, eta), quad ext{where }eta ightarrowinfty
  label{eq_10}
  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}$
  结合式( ef{eq_1})、式( ef{eq_5})及kernel trick可得最终轮廓线方程:
  egin{equation}
  h(x) = sqrt{K(x^{mathrm{T}}, x) - 2alpha^{mathrm{T}}K(A^{mathrm{T}}, x) + alpha^{mathrm{T}}K(A^{mathrm{T}}, A)alpha}
  label{eq_11}
  end{equation}
  
  Part Ⅱ: 聚类分配阶段 $Rightarrow$ 簇划分
  根据如下规则构造邻接矩阵$G$:
  egin{equation}
  G_{ij} =
  egin{cases}
  1quad ext{if $h(x) leq R$, where $x = lambda x^{(i)} + (1 - lambda)x^{(j)}$, $forall lambda in (0, 1)$}   \
  0quad ext{otherwise}
  end{cases}
  end{equation}
  对该矩阵所定义的图进行深度优先或广度优先遍历, 所得连通分量即为簇的划分.
  
  Part Ⅲ: 优化协助信息
  拟设式( ef{eq_10})之目标函数为$J$, 则有:
  egin{equation}
  J = frac{1}{2}R^2 + frac{c}{2}sum_i p^2(z_i - R, eta)
  end{equation}
  其中, $z_i = sqrt{K(x^{(i)mathrm{T}}, x^{(i)}) - 2K(x^{(i)mathrm{T}}, A)alpha + alpha^{mathrm{T}}K(A^{mathrm{T}}, A)alpha}$.
  相关一阶信息如下:
  egin{gather*}
  y_i = K(A^{mathrm{T}}, A)alpha - K^{mathrm{T}}(x^{(i)mathrm{T}}, A)   \
  frac{partial J}{partial alpha} = csum_ip(z_i - R, eta)s(z_i - R, eta)z_i^{-1}y_i   \
  frac{partial J}{partial R} = R - csum_ip(z_i - R, eta)s(z_i - R, eta)
  end{gather*}
  目标函数$J$之gradient如下:
  egin{gather*}
  abla J =
  egin{bmatrix}
  frac{partial J}{partial alpha} \
  frac{partial J}{partial R}
  end{bmatrix}
  end{gather*}
  
  Part Ⅳ: 聚类数据集
  现以如下环状数据集为例进行算法实施, 相关数据分布如下:
  

  # 生成聚类数据集 - 环状
  
  import numpy
  from matplotlib import pyplot as plt
  
  
  numpy.random.seed(0)
  
  
  def ring(type1Size=100, type2Size=100):
      theta1 = numpy.random.uniform(0, 2 * numpy.pi, type1Size)
      theta2 = numpy.random.uniform(0, 2 * numpy.pi, type2Size)
      
      R1 = numpy.random.uniform(0, 0.3, type1Size).reshape((-1, 1))
      R2 = numpy.random.uniform(0.7, 1, type2Size).reshape((-1, 1))
      X1 = R1 * numpy.hstack((numpy.cos(theta1).reshape((-1, 1)), numpy.sin(theta1).reshape((-1, 1))))
      X2 = R2 * numpy.hstack((numpy.cos(theta2).reshape((-1, 1)), numpy.sin(theta2).reshape((-1, 1))))
      
      X = numpy.vstack((X1, X2))
      return X
  
  
  def ring_plot():
      X = ring(300, 300)
      
      fig = plt.figure(figsize=(5, 3))
      ax1 = plt.subplot()
      
      ax1.scatter(X[:, 0], X[:, 1], marker="o", color="red", s=5)
      ax1.set(xlim=[-1.1, 1.1], ylim=[-1.1, 1.1], xlabel="$x_1$", ylabel="$x_2$")
      
      fig.tight_layout()
      fig.savefig("ring.png", dpi=100)
      
      
      
  if __name__ == "__main__":
      ring_plot()

  

• 代码实现:
  

  # Smooth Support Vector Clustering之实现
  
  import copy
  import numpy
  from matplotlib import pyplot as plt
  
  numpy.random.seed(0)
  
  
  def ring(type1Size=100, type2Size=100):
      theta1 = numpy.random.uniform(0, 2 * numpy.pi, type1Size)
      theta2 = numpy.random.uniform(0, 2 * numpy.pi, type2Size)
      
      R1 = numpy.random.uniform(0, 0.3, type1Size).reshape((-1, 1))
      R2 = numpy.random.uniform(0.7, 1, type2Size).reshape((-1, 1))
      X1 = R1 * numpy.hstack((numpy.cos(theta1).reshape((-1, 1)), numpy.sin(theta1).reshape((-1, 1))))
      X2 = R2 * numpy.hstack((numpy.cos(theta2).reshape((-1, 1)), numpy.sin(theta2).reshape((-1, 1))))
      
      # R1 = numpy.random.uniform(0, 0.3, type1Size).reshape((-1, 1))
      # X1 = R1 * numpy.hstack((numpy.cos(theta1).reshape((-1, 1)), numpy.sin(theta1).reshape((-1, 1)))) - 0.5
      # X2 = R1 * numpy.hstack((numpy.cos(theta2).reshape((-1, 1)), numpy.sin(theta2).reshape((-1, 1)))) + 0.5
      
      X = numpy.vstack((X1, X2))
      return X
      
      
  # 生成聚类数据
  X = ring(300, 300)
  
  
  class SSVC(object):
      
      def __init__(self, X, c=1, mu=1, beta=300, sampleCnt=30):
          '''
          X: 聚类数据集, 1行代表一个样本
          '''
          self.__X = X                                 # 聚类数据集
          self.__c = c                                 # 误差项权重系数
          self.__mu = mu                               # gaussian kernel参数
          self.__beta = beta                           # 光滑华参数
          self.__sampleCnt = sampleCnt                 # 邻接矩阵采样数
          
          self.__A = X.T
          
      
      def get_clusters(self, IDPs, SVs, alpha, R, KAA=None):
          '''
          构造邻接矩阵, 进行簇划分
          注意, 此处仅需关注IDPs、SVs
          '''
          if KAA is None:
              KAA = self.get_KAA()
          
          if IDPs.shape[0] == 0 and SVs.shape[0] == 0:
              raise Exception(">>> Both IDPs and SVs are empty! <<<")
          elif IDPs.shape[0] == 0 and SVs.shape[0] != 0:
              currX = SVs
          elif IDPs.shape[0] != 0 and SVs.shape[0] == 0:
              currX = IDPs
          else:
              currX = numpy.vstack((IDPs, SVs))
          adjMat = self.__build_adjMat(currX, alpha, R, KAA)
          
          idxList = list(range(currX.shape[0]))
          clusters = self.__get_clusters(idxList, adjMat)
          
          clustersOri = list()
          for idx, cluster in enumerate(clusters):
              print('Size of cluster {} is about {}'.format(idx, len(cluster)))
              clusterOri = list(currX[idx, :] for idx in cluster)
              clustersOri.append(numpy.array(clusterOri))
          return clustersOri
          
      
      def get_IDPs_SVs_BSVs(self, alpha, R):
          '''
          获取IDPs: Interior Data Points
          获取SVs: Support Vectors
          获取BSVs: Bounded Support Vectors
          '''
          X = self.__X
          KAA = self.get_KAA()
          
          epsilon = 3.e-3 * R
          IDPs, SVs, BSVs = list(), list(), list()
          
          for x in X:
              hVal = self.get_hVal(x, alpha, KAA)
              if hVal >= R + epsilon:
                  BSVs.append(x)
              elif numpy.abs(hVal - R) < epsilon:
                  SVs.append(x)
              else:
                  IDPs.append(x)
          return numpy.array(IDPs), numpy.array(SVs), numpy.array(BSVs)
          
          
      def get_hVal(self, x, alpha, KAA=None):
          '''
          计算x与超球球心在特征空间中的距离
          '''
          A = self.__A
          mu = self.__mu
          
          x = numpy.array(x).reshape((-1, 1))
          KAx = self.__get_KAx(A, x, mu)
          if KAA is None:
              KAA = self.__get_KAA(A, mu)
          term1 = self.__calc_gaussian(x, x, mu)
          term2 = -2 * numpy.matmul(alpha.T, KAx)[0, 0]
          term3 = numpy.matmul(numpy.matmul(alpha.T, KAA), alpha)[0, 0]
          hVal = numpy.sqrt(term1 + term2 + term3)
          return hVal
          
          
      def get_KAA(self):
          A = self.__A
          mu = self.__mu
          
          KAA = self.__get_KAA(A, mu)
          return KAA
      
      
      def optimize(self, maxIter=3000, epsilon=1.e-9):
          '''
          maxIter: 最大迭代次数
          epsilon: 收敛判据, 梯度趋于0则收敛
          '''
          A = self.__A
          c = self.__c
          mu = self.__mu
          beta = self.__beta
          
          alpha, R = self.__init_alpha_R((A.shape[1], 1))
          KAA = self.__get_KAA(A, mu)
          JVal = self.__calc_JVal(KAA, c, beta, alpha, R)
          grad = self.__calc_grad(KAA, c, beta, alpha, R)
          D = self.__init_D(KAA.shape[0] + 1)
          
          for i in range(maxIter):
              print("iterCnt: {:3d},   JVal: {}".format(i, JVal))
              if self.__converged1(grad, epsilon):
                  return alpha, R, True
              
              dCurr = -numpy.matmul(D, grad)
              ALPHA = self.__calc_ALPHA_by_ArmijoRule(alpha, R, JVal, grad, dCurr, KAA, c, beta)
              
              delta = ALPHA * dCurr
              alphaNew = alpha + delta[:-1, :]
              RNew = R + delta[-1, -1]
              JValNew = self.__calc_JVal(KAA, c, beta, alphaNew, RNew)
              if self.__converged2(delta, JValNew - JVal, epsilon ** 2):
                  return alpha, R, True
              
              gradNew = self.__calc_grad(KAA, c, beta, alphaNew, RNew)
              DNew = self.__update_D_by_BFGS(delta, gradNew - grad, D)
              
              alpha, R, JVal, grad, D = alphaNew, RNew, JValNew, gradNew, DNew
          else:
              if self.__converged1(grad, epsilon):
                  return alpha, R, True
          return alpha, R, False
          
      
      def __update_D_by_BFGS(self, sk, yk, D):
          rk = 1 / numpy.matmul(yk.T, sk)[0, 0]
          
          term1 = rk * numpy.matmul(sk, yk.T)
          term2 = rk * numpy.matmul(yk, sk.T)
          I = numpy.identity(term1.shape[0])
          term3 = numpy.matmul(I - term1, D)
          term4 = numpy.matmul(term3, I - term2)
          term5 = rk * numpy.matmul(sk, sk.T)
          
          DNew = term4 + term5
          return DNew
      
      
      def __calc_ALPHA_by_ArmijoRule(self, alphaCurr, RCurr, JCurr, gCurr, dCurr, KAA, c, beta, C=1.e-4, v=0.5):
          i = 0
          ALPHA = v ** i
          delta = ALPHA * dCurr
          alphaNext = alphaCurr + delta[:-1, :]
          RNext = RCurr + delta[-1, -1]
          JNext = self.__calc_JVal(KAA, c, beta, alphaNext, RNext)
          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, :]
              RNext = RCurr + delta[-1, -1]
              JNext = self.__calc_JVal(KAA, c, beta, alphaNext, RNext)
          return ALPHA
      
      
      def __converged1(self, grad, epsilon):
          if numpy.linalg.norm(grad, ord=numpy.inf) <= epsilon:
              return True
          return 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 __init_D(self, n):
          D = numpy.identity(n)
          return D
      
      
      def __init_alpha_R(self, shape):
          '''
          alpha、R之初始化
          '''
          alpha, R = numpy.zeros(shape), 0
          return alpha, R
      
          
      def __calc_grad(self, KAA, c, beta, alpha, R):
          grad_alpha = numpy.zeros((KAA.shape[0], 1))
          grad_R = 0
          
          term = numpy.matmul(numpy.matmul(alpha.T, KAA), alpha)[0, 0]
          z = list(numpy.sqrt(KAA[idx, idx] - 2 * numpy.matmul(KAA[idx:idx+1, :], alpha)[0, 0] + term) for idx in range(KAA.shape[0]))
          term1 = numpy.matmul(KAA, alpha)
          y = list(term1 - KAA[:, idx:idx+1] for idx in range(KAA.shape[0]))
          
          for i in range(KAA.shape[0]):
              term2 = self.__p(z[i] - R, beta) * self.__s(z[i] - R, beta)
              grad_alpha += term2 * y[i] / z[i]
              grad_R += term2
          grad_alpha *= c
          grad_R = R - c * grad_R
          
          grad = numpy.vstack((grad_alpha, [[grad_R]]))
          return grad
          
          
      def __calc_JVal(self, KAA, c, beta, alpha, R):
          term = numpy.matmul(numpy.matmul(alpha.T, KAA), alpha)[0, 0]
          z = list(numpy.sqrt(KAA[idx, idx] - 2 * numpy.matmul(KAA[idx:idx+1, :], alpha)[0, 0] + term) for idx in range(KAA.shape[0]))
          
          term1 = R ** 2 / 2
          term2 = sum(self.__p(zi - R, beta) ** 2 for zi in z) * c / 2
          JVal = term1 + term2
          return JVal
          
          
      def __p(self, x, beta):
          term = x * beta
          if term > 10:
              val = x + numpy.log(1 + numpy.exp(-term)) / beta
          else:
              val = numpy.log(numpy.exp(term) + 1) / beta
          return val
          
          
      def __s(self, x, beta):
          term = x * beta
          if term > 10:
              val = 1 / (numpy.exp(-beta * x) + 1)
          else:
              term1 = numpy.exp(term)
              val = term1 / (1 + term1)
          return val
          
          
      def __d(self, x, beta):
          term = x * beta
          if term > 10:
              term1 = numpy.exp(-term)
              val = beta * term1 / (1 + term1) ** 2
          else:
              term1 = numpy.exp(term)
              val = beta * term1 / (term1 + 1) ** 2
          return val
          
          
      def __calc_gaussian(self, x1, x2, mu):
          val = numpy.exp(-mu * numpy.linalg.norm(x1 - x2) ** 2)
          # val = numpy.sum(x1 * x2)
          return val
          
          
      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 __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 __get_segment(self, xi, xj, sampleCnt):
          lamdaList = numpy.linspace(0, 1, sampleCnt + 2)
          segmentList = list(lamda * xi + (1 - lamda) * xj for lamda in lamdaList[1:-1])
          return segmentList
          
          
      def __build_adjMat(self, X, alpha, R, KAA):
          sampleCnt = self.__sampleCnt
          
          adjMat = numpy.zeros((X.shape[0], X.shape[0]))
          for i in range(adjMat.shape[0]):
              for j in range(i+1, adjMat.shape[1]):
                  xi = X[i, :]
                  xj = X[j, :]
                  xs = self.__get_segment(xi, xj, sampleCnt)
                  for x in xs:
                      dist = self.get_hVal(x, alpha, KAA)
                      if dist > R:
                          adjMat[i, j] = adjMat[j, i] = 0
                          break
                  else:
                      adjMat[i, j] = adjMat[j, i] = 1
  
          numpy.savetxt("adjMat.txt", adjMat, fmt="%d")
          return adjMat
          
          
      def __get_clusters(self, idxList, adjMat):
          clusters = list()
          while idxList:
              idxCurr = idxList.pop(0)
              queue = [idxCurr]
              cluster = list()
              while queue:
                  idxPop = queue.pop(0)
                  cluster.append(idxPop)
                  for idx in copy.deepcopy(idxList):
                      if adjMat[idxPop, idx] == 1:
                          idxList.remove(idx)
                          queue.append(idx)
              clusters.append(cluster)
          
          return clusters
              
              
  class SSVCPlot(object):
      
      def profile_plot(self, X, ssvcObj, alpha, R):
          surfX1 = numpy.linspace(-1.1, 1.1, 100)
          surfX2 = numpy.linspace(-1.1, 1.1, 100)
          surfX1, surfX2 = numpy.meshgrid(surfX1, surfX2)
          
          KAA = ssvcObj.get_KAA()
          surfY = numpy.zeros(surfX1.shape)
          for rowIdx in range(surfY.shape[0]):
              for colIdx in range(surfY.shape[1]):
                  surfY[rowIdx, colIdx] = ssvcObj.get_hVal((surfX1[rowIdx, colIdx], surfX2[rowIdx, colIdx]), alpha, KAA)
          
          IDPs, SVs, BSVs = ssvcObj.get_IDPs_SVs_BSVs(alpha, R)
          clusters = ssvcObj.get_clusters(IDPs, SVs, alpha, R, KAA)
          
          fig = plt.figure(figsize=(10, 3))
          ax1 = plt.subplot(1, 2, 1)
          ax2 = plt.subplot(1, 2, 2)
          
          ax1.contour(surfX1, surfX2, surfY, levels=36)
          ax1.scatter(X[:, 0], X[:, 1], marker="o", color="red", s=5)
          ax1.set(xlim=[-1.1, 1.1], ylim=[-1.1, 1.1], xlabel="$x_1$", ylabel="$x_2$")
          
          if BSVs.shape[0] != 0:
              ax2.scatter(BSVs[:, 0], BSVs[:, 1], color="k", s=5, label="$BSVs$")
          for idx, cluster in enumerate(clusters):
              ax2.scatter(cluster[:, 0], cluster[:, 1], s=3, label="$cluster {}$".format(idx))
          ax2.set(xlim=[-1.1, 1.1], ylim=[-1.1, 1.1], xlabel="$x_1$", ylabel="$x_2$")
          ax2.legend()
          
          fig.tight_layout()
          fig.savefig("profile.png", dpi=100)
          
  
  
  if __name__ == "__main__":
      ssvcObj = SSVC(X, c=100, mu=7, beta=300)
      alpha, R, tab = ssvcObj.optimize()
      spObj = SSVCPlot()
      spObj.profile_plot(X, ssvcObj, alpha, R)
      

• 结果展示:
  左侧为聚类轮廓, 右侧为簇划分. 很显然, 此处在合适的超参数选取下, 将原始数据集划分为2个簇.
• 使用建议:
  ①. 本文的优化问题与参考文档的优化问题略有差异, 参考文档中关于SVC的原问题并非严格意义上的凸问题.
  ②. SSVC本质上求解的仍然是原问题, 此时$alpha$仅作为变数变换引入, 无需满足KKT条件中对偶可行性.
  ③. 数据集是否需要归一化, 应按实际情况予以考虑; 簇的划分敏感于超参数的选取.
• 参考文档:
  Ben-Hur A, Horn D, Siegelmann HT, Vapnik V. Support vector clustering. Journal of Machine Learning Research, 2001, 2(12): 125-137.