SOM网络2： 代码的实现

SOM自组织映射神经网络的原理，详见博客：SOM网络1:原理讲解

训练的主函数

train_SO代码如下：

def train_SOM(X,                                                  # 输出节点行数Y,                                                  # 输出节点列数N_epoch,                                            # epochdatas,                                              # 训练数据(N x D)  N个D维样本init_lr=0.5,                                        # 初始化学习率  lrsigma = 0.5,                                        # 初始化 sigma 用来更新领域节点权重dis_func = euclidean_distance,                      # 距离公式  默认欧拉距离neighborhood_func = gaussion_neighborhood,          # 邻域节点权重公式g  默认高斯函数init_weight_fun=None,                               #初始化权重函数seed=10):  			# 获取输入的特征维度N,D =np.shape(datas)# 训练的步数N_steps =N_epoch*N#对权重进行初始化rng = np.random.RandomState(seed)if init_weight_fun is None:weights =rng.rand(X,Y,D)*2-1                            #随机初始化weights /=np.linalg.norm(weights,axis=-1,keepdims=True) #标准化else:weights = init_weight_fun(X,Y,datas)                       # 一般使用PCA初始化

PCA 初始化权重

def weights_PCA(X,Y,data):N,D=np.shape(data)weights=np.zeros([X,Y,D])pc_value,pc=np.linalg.eig(np.cov(np.transpose(data)))   # pc_vale为特征值,pc 为特征向量 DXD维pc_order=np.argsort(-pc_value)  # 特征值从大到小排序，并返回Index# 对W:[X,Y,D]进行初始化for i,c1 in enumerate(np.linspace(-1,1,X)):for j,c2 in enumerate(np.linsapce(-1,1,Y)):weights[i,j]=c1*pc[pc_order[0]]+c2*pc[pc_order[1]]   #利用最大的2个特征值对应的特征向量加权组合成i,j位置的D维表征向量	

完整的训练代码

def train_SOM(X,                                                  # 输出节点行数Y,                                                  # 输出节点列数N_epoch,                                            # epochdatas,                                              # 训练数据(N x D)  N个D维样本init_lr=0.5,                                        # 初始化学习率  lrsigma = 0.5,                                        # 初始化 sigma 用来更新领域节点权重dis_func = euclidean_distance,                      # 距离公式  默认欧拉距离neighborhood_func = gaussion_neighborhood,          # 邻域节点权重公式g  默认高斯函数init_weight_func=weights_PCA,                        #初始化权重函数seed=10):  			# 获取输入的特征维度N,D =np.shape(datas)# 训练的步数N_steps =N_epoch*N#对权重进行初始化rng = np.random.RandomState(seed)if init_weight_func is None:weights =rng.rand(X,Y,D)*2-1                            #随机初始化weights /=np.linalg.norm(weights,axis=-1,keepdims=True) #标准化else:weights = init_weight_fun(X,Y,datas)                       # 一般使用PCA初始化for n_epoch in range(N_epoch):print("Epoch %d" %(n_epoch+1))#打乱样本次序index=rng.permulation(np.arange(N))for n_step,_id in enumerate(index):# 取一个样本x=datas[_id]#计算learning rate (eta)t=N*n_epoch + n_stepeta=get_learning_rate(init_lr,t,N_steps)#计算样本距离输出的每个节点的距离，并获取激活点的位置winner=get_winner_index(x,weights,dis_func)#根据激活点的位置计算临近点的权重   随着迭代的进行sigma也需要不断减少new_sigma=get_learning_rate(sigma,t,N_steps)  # sigma 更新的方式和学习率一样g=neighborhood_fun(X,Y,winner,new_sigma) g=g*eta#进行权重的更新weights =  weights + np.expand_dims(g,-1)*(x-weights)    # 打印量化误差print("quantization_error=%.4f" %(get_quantization_error(data,weights))) return weights#计算学习率
def get_learning_rate(lr,t,max_steps):  # t当前的steps   max_steps=N x epoch  (N样本数)		return lr/(1+t/(max_steps/2))	# 获取激活(获胜点)节点的位置，与x距离最小的输出节点位置
def get_winner_index(x,w,dis_func=euclidean_distance):# 计算输入样本和各个节点的距离dis = dis_func(x,w)#找到距离最小的位置index=np.where(dis ==np.min(dis))return (index[0][0],index[1][0])#利用高斯距离法计算临近点的权重
# X,Y模板大小,c中心点的位置  def gaussion_neighborhood(X,Y,c,sigma)xx,yy=np.meshgrid(np.arange(X),np.arange(Y))d=2*sigma*sigmaax=np.exp(-np.power(xx-xx.T[c],2)/d)ay=np.exp(-np.power(yy-yy.T[c],2)/d)return (ax*ay).T# 计算欧式距离
def euclidean_distance(x,w):dis=np.expand_dims(x,axis=(0,1))-w   # x:D w:[X,Y,D]  因此需要增加两维 x:D->x:[1,1,D]return np.linalg.norm(dis,axis=-1)                   # 输出[X,Y] 二范数 即为欧拉距离# 特征标准化 (x-mu)/std
def feature_normalization(data):mu=np.mean(data,axis=0,keepdims=True)sigma=np.std(data,axis=0,keepdims=True)return (data-mu)/sigmadef get_U_Matrix(weights):X,Y,D=np.shape(weights)um=na.nan * np.zeros((X,Y,8))  #8 领域ii=[0 ,-1,-1,-1,0,1,1, 1]jj=[-1,-1, 0, 1,1,1,0,-1]for x in range(X):for y in range(Y):w_2=weights[x,y]for k,(i,j) in enumerate(zip(ii,jj)):if(x+i >=0 and x+i<X and y+j>=0 and y+j <Y):w_1=weights[x+i,y+j]um[x,y,k]=np.linalg.norm(w_1-w_2)um=np.nansum(um,axis=2)return um/um.max()#计算量化误差   计算每个样本点和映射点之间的平均距离
def get_quantization_error(data,weights):w_x,w_y=zip(*[get_winner_index(d,weights) for d in datas])error=datas-weights[w_x,w_y]             # 数据域聚类中心的距离error=np.linalg.norm(error,axis=-1)  return np.mean(error)

训练完成后，返回输出节点的weights,维度为$[X,Y,D]$, 相当于固化了模型的权重weights, weights表征了当前的训练样本。

测试

if __name__ == "__main__":# seed 数据展示columns=['area','perimeter','compactness','length_kernel','width_kernel','asymmetry_coefficient','length_kernel_groove','target']data = pd.read_csv('seeds_dataset.txt',names=columns,sep='\t+',engine='python')labs=data['target'].valueslab_names={1:'Kama',2:'Rosa',3:'Canadian'}datas=data[data.columns[:-1]].valuesN,D=np.shape(datas)print(N,D)# 对训练数据进行标准化datas = feature_normalization(datas)#SOM的训练weights=train_SOM()X=9,Y=9,N_epoch=2,datas=datas,sigma=1.5,init_weight_func=weights_PCA)# 获取UMAP  用于可视化UM=get_U_Matrix(weights)plt.figure(figure=(9,9))plt.pcolor(UM.T,cmap='bone_r')  #plotting the distance map as backgroundplt.colorbar()

测试数据

U_Matrix

• 颜色越深说明与邻近点的关系越强烈，颜色越强说明与邻近点的关系越不强烈。

测试分类的效果

```python
if __name__ == "__main__":# seed 数据展示columns=['area','perimeter','compactness','length_kernel','width_kernel','asymmetry_coefficient','length_kernel_groove','target']data = pd.read_csv('seeds_dataset.txt',names=columns,sep='\t+',engine='python')labs=data['target'].valueslab_names={1:'Kama',2:'Rosa',3:'Canadian'}datas=data[data.columns[:-1]].valuesN,D=np.shape(datas)print(N,D)# 对训练数据进行标准化datas = feature_normalization(datas)#SOM的训练weights=train_SOM()X=9,Y=9,N_epoch=2,datas=datas,sigma=1.5,init_weight_func=weights_PCA)# 获取UMAP  用于可视化UM=get_U_Matrix(weights)plt.figure(figure=(9,9))plt.pcolor(UM.T,cmap='bone_r')  #plotting the distance map as backgroundplt.colorbar()# 查看分类的效果markers=['o','s','D']colors =['C0','C1','C2']for i in range(N):x =datas[i]w=get_winner_index(x,weights)i_lab=labs[i]-1plt.plot(w[0]+.5,w[1]+.5,markers[i_lab],markerfacecolor='None'markeredgecolor=colors[i_lab],markersize=12,markeredgewidth=2)plt.show()