保持结构的非监督特征选择（Structure Preserving Unsupervised Feature Selection）--- 机器学习大作业

最新推荐文章于 2023-01-21 15:23:29 发布

McQueen_LT

最新推荐文章于 2023-01-21 15:23:29 发布

阅读量825

点赞数 1

分类专栏：机器学习文章标签：算法聚类 python 机器学习

本文链接：https://blog.csdn.net/McQueen_LT/article/details/117114442

版权

机器学习专栏收录该内容

5 篇文章 2 订阅

订阅专栏

一、概述

机器学习大作业。机器学习课程最后要求实现一篇论文，这里做一下记录。

二、实验内容

这里只放上最重要的部分。

本次课程论文需要实现的算法如算法1所示：
在这里插入图片描述

2.1 权重矩阵预计算

根据算法1可知，输入数据为原始数据 $X\in\R^{n\times d}$ ，权重矩阵 $S$ ，参数 $\alpha$ 和 $\beta$ ，以及需要从原始数据中筛选出的特征数 $h$ ，而输出为 $h$ 维被筛选出来的特征，即 $X\in \R^{n\times h}$ 。此外，权重矩阵 $S$ 需要根据原始数据 $X$ 预先计算得到。

根据公式 $8$ 可知，权重矩阵 $S$ 为所有 $n$ 个样本数据循环遍历，依次经径向基核 $k(x_i,x_j)=e^{-\frac{\|x_i-x_j\|^2}{\sigma}}$ 计算得到的核矩阵，即：
$S=\left[\begin{matrix} k(x_1,x_1)&k(x_1,x_2)&k(x_1,x_3)&\cdots&k(x_1,x_n)\\ k(x_2,x_1)&k(x_2,x_2)&k(x_2,x_3)&\cdots&k(x_2,x_n)\\ k(x_3,x_1)&k(x_3,x_2)&k(x_3,x_3)&\cdots&k(x_3,x_n)\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ k(x_n,x_1)&k(x_n,x_2)&k(x_n,x_3)&\cdots&k(x_n,x_n) \end{matrix}\right]$

在实现过程中 $\sigma=1$ 。计算权重矩阵 $S$ 的代码如下：

def matrixGaussianKernel_S(X):
    num_samples = len(X)
    S = np.zeros((num_samples,num_samples))
    for i in range(num_samples):
        for j in range(num_samples):
            S[i][j] = math.exp(-math.pow(np.linalg.norm(X[i]-X[j]),2)/1.0)
    return S

其中参数 $X$ 为输入的原始数据 $X\in\R^{n\times d}$ ，返回值为核矩阵 $S$ 。

2.2 特征提取算法实现

算法1具体实现步骤如下：

步骤1：初始化矩阵 $Q = I$ ，其中 $I$ 为单位矩阵， $Q\in\R^{d\times d}$ ， $Q$ 初始化代码如下：

def initializeIdeMatrix_Q(num_feature):
    I = np.zeros((num_feature,num_feature))
    for i in range(num_feature):
        I[i][i] = 1
    return I

参数 $num\_feature$ 为原始数据 $X$ 所有特征的维度 $d$ 。

步骤2：计算 $L a p l a c e$ 矩阵 $L_S=D-S$ ，其中 $D$ 为对角矩阵，并且每个对角元素为权重矩阵 $S$ 的各行元素之和， $D$ 可如下表示：
$D=\left[\begin{matrix} sum(S_1)&0&0&\cdots&0\\ 0&sum(S_2)&0&\cdots&0\\ 0&0&sum(S_3)&\cdots&0\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 0&0&0&\cdots&sum(S_n) \end{matrix}\right]$
函数 $sum(S_j)$ 表示对矩阵 $S$ 第 $j$ 行所有元素求和，即
$D_{ii}=sum(S_i)=\sum_{j=1}^{n}S_{ij}$
由上文已知矩阵 $S$ 为径向基核矩阵，所以即可求得 $L a p l a c e$ 矩阵，实现代码如下：

def calculateMatrix_Ls(S):
    D = np.zeros((S.shape[0],S.shape[0]))
    for i in range(S.shape[0]):
        D[i][i] = np.sum(S[i])
    return D-S

步骤3：迭代计算表达矩阵 $W$ ，不断训练更新矩阵 $Q$ ，直到矩阵 $W$ 收敛。表达矩阵 $W$ 的计算公式为：
$W=(\beta X^TL_SX+XX^T+\alpha Q)^{-1}X^TX$
代码实现为：

