In mmfd99/StatComp21092: Final

knitr::opts_chunk$set(echo = TRUE)

1 K-SMOTE

1.1 SMOTE

SMOTE(Synthetic Minority Oversampling Technique)是一种基于线性插补的过采样方法。

其想法是对某一少数类样本点与少数类k近邻之一进行线性插补，得到新的少数类样本，从而达到数据增强的目的。

下面给出了数学表达式：

$$\boldsymbol{x}{syn}=\boldsymbol{x}_p+\delta(\boldsymbol{x}_q-\boldsymbol{x}_p),\quad \boldsymbol{x}_p,\boldsymbol{x}_q\in S{min},\quad \delta\sim U(0,1) $$

这里$\boldsymbol{x}_q$为$\boldsymbol{x}_p$的少数类k近邻之一。

1.2 核方法

核方法是被广泛应用的方法，其典型算法是支持向量机SVM与RBF-NN。

其想法是利用内积代替映射的显式表达，通过改变空间的度量结构达到简化学习任务的作用。

下面给出了核函数的定义：

$$K(\boldsymbol{x}_i,\boldsymbol{x}_j)=\langle \phi(\boldsymbol{x}_i),\phi(\boldsymbol{x}_j) \rangle$$ 这里$\boldsymbol{x}_i,\boldsymbol{x}_j$为特征空间中的点，$\boldsymbol{K}$为核函数，为再生核希尔伯特空间上的内积。

再生核希尔伯特空间上由核函数诱导的距离如下：

$$ d^{\phi}(\boldsymbol{x}_i,\boldsymbol{x}_j)=\Vert\phi(\boldsymbol{x}_i)-\phi(\boldsymbol{x}_j)\Vert^2=K(\boldsymbol{x}_i,\boldsymbol{x}_i)+K(\boldsymbol{x}_j,\boldsymbol{x}_j)-2K(\boldsymbol{x}_i,\boldsymbol{x}_j) $$

1.3 基于核函数的SMOTE

由于样本点间往往不具备线性结构，所以做简单的线性插补不能达到优化学习模型的效果。

基于这样的想法，我们可以构建基于核函数的SMOTE。

$$ \phi(\boldsymbol{x}_{syn})=\phi(\boldsymbol{x}_p)+\delta(\phi(\boldsymbol{x}_p)-\phi(\boldsymbol{x}_q)),\quad \delta\sim U(0,1) $$

这里为了方便表达，我们利用了$\phi$这一显式映射，事实上，在后续的算法实现过程中，$\phi$是不会被表示出来的，我们会利用其内积，即核函数$\boldsymbol{K}$来达到构建再生核希尔伯特空间的作用。

2 SVM

SVM是核方法的一个典型应用。

在经典的SVM算法中，我们需要解决如下优化问题：

$$ \begin{align} &\min_{\alpha}\quad \frac{1}{2}\sum_{i,j=1}^N\alpha_i\alpha_jy_iy_j K(x_i,x_j)-\sum_{i=1}^N\alpha_i\ &s.t \quad \sum_{i=1}^N\alpha_iy_i=0\ &\quad \quad 0\le\alpha_i\le C,\quad i=1,2,\cdots,N \end{align} $$

利用KKT条件，问题的解为：

$$ \begin{align} &w^{\ast}=\sum_{i=1}^{N}\alpha_i^{\ast}y_ix_i\ &b^{\ast}=y_j-\sum_{i=1}^N\alpha_i^{\ast}y_i K(x_i, x_j) \end{align} $$

分类决策函数为：

$$ \begin{align} f(x)&=sign(\sum_{i=1}^{N_s}\alpha_i^{\ast}y_iK(x_i,x)+b^{\ast}) \end{align} $$

无论是求解过程还是分类决策函数都与映射$\phi$无关，这意味着我们关心核矩阵即可。

3 计算核矩阵

计算Gram matrix K。这一步是核心，连接SMOTE和SVM。

