基于SOM算法的Iris数据分类
自组织特征映射神经网络SOM(Self-Organizing Feature Map)是一种无监督学习算法，不同于一般神经网络基于损失函数的优化训练，SOM是运用竞争学习策略来逐步优化网络的。SOM算法作为一种优良的聚类工具，具有无需监督，能自动对输入模式进行聚类的优点，目前已经得到了广泛的应用。本文利用SOM算法，最终实现了对Iris数据集的分类。
一、 算法介绍
1.1 SOM算法
自组织映射神经网络(self-organizing map，SOM)是20世纪80年代 Kohonen教授在研究联想记忆和自适应学习机器的基础上发展起来的，早期叫做自组织特征映射(self-organizing feature map即SOFM)。后来，人们发现这种学习机器和大脑皮层上的分区组织现象有很多相似性，所以也经常从对大脑功能数学模拟的角度来讨论SOM网络。
SOM是一种无监督的人工神经网络。不同于一般神经网络基于损失函数的反向传递来训练，它运用竞争学习(competitive learning)策略,依靠神经元之间互相竞争逐步优化网络，且使用近邻关系函数(neighborhood function)来维持输入空间的拓扑结构。维持输入空间的拓扑结构：意味着二维映射包含了数据点之间的相对距离。输入空间中相邻的样本会被映射到相邻的输出神经元。
由于SOM是基于无监督学习的,这就意味着训练阶段将不需要人工介入(即不需要样本标签)，我们可以在不知道类别的情况下，对数据进行聚类；可以识别针对某问题具有内在关联的特征。
1.2 算法特点
SOM作为一种广泛应用的算法，具有着其自独特的优势：
（1） 神经网络，竞争学习策略
（2） 无监督学习，不需要额外标签
（3） 适合高维数据的可视化，能够维持输入空间的拓扑结构
（4） 具有很高的泛化能力，甚至能识别之前从没遇过的输入样本
二、 相关理论概述
2.1 SOM网络结构
SOM网络结构比较简单，只有输入层和输出层（输出层我们通常也称为竞争层），输入层神经元的数量是由输入向量的维度决定的，一个神经元对应一个特征，SOM网络结构的区别主要在竞争层：可以的有一维和二维的（竞争层也可以有更高的维度。不过出于可视化的目的，高维竞争层用的比较少）。其中，二维平面有2种平面结构：①Rectangular；②Hexagonal
在竞争层，神经元可以排列形成多种形式。每个样本都可以映射到与其距离最小的竞争层神经元上，最终在竞争层生成一个低维、离散的映射。竞争层SOM神经元的数量决定了最终模型的粒度与规模，这对最终模型的准确性与泛化能力影响很大。结构如下图所示：
与前馈型神经网路不同，SOM网络的神经元节点都在同一层上，在一个平面上呈规则排列。常见的排列形式包括方形网格排列或蜂窝状排列。样本特征向量的每一维都通过一定的权重输入到SOM网络的每一个节点上。
自组织映射网络的神经元节点之间并没有直接的连接，但是，在神经元平面上相邻的节点间在学习(训练)过程中有一定的相互影响，构成领域相互作用。
神经元节点的计算功能就是对输入的样本给出响应。输入向量连接到某个节点的权值组成的向量称作该节点的权值向量。一个节点对输入样本的响应强度，就是该节点的权值向量与输入向量的匹配程度，可以用欧式距离或内积来计算，如果距离小或内积大则响应强度大。对一个输入样本，在神经元平面上所有的节点中响应最大的节点称作获胜节点(winner)。
最终，可以将SOM算法的结构归纳为：①输入层：接收外界信息，将输入模式向竞争层传递，起“观察”作用；②竞争层：负责对输入模式进行“分析比较”，寻找规律并归类。SOM完整网络结构如下图所示：
2.2 算法原理
自组织映射网络的学习过程比较简单，基本算法如下：
设X={x∈R^d}是d维样本向量集合，记所有神经元集合为A，第i个每个神经元的权重为m_i：
（1）权值初始化：用小随机数初始化权重向量。注意各个节点的初始权值不能相等；
（2）在时刻t，按照给定的顺序或随机顺序加入一个样本。记为x(t)；
（3）计算神经元响应，找到当前获胜节点c。如用欧氏距离作为匹配准则，则获胜节点为：
（4）权值竞争学习。对所有神经元节点，用下述准则更新各自的权值：
其中，α(t)是学习的步长，d[.,.]是两个向量间的欧式距离，h_ci (t)是节点i与c之间的近邻函数值，如果采用方形网格结构，则相当于在节点c的周围定义一个矩形领域范围N_c (t)，在该领域内h_ci (t)为1、否则为0，即权值按照以下规则更新：
（5）更新步长α(t)和领域N_c (t)，如达到终止条件，则算法停止；否则置t=t+1，继续步骤（2）
因此，总的来说竞争学习的步骤是：①向量归一化②寻找获胜神经元③网络输出与权值调整步骤④完成后回到步骤①继续训练，直到学习率衰减到0。其中，学习率αϵ(0,1]，一般随着学习的进展而减小，即调整的程度越来越小，神经元（权重）趋于聚类中心。
2.3 理论推导
（1）寻找获胜神经元
网络的输出神经元之间相互竞争以求被激活，结果在每一时刻只有一个输出神经元被激活。这个被激活的神经元称为竞争获胜神经元，而其它神经元的状态被抑制，故称为Winner Take All。
首先，对网络当前输入模式向量X和竞争层中各神经元对应的权重向量W_j全部进行归一化，使得X和W_j模为1。当网络得到一个输入模式向量X时，竞争层的所有神经元对应的权重向量均与其进行相似性比较，并将最相似的权重向量判为竞争获胜神经元。前面刚说过，归一化后，相似度最大就是内积最大：
也就是在单位圆（2D情况）中找到夹角最小的点，即：
（2）neighborhood function
neighborhood函数用来确定优胜节点对其近邻节点的影响强弱，即优胜邻域中每个节点的更新幅度。最常见的选择是高斯函数，它可以表征优胜邻域内，影响强弱与距离的关系。
其中，winner是优胜节点在输出平面的坐标；sigma确定邻域范围，且sigma越大，邻域范围越大。但是sigma必须大于0，否则没有神经元会被更新，且sigma不能大于2维输出平面的边长。
（3）Bubble近邻函数
Bubble函数能很好的在计算成本与高斯核的准确性之间平衡，我们可以用bubble函数很好地去近似估计高斯，bubble在优胜节点的邻域范围内是个常数，因此，邻域内的所有神经元更新的幅度是相同的。
三、 实验结果
我们本次实验是对Iris数据集利用SOM算法进行分类，下图为SOM算法在Iris数据集上的分类情况：
为了能够更加清晰的观察到SOM算法对Iris数据集的分类效果，我们对最后的输出结果进行了图层可视化。
首先，我们利用了U-Matrix矩阵来对结果就行可视化。根据权重矩阵W,我们可以计算每个神经元距离它的邻近神经元们的距离，计算好的矩阵就是U-Matrix矩阵。这样我们就能够看到训练样本都分别落在输出层的哪些格子和它们对应的标签信息。可视化结果如下图所示：
我们可以看到，三种类别的样本落在了输出平面上的不同位置，并且这条分界线大致将蓝色样本与橙色、绿色样本划分开了。
虽然上述的可视化图案能让我们看到样本的大致分布，并且具有较好视觉效果。但仍有两个不足之处：①不清楚每个格子的样本数量具体是多少②如果有个类别落到同一个格子，看不出该类别的比例。
因此，我们对上述可视化进行了改进，我们决定在每个格子里画饼图，且用颜色表示类别，用数字表示总样本数量，这样可以更为清楚的可视化Iris数据集的分类结果。改进的可视化结果如下图所示：
四、 总结
Iris数据集比较简单，我们的分类器基本上可以在测试集上达得96的F1-score。从上述实验结果看来，它将相邻关系强加在簇质心上，所以，互为邻居的簇之间比非邻居的簇之间更相关。这种联系有利于聚类结果的解释和可视化。最后，由于SOM算法提出的比较早，因此我个人认为它还是具有很大的发展空间的。
五、 代码实现
GitHub：https://github.com/DeepVegChicken/Learning-Som_IrisClassification
运行环境：
python 3.7.9
numpy 1.19.1
minisom 2.2.7
scikit-learn 0.23.2
matplotlib 3.3.1
SOM算法的单独实现：
import numpy as np
import pylab as pl
class SOM(object):
def __init__(self, X, output, iteration, batch_size):
"""
:param X: 形状是N*D， 输入样本有N个,每个D维
:param output: (n,m)一个元组，为输出层的形状是一个n*m的二维矩阵
:param iteration:迭代次数
:param batch_size:每次迭代时的样本数量
初始化一个权值矩阵，形状为D*(n*m)，即有n*m权值向量，每个D维
"""
self.X = X
self.output = output
self.iteration = iteration
self.batch_size = batch_size
self.W = np.random.rand(X.shape[1], output[0] * output[1])
print(self.W.shape)
def GetN(self, t):
"""
:param t:时间t, 这里用迭代次数来表示时间
:return: 返回一个整数，表示拓扑距离，时间越大，拓扑邻域越小
"""
a = min(self.output)
return int(a - float(a) * t / self.iteration)
def Geteta(self, t, n):
"""
:param t: 时间t, 这里用迭代次数来表示时间
:param n: 拓扑距离
:return: 返回学习率，
"""
return np.power(np.e, -n) / (t + 2)
def updata_W(self, X, t, winner):
N = self.GetN(t)
for x, i in enumerate(winner):
to_update = self.getneighbor(i[0], N)
for j in range(N + 1):
e = self.Geteta(t, j)
for w in to_update[j]:
self.W[:, w] = np.add(self.W[:, w], e * (X[x, :] - self.W[:, w]))
def getneighbor(self, index, N):
"""
:param index:获胜神经元的下标
:param N: 邻域半径
:return ans: 返回一个集合列表，分别是不同邻域半径内需要更新的神经元坐标
"""
a, b = self.output
length = a * b
def distence(index1, index2):
i1_a, i1_b = index1 // a, index1 % b
i2_a, i2_b = index2 // a, index2 % b
return np.abs(i1_a - i2_a), np.abs(i1_b - i2_b)
ans = [set() for i in range(N + 1)]
for i in range(length):
dist_a, dist_b = distence(i, index)
if dist_a <= N and dist_b <= N: ans[max(dist_a, dist_b)].add(i)
return ans
def train(self):
"""
train_Y:训练样本与形状为batch_size*(n*m)
winner:一个一维向量，batch_size个获胜神经元的下标
:return:返回值是调整后的W
"""
count = 0
while self.iteration > count:
train_X = self.X[np.random.choice(self.X.shape[0], self.batch_size)]
normal_W(self.W)
normal_X(train_X)
train_Y = train_X.dot(self.W)
winner = np.argmax(train_Y, axis=1).tolist()
self.updata_W(train_X, count, winner)
count += 1
return self.W
def train_result(self):
normal_X(self.X)
train_Y = self.X.dot(self.W)
winner = np.argmax(train_Y, axis=1).tolist()
print(winner)
return winner
def normal_X(X):
"""
:param X:二维矩阵，N*D，N个D维的数据
:return: 将X归一化的结果
"""
N, D = X.shape
for i in range(N):
temp = np.sum(np.multiply(X[i], X[i]))
X[i] /= np.sqrt(temp)
return X
def normal_W(W):
"""
:param W:二维矩阵，D*(n*m)，D个n*m维的数据
:return: 将W归一化的结果
"""
for i in range(W.shape[1]):
temp = np.sum(np.multiply(W[:, i], W[:, i]))
W[:, i] /= np.sqrt(temp)
return W
# 画图
def draw(C):
colValue = ['r', 'y', 'g', 'b', 'c', 'k', 'm']
for i in range(len(C)):
coo_X = [] # x坐标列表
coo_Y = [] # y坐标列表
for j in range(len(C[i])):
coo_X.append(C[i][j][0])
coo_Y.append(C[i][j][1])
pl.scatter(coo_X, coo_Y, marker='x', color=colValue[i % len(colValue)], label=i)
pl.legend(loc='upper right')
pl.show()
# 数据集：每三个是一组分别是西瓜的编号，密度，含糖量
data = """
1,0.697,0.46,2,0.774,0.376,3,0.634,0.264,4,0.608,0.318,5,0.556,0.215,
6,0.403,0.237,7,0.481,0.149,8,0.437,0.211,9,0.666,0.091,10,0.243,0.267,
11,0.245,0.057,12,0.343,0.099,13,0.639,0.161,14,0.657,0.198,15,0.36,0.37,
16,0.593,0.042,17,0.719,0.103,18,0.359,0.188,19,0.339,0.241,20,0.282,0.257,
21,0.748,0.232,22,0.714,0.346,23,0.483,0.312,24,0.478,0.437,25,0.525,0.369,
26,0.751,0.489,27,0.532,0.472,28,0.473,0.376,29,0.725,0.445,30,0.446,0.459"""
a = data.split(',')
dataset = np.mat([[float(a[i]), float(a[i + 1])] for i in range(1, len(a) - 1, 3)])
dataset_old = dataset.copy()
som = SOM(dataset, (5, 5), 1, 30)
som.train()
res = som.train_result()
classify = {}
for i, win in enumerate(res):
if not classify.get(win[0]):
classify.setdefault(win[0], [i])
else:
classify[win[0]].append(i)
C = [] # 未归一化的数据分类结果
D = [] # 归一化的数据分类结果
for i in classify.values():
C.append(dataset_old[i].tolist())
D.append(dataset[i].tolist())
draw(C)
draw(D)
基于SOM的Iris数据集分类：
import math
import numpy as np
from minisom import MiniSom
from sklearn import datasets
from numpy import sum as npsum
from sklearn.metrics import classification_report
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from matplotlib.gridspec import GridSpec
# 分类函数
def classify(som,data,winmap):
default_class = npsum(list(winmap.values())).most_common()[0][0]
result = []
for d in data:
win_position = som.winner(d)
if win_position in winmap:
result.append(winmap[win_position].most_common()[0][0])
else:
result.append(default_class)
return result
# 可视化
def show(som):
""" 在输出层画标签图案 """
plt.figure(figsize=(16, 16))
# 定义不同标签的图案标记
markers = ['o', 's', 'D']
colors = ['C0', 'C1', 'C2']
category_color = {'setosa': 'C0', 'versicolor': 'C1', 'virginica': 'C2'}
# 背景上画U-Matrix
heatmap = som.distance_map()
# 画背景图
plt.pcolor(heatmap, cmap='bone_r')
for cnt, xx in enumerate(X_train):
w = som.winner(xx)
# 在样本Heat的地方画上标记
plt.plot(w[0] + .5, w[1] + .5, markers[Y_train[cnt]], markerfacecolor='None',
markeredgecolor=colors[Y_train[cnt]], markersize=12, markeredgewidth=2)
plt.axis([0, size, 0, size])
ax = plt.gca()
# 颠倒y轴方向
ax.invert_yaxis()
legend_elements = [Patch(facecolor=clr, edgecolor='w', label=l) for l, clr in category_color.items()]
plt.legend(handles=legend_elements, loc='center left', bbox_to_anchor=(1, .95))
plt.show()
""" 在每个格子里画饼图，且用颜色表示类别，用数字表示总样本数量 """
plt.figure(figsize=(16, 16))
the_grid = GridSpec(size, size)
for position in winmap.keys():
label_fracs = [winmap[position][label] for label in [0, 1, 2]]
plt.subplot(the_grid[position[1], position[0]], aspect=1)
patches, texts = plt.pie(label_fracs)
plt.text(position[0] / 100, position[1] / 100, str(len(list(winmap[position].elements()))),
color='black', fontdict={'weight': 'bold', 'size': 15}, va='center', ha='center')
plt.legend(patches, class_names, loc='center right', bbox_to_anchor=(-1, 9), ncol=3)
plt.show()
if __name__ == '__main__':
# 导入数据集
iris = datasets.load_iris()
# 提取iris数据集的标签与数据
feature_names = iris.feature_names
class_names = iris.target_names
X = iris.data
Y = iris.target
# 划分训练集、测试集 7:3
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.3, random_state=0)
# 样本数量
N = X_train.shape[0]
# 维度/特征数量
M = X_train.shape[1]
# 最大迭代次数
max_iter = 200
# 经验公式：决定输出层尺寸
size = math.ceil(np.sqrt(5 * np.sqrt(N)))
print("训练样本个数:{} 测试样本个数:{}".format(N, X_test.shape[0]))
print("输出网格最佳边长为:", size)
# 初始化模型
som = MiniSom(size, size, M, sigma=3, learning_rate=0.5, neighborhood_function='bubble')
# 初始化权值
som.pca_weights_init(X_train)
# 模型训练
som.train_batch(X_train, max_iter, verbose=False)
# 利用标签信息，标注训练好的som网络
winmap = som.labels_map(X_train, Y_train)
# 进行分类预测
y_pred = classify(som, X_test, winmap)
# 展示在测试集上的效果
print(classification_report(Y_test, np.array(y_pred)))
# 可视化
show(som)