def iterationUntilConvergence(alpha,beta,epsilon,Ls,Q,X):
    #temp_W = np.zeros((X.shape[1],X.shape[1]))
    '''W = np.matmul(np.matmul(
        np.linalg.inv(np.matmul(
            np.matmul(beta * X.T, Ls), X
        ) + np.matmul(X.T, X) + alpha * Q), X.T), X)'''
    #Q = updateQ(epsilon, W)
    
    #itetimes = 0
    #while(np.sum(np.abs(temp_W-W)) > 0.01):
    for i in range(50):
        #print(np.sum(np.abs(temp_W-W)))	#输出每迭代一次W与上一次W之间的差
        #temp_W = W
        #itetimes += 1  #统计W达到收敛精度0.01时的迭代次数
        W = np.matmul(np.matmul(
            np.linalg.inv(np.matmul(
                np.matmul(beta*X.T,Ls),X
            )+np.matmul(X.T,X)+alpha*Q),X.T),X)
        Q = updateQ(epsilon, W)

    #print(itetimes)
    return W

根据论文实验结果，迭代训练50次即可，其中，利用 $u p d a t e Q ()$ 方法来训练更新矩阵 $Q$ ，计算 $Q$ 的公式为：
$Q_{ii}=\frac{1}{2\sqrt{W^T_iW_i+\varepsilon}}$
其中 $\varepsilon\rightarrow0$ ，实现过程中取 $\varepsilon=0.0001$ 。 $u p d a t e Q ()$ 定义如下：

def updateQ(epsilon,W):
    Q = np.zeros((W.shape[0],W.shape[0]))
    for i in range(W.shape[0]):
        Q[i][i] = 1/(2*math.sqrt(np.sum(np.square(W[i]))+epsilon))
    return Q

步骤4：令 $r^{d\times 1}=[r_1,r_2,\cdots,r_d]^T$ ，分别计算表达矩阵 $W$ 各行元素的2范数 $r_i=\|W_i\|_2$ ，并将 $r^{d\times 1}$ 降序排列，取最大的前 $h$ 个 $r_i$ 的索引 $index=\{i_1,i_2,\cdots,i_h\}$ ，则原始数据 $X$ 被提取后的特征数据为 $feature\_X=\{X^T_{i_1},X^T_{i_2},\cdots,X^T_{i_h}\}^T$ 。实现过程如下：

def rankBasedW(h,W,X):
    norm_W = np.linalg.norm(W,axis=1)
    index = heapq.nlargest(h,range(len(norm_W)),norm_W.take) #取前h个最大值的索引
    print(index)
    feature_X = np.zeros((X.shape[0],h))
    for i in range(X.shape[0]):
        feature_X[i] = X[i][index] #取每个数据中被提取的特征数据
    return feature_X

以上为算法1的实现过程，还需要对 $X$ 经过提取特征之后的数据 $feature\_X$ 进行实验结果验证。

2.3 标签重排与精度计算

采用 $K M e a n s$ 聚类算法，对 $feature\_X$ 进行聚类，聚类方法可以直接利用 $s k l e a r n$ 机器学习库方法：

def k_means(data, clusters):
    return KMeans(n_clusters=clusters,random_state=0).fit(data).predict(data)

lable_pred = k_means(feature_X,len(np.unique(Y)))

其中 $label\_pred$ 为经过聚类后得到的数据标签。由于数据标签的排列标准不一致，需要利用 $M u n k r e s$ 算法对 $label\_pred$ 进行重排，重排过程实现如下：

from munkres import Munkres
import numpy as np

def maplabels(L1, L2):
    L2 = L2+1
    Label1 = np.unique(L1)
    Label2 = np.unique(L2)
    nClass1 = len(Label1)
    nClass2 = len(Label2)
    nClass = np.maximum(nClass1, nClass2)
    G = np.zeros((nClass, nClass))
    for i in range(nClass1):
        ind_cla1 = L1 == Label1[i]
        ind_cla1 = ind_cla1.astype(float)
        for j in range(nClass2):
            ind_cla2 = L2 == Label2[j]
            ind_cla2 = ind_cla2.astype(float)
            G[i, j] = np.sum(ind_cla2*ind_cla1)

    m = Munkres()
    index = m.compute(-G.T)
    index = np.array(index)
    index = index+1
    print(-G.T)
    print(index)
    newL2 = np.zeros(L2.shape, dtype=int)
    for i in range(nClass2):
        for j in range(len(L2)):
            if L2[j] == index[i, 0]:
                newL2[j] = index[i, 1]

    return newL2

传入参数 $L_1,L_2$ 分别为原始数据标签和经过聚类得到的预测标签，返回值为经过重排后的标签。分别使用 $A C$ 和 $N M I$ 指标评价聚类精度， $A C$ 实现如下：

import numpy as np

def acc(L1, L2):
    sum = np.sum(L1[:]==L2[:])
    return sum/len(L2)