$$ K=(K(\boldsymbol{x}i,\boldsymbol{x}_j)){(N+P)\times(N+P)}= \begin{bmatrix} K_1&K_2\ K_2^T&K_3 \end{bmatrix} $$

这里，$K_1=(K(\boldsymbol{x}i,\boldsymbol{x}_j)){N\times N},\quad \boldsymbol{x}_i,\boldsymbol{x}_j\in S$

$K_2=(K(\boldsymbol{x}i,\boldsymbol{x}_j)){N\times P},\quad \boldsymbol{x}_i\in S,\boldsymbol{x}_j\in S^{syn}$

$K_3=(K(\boldsymbol{x}i,\boldsymbol{x}_j)){P\times P},\quad \boldsymbol{x}_i,\boldsymbol{x}_j\in S^{syn}$

此外，为了方便后面的表示，我们用$x_i^{pq}$来表示$\boldsymbol{x}_p,\boldsymbol{x}_q$生成的样本，$\delta^{pq}$表示对应的随机数。

下面，我们可以给出$K_2,K_3$中元素的具体表达。(我们需要将其表示成原样本集与核函数的形式，这样我们就避免了显式生成新样本，说白了，我们生成的不是样本，而是$K(\boldsymbol{x}_i,\boldsymbol{x}_j)$，这就是核方法的精彩了)

$K_2$中元素：

$$ K(\boldsymbol{x}_i,\boldsymbol{x}_j^{pq})=(1-\delta^{pq})K(\boldsymbol{x}_i,\boldsymbol{x}_p)+\delta^{pq}K(\boldsymbol{x}_i,\boldsymbol{x}_q),\quad \boldsymbol{x}_i\in S,\boldsymbol{x}_j^{pq}\in S^{syn} $$

$K_3$中元素：

$$ K(\boldsymbol{x}_i^{lm},\boldsymbol{x}_j^{pq})=(1-\delta^{pq})(1-\delta^{lm})K(\boldsymbol{x}_p,\boldsymbol{x}_l)+(1-\delta^{pq})(\delta^{lm})K(\boldsymbol{x}_p,\boldsymbol{x}_m)+(\delta^{pq})(1-\delta^{lm})K(\boldsymbol{x}_q,\boldsymbol{x}_l)+(\delta^{pd})(\delta^{lm})K(\boldsymbol{x}_q,\boldsymbol{x}_m) $$

这里，值得一提的是，我们的样本生成与核函数的具体表示无关，所以这一SMOTE方法适用于所有的核函数。

下面是实现，分两步，一个是K-SMOTE，还有一个是基于SMO算法的SVM分类器。

4 K-SMOTE

首先加载训练集和测试集。

library(StatComp21092)
data("yeast_tra")
data("yeast_test")
train_set=yeast_tra[1:500,]
test_set=yeast_test

接下来我们利用前面的公式计算核矩阵。

#定义三阶多项式核函数
kernel_poly<-function(x,y)
{
  #x，y为特征空间上的向量
  return((x%*%y+1)^3)
}


#定义再生核希尔伯特空间上的距离平方
distance<-function(x,y)
{
  #x，y为特征空间上的向量
  d_square=kernel_poly(x,x)+kernel_poly(y,y)-2*kernel_poly(x,y)
  return(d_square)
}

