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 Ⅳ: 聚类数据集
    现以如下环状数据集为例进行算法实施, 相关数据分布如下:
     1 # 生成聚类数据集 - 环状
     2 
     3 import numpy
     4 from matplotlib import pyplot as plt
     5 
     6 
     7 numpy.random.seed(0)
     8 
     9 
    10 def ring(type1Size=100, type2Size=100):
    11     theta1 = numpy.random.uniform(0, 2 * numpy.pi, type1Size)
    12     theta2 = numpy.random.uniform(0, 2 * numpy.pi, type2Size)
    13     
    14     R1 = numpy.random.uniform(0, 0.3, type1Size).reshape((-1, 1))
    15     R2 = numpy.random.uniform(0.7, 1, type2Size).reshape((-1, 1))
    16     X1 = R1 * numpy.hstack((numpy.cos(theta1).reshape((-1, 1)), numpy.sin(theta1).reshape((-1, 1))))
    17     X2 = R2 * numpy.hstack((numpy.cos(theta2).reshape((-1, 1)), numpy.sin(theta2).reshape((-1, 1))))
    18     
    19     X = numpy.vstack((X1, X2))
    20     return X
    21 
    22 
    23 def ring_plot():
    24     X = ring(300, 300)
    25     
    26     fig = plt.figure(figsize=(5, 3))
    27     ax1 = plt.subplot()
    28     
    29     ax1.scatter(X[:, 0], X[:, 1], marker="o", color="red", s=5)
    30     ax1.set(xlim=[-1.1, 1.1], ylim=[-1.1, 1.1], xlabel="$x_1$", ylabel="$x_2$")
    31     
    32     fig.tight_layout()
    33     fig.savefig("ring.png", dpi=100)
    34     
    35     
    36     
    37 if __name__ == "__main__":
    38     ring_plot()
    View Code

  • 代码实现:
      1 # Smooth Support Vector Clustering之实现
      2 
      3 import copy
      4 import numpy
      5 from matplotlib import pyplot as plt
      6 
      7 numpy.random.seed(0)
      8 
      9 
     10 def ring(type1Size=100, type2Size=100):
     11     theta1 = numpy.random.uniform(0, 2 * numpy.pi, type1Size)
     12     theta2 = numpy.random.uniform(0, 2 * numpy.pi, type2Size)
     13     
     14     R1 = numpy.random.uniform(0, 0.3, type1Size).reshape((-1, 1))
     15     R2 = numpy.random.uniform(0.7, 1, type2Size).reshape((-1, 1))
     16     X1 = R1 * numpy.hstack((numpy.cos(theta1).reshape((-1, 1)), numpy.sin(theta1).reshape((-1, 1))))
     17     X2 = R2 * numpy.hstack((numpy.cos(theta2).reshape((-1, 1)), numpy.sin(theta2).reshape((-1, 1))))
     18     
     19     # R1 = numpy.random.uniform(0, 0.3, type1Size).reshape((-1, 1))
     20     # X1 = R1 * numpy.hstack((numpy.cos(theta1).reshape((-1, 1)), numpy.sin(theta1).reshape((-1, 1)))) - 0.5
     21     # X2 = R1 * numpy.hstack((numpy.cos(theta2).reshape((-1, 1)), numpy.sin(theta2).reshape((-1, 1)))) + 0.5
     22     
     23     X = numpy.vstack((X1, X2))
     24     return X
     25     
     26     
     27 # 生成聚类数据
     28 X = ring(300, 300)
     29 
     30 
     31 class SSVC(object):
     32     
     33     def __init__(self, X, c=1, mu=1, beta=300, sampleCnt=30):
     34         '''
     35         X: 聚类数据集, 1行代表一个样本
     36         '''
     37         self.__X = X                                 # 聚类数据集
     38         self.__c = c                                 # 误差项权重系数
     39         self.__mu = mu                               # gaussian kernel参数
     40         self.__beta = beta                           # 光滑华参数
     41         self.__sampleCnt = sampleCnt                 # 邻接矩阵采样数
     42         
     43         self.__A = X.T
     44         
     45     
     46     def get_clusters(self, IDPs, SVs, alpha, R, KAA=None):
     47         '''
     48         构造邻接矩阵, 进行簇划分
     49         注意, 此处仅需关注IDPs、SVs
     50         '''
     51         if KAA is None:
     52             KAA = self.get_KAA()
     53         
     54         if IDPs.shape[0] == 0 and SVs.shape[0] == 0:
     55             raise Exception(">>> Both IDPs and SVs are empty! <<<")
     56         elif IDPs.shape[0] == 0 and SVs.shape[0] != 0:
     57             currX = SVs
     58         elif IDPs.shape[0] != 0 and SVs.shape[0] == 0:
     59             currX = IDPs
     60         else:
     61             currX = numpy.vstack((IDPs, SVs))
     62         adjMat = self.__build_adjMat(currX, alpha, R, KAA)
     63         
     64         idxList = list(range(currX.shape[0]))
     65         clusters = self.__get_clusters(idxList, adjMat)
     66         
     67         clustersOri = list()
     68         for idx, cluster in enumerate(clusters):
     69             print('Size of cluster {} is about {}'.format(idx, len(cluster)))
     70             clusterOri = list(currX[idx, :] for idx in cluster)
     71             clustersOri.append(numpy.array(clusterOri))
     72         return clustersOri
     73         
     74     
     75     def get_IDPs_SVs_BSVs(self, alpha, R):
     76         '''
     77         获取IDPs: Interior Data Points
     78         获取SVs: Support Vectors
     79         获取BSVs: Bounded Support Vectors
     80         '''
     81         X = self.__X
     82         KAA = self.get_KAA()
     83         
     84         epsilon = 3.e-3 * R
     85         IDPs, SVs, BSVs = list(), list(), list()
     86         
     87         for x in X:
     88             hVal = self.get_hVal(x, alpha, KAA)
     89             if hVal >= R + epsilon:
     90                 BSVs.append(x)
     91             elif numpy.abs(hVal - R) < epsilon:
     92                 SVs.append(x)
     93             else:
     94                 IDPs.append(x)
     95         return numpy.array(IDPs), numpy.array(SVs), numpy.array(BSVs)
     96         
     97         
     98     def get_hVal(self, x, alpha, KAA=None):
     99         '''
    100         计算x与超球球心在特征空间中的距离
    101         '''
    102         A = self.__A
    103         mu = self.__mu
    104         
    105         x = numpy.array(x).reshape((-1, 1))
    106         KAx = self.__get_KAx(A, x, mu)
    107         if KAA is None:
    108             KAA = self.__get_KAA(A, mu)
    109         term1 = self.__calc_gaussian(x, x, mu)
    110         term2 = -2 * numpy.matmul(alpha.T, KAx)[0, 0]
    111         term3 = numpy.matmul(numpy.matmul(alpha.T, KAA), alpha)[0, 0]
    112         hVal = numpy.sqrt(term1 + term2 + term3)
    113         return hVal
    114         
    115         
    116     def get_KAA(self):
    117         A = self.__A
    118         mu = self.__mu
    119         
    120         KAA = self.__get_KAA(A, mu)
    121         return KAA
    122     
    123     
    124     def optimize(self, maxIter=3000, epsilon=1.e-9):
    125         '''
    126         maxIter: 最大迭代次数
    127         epsilon: 收敛判据, 梯度趋于0则收敛
    128         '''
    129         A = self.__A
    130         c = self.__c
    131         mu = self.__mu
    132         beta = self.__beta
    133         
    134         alpha, R = self.__init_alpha_R((A.shape[1], 1))
    135         KAA = self.__get_KAA(A, mu)
    136         JVal = self.__calc_JVal(KAA, c, beta, alpha, R)
    137         grad = self.__calc_grad(KAA, c, beta, alpha, R)
    138         D = self.__init_D(KAA.shape[0] + 1)
    139         
    140         for i in range(maxIter):
    141             print("iterCnt: {:3d},   JVal: {}".format(i, JVal))
    142             if self.__converged1(grad, epsilon):
    143                 return alpha, R, True
    144             
    145             dCurr = -numpy.matmul(D, grad)
    146             ALPHA = self.__calc_ALPHA_by_ArmijoRule(alpha, R, JVal, grad, dCurr, KAA, c, beta)
    147             
    148             delta = ALPHA * dCurr
    149             alphaNew = alpha + delta[:-1, :]
    150             RNew = R + delta[-1, -1]
    151             JValNew = self.__calc_JVal(KAA, c, beta, alphaNew, RNew)
    152             if self.__converged2(delta, JValNew - JVal, epsilon ** 2):
    153                 return alpha, R, True
    154             
    155             gradNew = self.__calc_grad(KAA, c, beta, alphaNew, RNew)
    156             DNew = self.__update_D_by_BFGS(delta, gradNew - grad, D)
    157             
    158             alpha, R, JVal, grad, D = alphaNew, RNew, JValNew, gradNew, DNew
    159         else:
    160             if self.__converged1(grad, epsilon):
    161                 return alpha, R, True
    162         return alpha, R, False
    163         
    164     
    165     def __update_D_by_BFGS(self, sk, yk, D):
    166         rk = 1 / numpy.matmul(yk.T, sk)[0, 0]
    167         
    168         term1 = rk * numpy.matmul(sk, yk.T)
    169         term2 = rk * numpy.matmul(yk, sk.T)
    170         I = numpy.identity(term1.shape[0])
    171         term3 = numpy.matmul(I - term1, D)
    172         term4 = numpy.matmul(term3, I - term2)
    173         term5 = rk * numpy.matmul(sk, sk.T)
    174         
    175         DNew = term4 + term5
    176         return DNew
    177     
    178     
    179     def __calc_ALPHA_by_ArmijoRule(self, alphaCurr, RCurr, JCurr, gCurr, dCurr, KAA, c, beta, C=1.e-4, v=0.5):
    180         i = 0
    181         ALPHA = v ** i
    182         delta = ALPHA * dCurr
    183         alphaNext = alphaCurr + delta[:-1, :]
    184         RNext = RCurr + delta[-1, -1]
    185         JNext = self.__calc_JVal(KAA, c, beta, alphaNext, RNext)
    186         while True:
    187             if JNext <= JCurr + C * ALPHA * numpy.matmul(dCurr.T, gCurr)[0, 0]: break
    188             i += 1
    189             ALPHA = v ** i
    190             delta = ALPHA * dCurr
    191             alphaNext = alphaCurr + delta[:-1, :]
    192             RNext = RCurr + delta[-1, -1]
    193             JNext = self.__calc_JVal(KAA, c, beta, alphaNext, RNext)
    194         return ALPHA
    195     
    196     
    197     def __converged1(self, grad, epsilon):
    198         if numpy.linalg.norm(grad, ord=numpy.inf) <= epsilon:
    199             return True
    200         return False
    201         
    202         
    203     def __converged2(self, delta, JValDelta, epsilon):
    204         val1 = numpy.linalg.norm(delta, ord=numpy.inf)
    205         val2 = numpy.abs(JValDelta)
    206         if val1 <= epsilon or val2 <= epsilon:
    207             return True
    208         return False
    209     
    210 
    211     def __init_D(self, n):
    212         D = numpy.identity(n)
    213         return D
    214     
    215     
    216     def __init_alpha_R(self, shape):
    217         '''
    218         alpha、R之初始化
    219         '''
    220         alpha, R = numpy.zeros(shape), 0
    221         return alpha, R
    222     
    223         
    224     def __calc_grad(self, KAA, c, beta, alpha, R):
    225         grad_alpha = numpy.zeros((KAA.shape[0], 1))
    226         grad_R = 0
    227         
    228         term = numpy.matmul(numpy.matmul(alpha.T, KAA), alpha)[0, 0]
    229         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]))
    230         term1 = numpy.matmul(KAA, alpha)
    231         y = list(term1 - KAA[:, idx:idx+1] for idx in range(KAA.shape[0]))
    232         
    233         for i in range(KAA.shape[0]):
    234             term2 = self.__p(z[i] - R, beta) * self.__s(z[i] - R, beta)
    235             grad_alpha += term2 * y[i] / z[i]
    236             grad_R += term2
    237         grad_alpha *= c
    238         grad_R = R - c * grad_R
    239         
    240         grad = numpy.vstack((grad_alpha, [[grad_R]]))
    241         return grad
    242         
    243         
    244     def __calc_JVal(self, KAA, c, beta, alpha, R):
    245         term = numpy.matmul(numpy.matmul(alpha.T, KAA), alpha)[0, 0]
    246         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]))
    247         
    248         term1 = R ** 2 / 2
    249         term2 = sum(self.__p(zi - R, beta) ** 2 for zi in z) * c / 2
    250         JVal = term1 + term2
    251         return JVal
    252         
    253         
    254     def __p(self, x, beta):
    255         term = x * beta
    256         if term > 10:
    257             val = x + numpy.log(1 + numpy.exp(-term)) / beta
    258         else:
    259             val = numpy.log(numpy.exp(term) + 1) / beta
    260         return val
    261         
    262         
    263     def __s(self, x, beta):
    264         term = x * beta
    265         if term > 10:
    266             val = 1 / (numpy.exp(-beta * x) + 1)
    267         else:
    268             term1 = numpy.exp(term)
    269             val = term1 / (1 + term1)
    270         return val
    271         
    272         
    273     def __d(self, x, beta):
    274         term = x * beta
    275         if term > 10:
    276             term1 = numpy.exp(-term)
    277             val = beta * term1 / (1 + term1) ** 2
    278         else:
    279             term1 = numpy.exp(term)
    280             val = beta * term1 / (term1 + 1) ** 2
    281         return val
    282         
    283         
    284     def __calc_gaussian(self, x1, x2, mu):
    285         val = numpy.exp(-mu * numpy.linalg.norm(x1 - x2) ** 2)
    286         # val = numpy.sum(x1 * x2)
    287         return val
    288         
    289         
    290     def __get_KAA(self, A, mu):
    291         KAA = numpy.zeros((A.shape[1], A.shape[1]))
    292         for rowIdx in range(KAA.shape[0]):
    293             for colIdx in range(rowIdx + 1):
    294                 x1 = A[:, rowIdx:rowIdx+1]
    295                 x2 = A[:, colIdx:colIdx+1]
    296                 val = self.__calc_gaussian(x1, x2, mu)
    297                 KAA[rowIdx, colIdx] = KAA[colIdx, rowIdx] = val
    298         return KAA
    299         
    300         
    301     def __get_KAx(self, A, x, mu):
    302         KAx = numpy.zeros((A.shape[1], 1))
    303         for rowIdx in range(KAx.shape[0]):
    304             x1 = A[:, rowIdx:rowIdx+1]
    305             val = self.__calc_gaussian(x1, x, mu)
    306             KAx[rowIdx, 0] = val
    307         return KAx
    308         
    309         
    310     def __get_segment(self, xi, xj, sampleCnt):
    311         lamdaList = numpy.linspace(0, 1, sampleCnt + 2)
    312         segmentList = list(lamda * xi + (1 - lamda) * xj for lamda in lamdaList[1:-1])
    313         return segmentList
    314         
    315         
    316     def __build_adjMat(self, X, alpha, R, KAA):
    317         sampleCnt = self.__sampleCnt
    318         
    319         adjMat = numpy.zeros((X.shape[0], X.shape[0]))
    320         for i in range(adjMat.shape[0]):
    321             for j in range(i+1, adjMat.shape[1]):
    322                 xi = X[i, :]
    323                 xj = X[j, :]
    324                 xs = self.__get_segment(xi, xj, sampleCnt)
    325                 for x in xs:
    326                     dist = self.get_hVal(x, alpha, KAA)
    327                     if dist > R:
    328                         adjMat[i, j] = adjMat[j, i] = 0
    329                         break
    330                 else:
    331                     adjMat[i, j] = adjMat[j, i] = 1
    332 
    333         numpy.savetxt("adjMat.txt", adjMat, fmt="%d")
    334         return adjMat
    335         
    336         
    337     def __get_clusters(self, idxList, adjMat):
    338         clusters = list()
    339         while idxList:
    340             idxCurr = idxList.pop(0)
    341             queue = [idxCurr]
    342             cluster = list()
    343             while queue:
    344                 idxPop = queue.pop(0)
    345                 cluster.append(idxPop)
    346                 for idx in copy.deepcopy(idxList):
    347                     if adjMat[idxPop, idx] == 1:
    348                         idxList.remove(idx)
    349                         queue.append(idx)
    350             clusters.append(cluster)
    351         
    352         return clusters
    353             
    354             
    355 class SSVCPlot(object):
    356     
    357     def profile_plot(self, X, ssvcObj, alpha, R):
    358         surfX1 = numpy.linspace(-1.1, 1.1, 100)
    359         surfX2 = numpy.linspace(-1.1, 1.1, 100)
    360         surfX1, surfX2 = numpy.meshgrid(surfX1, surfX2)
    361         
    362         KAA = ssvcObj.get_KAA()
    363         surfY = numpy.zeros(surfX1.shape)
    364         for rowIdx in range(surfY.shape[0]):
    365             for colIdx in range(surfY.shape[1]):
    366                 surfY[rowIdx, colIdx] = ssvcObj.get_hVal((surfX1[rowIdx, colIdx], surfX2[rowIdx, colIdx]), alpha, KAA)
    367         
    368         IDPs, SVs, BSVs = ssvcObj.get_IDPs_SVs_BSVs(alpha, R)
    369         clusters = ssvcObj.get_clusters(IDPs, SVs, alpha, R, KAA)
    370         
    371         fig = plt.figure(figsize=(10, 3))
    372         ax1 = plt.subplot(1, 2, 1)
    373         ax2 = plt.subplot(1, 2, 2)
    374         
    375         ax1.contour(surfX1, surfX2, surfY, levels=36)
    376         ax1.scatter(X[:, 0], X[:, 1], marker="o", color="red", s=5)
    377         ax1.set(xlim=[-1.1, 1.1], ylim=[-1.1, 1.1], xlabel="$x_1$", ylabel="$x_2$")
    378         
    379         if BSVs.shape[0] != 0:
    380             ax2.scatter(BSVs[:, 0], BSVs[:, 1], color="k", s=5, label="$BSVs$")
    381         for idx, cluster in enumerate(clusters):
    382             ax2.scatter(cluster[:, 0], cluster[:, 1], s=3, label="$cluster {}$".format(idx))
    383         ax2.set(xlim=[-1.1, 1.1], ylim=[-1.1, 1.1], xlabel="$x_1$", ylabel="$x_2$")
    384         ax2.legend()
    385         
    386         fig.tight_layout()
    387         fig.savefig("profile.png", dpi=100)
    388         
    389 
    390 
    391 if __name__ == "__main__":
    392     ssvcObj = SSVC(X, c=100, mu=7, beta=300)
    393     alpha, R, tab = ssvcObj.optimize()
    394     spObj = SSVCPlot()
    395     spObj.profile_plot(X, ssvcObj, alpha, R)
    396     
    View Code
  • 结果展示:
    左侧为聚类轮廓, 右侧为簇划分. 很显然, 此处在合适的超参数选取下, 将原始数据集划分为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.
原文地址:https://www.cnblogs.com/xxhbdk/p/12586604.html