传入参数 $L_1,L_2$ 分别为原始数据标签和经过 $M u n k r e s$ 算法重排后标签，返回值为 $A C$ 精度。调用 $s k l e a r n$ 机器学习库方法计算 $N M I$ 评价指标：

from sklearn import metrics

def nmi(L1, L2):
    return metrics.normalized_mutual_info_score(L1, L2)

传入参数同上，返回值为 $N M I$ 精度。

三、实验代码

README：

运行环境
Python 3.7 64位
pycharm
package：sklearn, numpy, math，heapq等
将code文件夹中所有.py文件复制到pycharm中，运行spufs.py文件即可。
注：运行之前，需要安装相关包，并将代码中的相关路径名进行修改。

Acc.py

import numpy as np

def acc(L1, L2):
    sum = np.sum(L1[:]==L2[:])
    return sum/len(L2)

datadvi.py

from scipy.io import loadmat
import numpy as np

def divdata():
    filename = '.../作业/机器学习/datasets/' + input("input name of data file: ")
    data = loadmat(filename)
#    print(data['X'])


    if filename == '.../作业/机器学习/datasets/COIL20.mat':
        dataX = data['fea']
        dataY = data['gnd'][0]

    else:
        dataX = data['X']
        dataY = data['Y'].T[0]

    print(len(dataX[0]))

    divideornot = input("divide data or not?(Yes/No): ")
    if divideornot == 'Yes':
        dataX_train = []
        dataX_predict = []
        dataY_train = []
        dataY_predict = []
        num_Y = np.unique(dataY).astype(int)
        for i in range(len(num_Y)):
            temp = dataY == num_Y[i]
            temp.astype(float)
            num_Y[i] = np.sum(temp)
            flag = 0
            for j in range(len(dataY)):
                if temp[j] == 1:
                    if flag < int(round(0.9 * num_Y[i])):
                        dataX_train.append(dataX[j])
                        dataY_train.append(dataY[j])
                        flag += 1
                    else:
                        dataX_predict.append(dataX[j])
                        dataY_predict.append(dataY[j])

        dataX_train = np.array(dataX_train)
        dataX_predict = np.array(dataX_predict)
        dataY_train = np.array(dataY_train)
        dataY_predict = np.array(dataY_predict)
        return dataX_train,dataX_predict,dataY_train,dataY_predict
    else:
        return dataX,dataX,dataY,dataY

def decreaseData(dataX,dataY):
    dataX_train = []
    dataY_train = []
    num_Y = np.unique(dataY).astype(int)
    print("this data has {} samples".format(len(dataX)))
    ratio = float(input("input the ratio you want to decrease: "))
    for i in range(len(num_Y)):
        temp = dataY == num_Y[i]
        temp.astype(float)
        num_Y[i] = np.sum(temp)
        flag = 0
        for j in range(len(dataY)):
            if temp[j] == 1:
                if flag < round(ratio * num_Y[i]):
                    dataX_train.append(dataX[j])
                    dataY_train.append(dataY[j])
                    flag += 1

    dataX_train = np.array(dataX_train)
    dataY_train = np.array(dataY_train)
    print(dataX_train)

    return dataX_train,dataY_train

kmeans.py

from sklearn.cluster import KMeans

def k_means(data, clusters):
    return KMeans(n_clusters=clusters,random_state=0).fit(data).predict(data)

maplabels.py

from munkres import Munkres, print_matrix
import numpy as np

def maplabels(L1, L2):
    L2 = L2+1
    Label1 = np.unique(L1)
    Label2 = np.unique(L2)
    nClass1 = len(Label1)
    nClass2 = len(Label2)
    nClass = np.maximum(nClass1, nClass2)
    G = np.zeros((nClass, nClass))
    for i in range(nClass1):
        ind_cla1 = L1 == Label1[i]
        ind_cla1 = ind_cla1.astype(float)
        for j in range(nClass2):
            ind_cla2 = L2 == Label2[j]
            ind_cla2 = ind_cla2.astype(float)
            G[i, j] = np.sum(ind_cla2*ind_cla1)

    m = Munkres()
    index = m.compute(-G.T)
    index = np.array(index)
    index = index+1
    print(-G.T)
    print(index)
    newL2 = np.zeros(L2.shape, dtype=int)
    for i in range(nClass2):
        for j in range(len(L2)):
            if L2[j] == index[i, 0]:
                newL2[j] = index[i, 1]

    return newL2

NMI.py

from sklearn import metrics

def nmi(L1, L2):
    return metrics.normalized_mutual_info_score(L1, L2)

spufs.py

import numpy as np
import math
import heapq
import datadvi
import kmeans
import maplabels
import Acc
import NMI

def matrixGaussianKernel_S(X):
    num_samples = len(X)
    S = np.zeros((num_samples,num_samples))
    for i in range(num_samples):
        for j in range(num_samples):
            S[i][j] = math.exp(-math.pow(np.linalg.norm(X[i]-X[j]),2)/1.0)
    return S