#K_SMOTE，同时计算新的核矩阵与分类决策函数的核因子
K_SMOTE<-function(train_set,test_set,P,k)
{
  #train_set,test_set分别为训练集和测试集，类型为dataframe
  #P为需要生成的少数类样本数
  #k为kNN中近邻的个数


  #记录训练集中少数类和多数类的index
  min_index=(which(train_set['Class']==1))
  max_index=(which(train_set['Class']==0))


  #初始化核矩阵
  N = nrow(train_set)
  K = matrix(0,N+P,N+P)
  D = matrix(0,N+P,N+P)

  #初始化生成对与随机数
  p_index=numeric(length = P)
  q_index=numeric(length = P)
  delta=runif(P,0,1)

  #计算原核矩阵
  for (i in 1:N)
     for (j in 1:N)
       K[i,j]=kernel_poly(data.matrix(train_set)[i,3:ncol(train_set)-1],data.matrix(train_set)[j,3:ncol(train_set)-1])

  #计算距离矩阵(我们这里只关心少数类)
  for (i in min_index)
    for (j in min_index)
       D[i,j]=distance(data.matrix(train_set)[i,3:ncol(train_set)-1],data.matrix(train_set)[j,3:ncol(train_set)-1])



  #用少数类的距离矩阵结合KNN得到生成对
  for (i in 1:P)
  {
    p=sample(min_index,1,replace = TRUE)
    r=min_index[rank(D[p,min_index])]
    q=sample(r[1:k],1,replace = TRUE)
    p_index[i]=p
    q_index[i]=q
  }



  #生成，即填补核矩阵

  #填K2和K2转置：一个为原样本，一个为生成样本
  for(i in 1:P)
    for(j in 1:N)
    {
      p=p_index[i]
      q=q_index[i]
      K[i+N,j]=(1-delta[i])*K[p,j]+delta[i]*K[q,j]
      K[j,i+N]=K[i+N,j]
    }

  #填K3：均为生成样本
  for(i in 1:P)
    for(j in 1:P)
    {
      l=p_index[i]
      m=q_index[i]
      p=p_index[j]
      q=q_index[j]
      K[i+N,j+N]=(1-delta[j])*(1-delta[i])*K[p,l]+(1-delta[j])*delta[i]*K[p,m]+delta[j]*(1-delta[i])*K[q,l]+delta[j]*delta[i]*K[q,m]
    }


  # 下面填测试样本与训练集的核矩阵
  n=nrow(test_set)
  gram_test=matrix(0,n,N+P)

  for(i in 1:n)
    for(j in 1:N)
      gram_test[i,j]=kernel_poly(data.matrix(test_set)[i,3:ncol(train_set)-1],data.matrix(train_set)[j,3:ncol(train_set)-1])

  for(i in 1:n)
    for(j in 1:P)
    {
      p=p_index[j]
      q=q_index[j]
      gram_test[i,j+N]=(1-delta[j])*gram_test[i,p]+delta[j]*gram_test[i,q]
    }

  mat=list("train_km"=K,"test_km"=gram_test)

  return(mat)
}

由于R包e1071的svm函数的核函数参数中，只有rbf，poly，linear和nn这四种核函数，没有提供我们需要的核矩阵选项，所以我们还需要写一个基于SMO算法的SVM分类器。

5 SMO算法

我们之前提到了SVM算法是解决如下凸优化问题：

$$ \begin{align} &\min_{\alpha}\quad \frac{1}{2}\sum_{i,j=1}^N\alpha_i\alpha_jy_iy_j K(x_i,x_j)-\sum_{i=1}^N\alpha_i\ &s.t \quad \sum_{i=1}^N\alpha_iy_i=0\ &\quad \quad 0\le\alpha_i\le C,\quad i=1,2,\cdots,N \end{align} $$

但随着样本量的增加，参数量会变得异常庞大。因此，我们不能直接对这个问题用梯度下降进行求解。

SMO的想法是我们每次只优化两个变量$\alpha_i,\alpha_j$，将其余变量视为常数，同时由于有等式约束，我们可以得到$\alpha_i$与$\alpha_j$之间的关系，进而可以求解。基于这样的想法，SMO算法将一个多变量的优化问题转化为多个简单的二次规划问题。

#预处理
#得到两个核矩阵
P=400
mat=K_SMOTE(train_set,test_set,P,3)
train_matrix=mat$train_km
test_matrix=mat$test_km

print(train_matrix[1:5,1:5])
print(test_matrix[1:5,1:5])

#提取label并修正为1和-1
modify_label=function(train_set, test_set, P)
{
  train_target=array(1, dim = nrow(train_set)+P)
  test_target=array(1, dim = nrow(test_set))

  for(i in 1:nrow(train_set))
    if(train_set[i,'Class']==0)
    {
      train_target[i]=-1
    }
    if(train_set[i,'Class']==1)
    {
      train_target[i]=1
    }

  for(i in 1:nrow(test_set))
    if(test_set[i,'Class']==0)
    {
      test_target[i]=-1
    }
    if(test_set[i,'Class']==1)
    {
      test_target[i]=1
    }
  target=list("train"=train_target,"test"=test_target)

  return(target)
}