def initializeIdeMatrix_Q(num_feature):
    I = np.zeros((num_feature,num_feature))
    for i in range(num_feature):
        I[i][i] = 1
    return I

def calculateMatrix_Ls(S):
    D = np.zeros((S.shape[0],S.shape[0]))
    for i in range(S.shape[0]):
        D[i][i] = np.sum(S[i])
    return D-S

def updateQ(epsilon,W):
    Q = np.zeros((W.shape[0],W.shape[0]))
    for i in range(W.shape[0]):
        Q[i][i] = 1/(2*math.sqrt(np.sum(np.square(W[i]))+epsilon))
    return Q

def objectiveFunc(alpha, beta, epsilon, W, Ls, X):
    medTermValue = 0
    for i in range(W.shape[0]):
        medTermValue += math.sqrt(np.sum(np.square(W[i]))+epsilon)
    medTermValue *= alpha
    lastTerm = np.matmul(np.matmul(np.matmul(np.matmul(W.T,X.T),Ls),X),W)
    valueOfObjFun = math.pow(np.linalg.norm(X-np.matmul(X,W)),2) + medTermValue + beta*np.trace(lastTerm)
    return valueOfObjFun

def iterationUntilConvergence(alpha,beta,epsilon,Ls,Q,X):
    temp_W = np.zeros((X.shape[1],X.shape[1]))
    W = np.matmul(np.matmul(
        np.linalg.inv(np.matmul(
            np.matmul(beta * X.T, Ls), X
        ) + np.matmul(X.T, X) + alpha * Q), X.T), X)
    Q = updateQ(epsilon, W)
    #for i in range(50):
    itetimes = 0
    #while(np.sum(np.abs(temp_W-W)) > 0.01):
    for i in range(50):
        print(objectiveFunc(alpha,beta,epsilon,W,Ls,X))
        #print(np.sum(np.abs(temp_W-W)))
        temp_W = W

        itetimes += 1

        W = np.matmul(np.matmul(
            np.linalg.inv(np.matmul(
                np.matmul(beta*X.T,Ls),X
            )+np.matmul(X.T,X)+alpha*Q),X.T),X)
        Q = updateQ(epsilon, W)
        #print(np.sum(np.abs(temp_W - W)))
        #print(Q)
    print(itetimes)
    return W

def rankBasedW(h,W,X):
    norm_W = np.linalg.norm(W,axis=1)
    index = heapq.nlargest(h,range(len(norm_W)),norm_W.take)
    print(index)
    feature_X = np.zeros((X.shape[0],h))
    for i in range(X.shape[0]):
        feature_X[i] = X[i][index]
    return feature_X


if __name__ == '__main__':
    alpha = float(input("input parameter alpha: "))
    beta = float(input("input parameter beta: "))
    epsilon = float(input("input parameter epsilon: "))
    h = int(input("input number of features h: "))
    X,X_pred,Y,Y_pred = datadvi.divdata()
    Q = initializeIdeMatrix_Q(X.shape[1])
    S = matrixGaussianKernel_S(X)
    Ls = calculateMatrix_Ls(S)
    W = iterationUntilConvergence(alpha,beta,epsilon,Ls,Q,X)
    feature_X = rankBasedW(h,W,X)

    print(feature_X.shape)
    lable_pred = kmeans.k_means(feature_X,len(np.unique(Y)))
    lable_pred = maplabels.maplabels(Y,lable_pred)
    print(Acc.acc(Y,lable_pred))
    print(NMI.nmi(Y,lable_pred))

'''
a = np.array([1,2,3])
b = np.array([4,5,6])
print(math.pow(np.linalg.norm(a),2))
print(np.sum(np.square(a)))
print(heapq.nlargest(2,range(len(a)),a.take))
print(a[[1,2]])
'''

McQueen_LT

关注

1
点赞
踩
7

收藏

觉得还不错? 一键收藏
7
评论
保持结构的非监督特征选择（Structure Preserving Unsupervised Feature Selection）--- 机器学习大作业

一、概述机器学习大作业。机器学习课程最后要求实现一篇论文，这里做一下记录。二、实验内容这里只放上最重要的部分。本次课程论文需要实现的算法如算法1所示：2.1 权重矩阵预计算根据算法1可知，输入数据为原始数据X∈Rn×dX\in\R^{n\times d}X∈Rn×d，权重矩阵SSS，参数α\alphaα和β\betaβ，以及需要从原始数据中筛选出的特征数hhh，而输出为hhh维被筛选出来的特征，即X∈Rn×hX\in \R^{n\times h}X∈Rn×h。此外，权重矩阵SSS需要根据原始数
复制链接

扫一扫