all_target=modify_label(train_set,test_set,P)

update_alpha_2=function(alpha_2,L,H)
{
  if(alpha_2>H)
    return(H)
  else
    if(alpha_2<L)
      return(L)
    else
      return(alpha_2)
}

SMO_SVC=function(train_mat,test_mat,train_target,test_target,iter_max=10,epsilon=0.01,C=1)
{
  #初始化
  N=nrow(train_mat)
  alpha=numeric(length = nrow(train_mat))
  b=1
  k=0#迭代次数
  #定义函数g
  g=function(i)
  {
    return(sum(alpha*train_target*train_mat[i,])+b)
  }

  sample_without_i=function(i)
  {
    j=sample(1:N,1)
    while(j==i)
      j=sample(1:N,1)
    return(j)
  }

  while (k<iter_max) 
  {
    change=0
    #取出i，j
    for(i in 1:N)
    {
      a_i=alpha[i]
      y_i=train_target[i]
      x_i=i
      E_i=g(x_i)-y_i

      j=sample_without_i(i)
      a_j=alpha[j]
      y_j=train_target[j]
      x_j=j
      E_j=g(x_j)-y_j

      eta=train_mat[i,i]+train_mat[j,j]-2*train_mat[i,j]
      if(eta<=0)
        next

      #未经剪辑的new_a_j
      new_a_j=a_j+y_j*(E_i-E_j)/eta

      #剪辑a_j_new
      if(y_i!=y_j)
      {
        L=max(0,a_j-a_i)
        H=min(C,C+a_j-a_i)
      }
      else
      {
        L=max(0,a_j+a_i-C)
        H=min(C,a_i+a_j)
      }

      a_j_new=update_alpha_2(new_a_j,L,H)
      a_i_new=a_i+y_i*y_j*(a_j-a_j_new)

      if(abs(a_j_new-a_j)<epsilon)
        next

      #更新alpha
      alpha[i]=a_i_new
      alpha[j]=a_j_new

      #更新b
      b_i=-E_i-y_i*train_mat[i,i]*(a_i_new - a_i) - y_j*train_mat[i,j]*(a_j_new - a_j) + b
      b_j=-E_j-y_i*train_mat[i,j]*(a_i_new - a_i) - y_j*train_mat[j,j]*(a_j_new - a_j) + b

      if(abs(a_i_new-C/2)<C/2)
        b=b_i
      else
        if(abs(a_j_new-C/2)<C/2)
          b=b_j
        else
          b=(b_i+b_j)/2

      change=change+1
      }
      if(change==0)
        k=k+1
      else
        k=0
  }

  #计算测试集结果

  result=numeric(length = nrow(test_mat))
  for (i in 1:nrow(test_mat)) 
  {
    f=sum(alpha*train_target*test_mat[i,])+b
    if(f>0)
      result[i]=1
    else
      result[i]=-1
  }
  return(result)
}
pred=SMO_SVC(train_matrix,test_matrix,all_target$train,all_target$test,1e1,1e-2,1)
target=all_target$test

ConfusionMatrix=function(pred,target)
{
  N=length(pred)
  TP=FP=FN=TN=0
  for(i in 1:N)
  {
    if(pred[i]==1 &&target[i]==1)
      TP=TP+1
    if(pred[i]==-1&&target[i]==-1)
      TN=TN+1
    if(pred[i]==1&&target[i]==-1)
      FP=FP+1
    if(pred[i]==-1&&target[i]==1)
      FN=FN+1
  }

  mat=matrix(c(TP,FP,FN,TN),ncol=2)
  return(mat)
}
ConfusionMatrix(pred,target)

mmfd99/StatComp21092 documentation built on Dec. 23, 2021, 10:14 p.m.