t-SNE 原理及Python实例

由于毕业设计有部分工作需要对比两个图像数据集合的差异，为了可视化差异，利用了目前降维首选的t-SNE。我花了点时间看了sklearn里面关于这部分的文档，也查阅了相关博客，最终成功的将两种图片数据集作了一个可视化的对比。我觉得这个方法很硬核而且还蛮有意思的，利用python sklearn库也很容易实现，加上很多教程都是仅仅列出了一个"手写数据集识别"的例子,而我们学习该法的目的在于学习如何将t-SNE利用到自己的特定任务中，本文就记录了关于该方法如何实现两个图片数据集的降维可视化并且与大家作一个分享。初学者，如果讲述有误，恳请指正！文末给出参考博客和官方文档！！

任务陈述：在Images目录下有两种图片，各500张，前500张为第一类，后500张为第二类如下图，可以看到两类的风格不一样，需要想办法将这两类降到2维或者3维聚类可视化。

两类图片数据集

在看下面内容之前，可以先看一下利用t-SEN手写数据集分类的例子,（文末也有实现）这也是各大博客最喜欢贴的例子，或者结合起来看。

直接上代码

确保python安装好了各个依赖库，特别是sklearn

第一部分是获取数据的代码：

无论你是什么类型的数据（图片，音频，文本等等），先将它先处理成np.array类型，确保data的形状为(数据条数，每条数据维度)，label的形状(数据条数,) 方便后面调用t-SNE方法

import os
import numpy as np
import cv2
from time import time
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn import datasets #手写数据集要用到
from sklearn.manifold import TSNE
#该函数是关键，需要根据自己的数据加以修改，将图片存到一个np.array里面，并且制作标签
#因为是两类数据，所以我分别用0,1来表示
def get_data(Input_path): #Input_path为你自己原始数据存储路径，我的路径就是上面的'./Images'Image_names=os.listdir(Input_path) #获取目录下所有图片名称列表data=np.zeros((len(Image_names),40000)) #初始化一个np.array数组用于存数据label=np.zeros((len(Image_names),)) #初始化一个np.array数组用于存数据#为前500个分配标签1，后500分配0for k in range(500):label[k]=1#读取并存储图片数据，原图为rgb三通道，而且大小不一，先灰度化，再resize成200x200固定大小for i in range(len(Image_names)):image_path=os.path.join(Input_path,Image_names[i])img=cv2.imread(image_path)img_gray=cv2.cvtColor(img,cv2.COLOR_RGB2GRAY)img=cv2.resize(img_gray,(200,200))img=img.reshape(1,40000)data[i]=imgn_samples, n_features = data.shapereturn data, label, n_samples, n_features‘’‘下面的两个函数，
一个定义了二维数据，一个定义了3维数据的可视化
不作详解，也无需再修改感兴趣可以了解matplotlib的常见用法
’‘’
def plot_embedding_2D(data, label, title):x_min, x_max = np.min(data, 0), np.max(data, 0)data = (data - x_min) / (x_max - x_min)fig = plt.figure()for i in range(data.shape[0]):plt.text(data[i, 0], data[i, 1], str(label[i]),color=plt.cm.Set1(label[i]),fontdict={'weight': 'bold', 'size': 9})plt.xticks([])plt.yticks([])plt.title(title)return fig
def plot_embedding_3D(data,label,title): x_min, x_max = np.min(data,axis=0), np.max(data,axis=0) data = (data- x_min) / (x_max - x_min) ax = plt.figure().add_subplot(111,projection='3d') for i in range(data.shape[0]): ax.text(data[i, 0], data[i, 1], data[i,2],str(label[i]), color=plt.cm.Set1(label[i]),fontdict={'weight': 'bold', 'size': 9}) return fig#主函数
def main():data, label, n_samples, n_features = get_data('./Images') #根据自己的路径合理更改print('Begining......') #时间会较长，所有处理完毕后给出finished提示tsne_2D = TSNE(n_components=2, init='pca', random_state=0) #调用TSNEresult_2D = tsne_2D.fit_transform(data)tsne_3D = TSNE(n_components=3, init='pca', random_state=0)result_3D = tsne_3D.fit_transform(data)print('Finished......')#调用上面的两个函数进行可视化fig1 = plot_embedding_2D(result_2D, label,'t-SNE')plt.show(fig1)fig2 = plot_embedding_3D(result_3D, label,'t-SNE')plt.show(fig2)
if __name__ == '__main__':main()

下面给出结果图：

降到2d可视化

降到3d的可视化

可以明显看出两类数据的差异巨大。t-SNE的降维能力实在是强悍！上面提到了”手写数据集“的例子，下面也作一个说明：

将上述get_data改成：

def get_data(Input_path)digits = datasets.load_digits(n_class=6)data = digits.datalabel = digits.targetn_samples, n_features = data.shape
return data, label, n_samples, n_features

这样就能实现手写数据集的分类可视化。

手写数据2d

手写数字3d

参考：

73 人赞同了该文章

t-SNE Python 例子

t-Distributed Stochastic Neighbor Embedding (t-SNE)是一种降维技术，用于在二维或三维的低维空间中表示高维数据集，从而使其可视化。与其他降维算法(如PCA)相比，t-SNE创建了一个缩小的特征空间，相似的样本由附近的点建模，不相似的样本由高概率的远点建模。

在高水平上，t-SNE为高维样本构建了一个概率分布，相似的样本被选中的可能性很高，而不同的点被选中的可能性极小。然后，t-SNE为低维嵌入中的点定义了相似的分布。最后，t-SNE最小化了两个分布之间关于嵌入点位置的Kullback-Leibler（KL)散度。

t-SNE是一种集降维与可视化于一体的技术，它是基于SNE可视化的改进，解决了SNE在可视化后样本分布拥挤、边界不明显的特点，是目前较好的降维可视化手段。

算法

如前所述，t-SNE采用一个高维数据集，并将其简化为一个保留了大量原始信息的低维图。

假设我们有一个由3个不同的类组成的数据集。

我们希望将2D地块缩减为1D地块，同时保持集群之间清晰的边界。

回想一下，简单地将数据投射到一个轴上是一种降低维数的糟糕方法，因为我们会丢失大量的信息。

相反，我们可以使用降维技术(提示:t-SNE)来实现我们想要的。t-SNE算法的第一步是测量一个点相对于其他点的距离。我们不是直接处理这些距离，而是将它们映射到一个概率分布。

在分布中，相对于当前点距离最小的点有很高的可能性，而远离当前点的点有很低的可能性。

再看一下2D图，请注意蓝色的点团比绿色的点团更分散。如果我们不解决比例上的差异，绿色点的可能性将大于蓝色点的可能性。为了解释这一事实，我们除以概率的总和。

因此，尽管两点之间的绝对距离不同，但它们被认为是相似的。

让我们试着把这些概念与基本的理论联系起来。数学上，正态分布的方程如下。

如果我们把指数前面的所有东西都放下，用另一个点代替均值，同时解决前面讨论的规模问题，我们得到了方程：

Visualizing Data using t-SNE

接下来，让我们讨论一下如何使用缩小的特征空间。首先，我们创建一个n_samples x n_components矩阵(在本例中为9x1)，并用随机值(即位置)填充它。

如果我们采用与上面相似的方法(测量点之间的距离并将它们映射到一个概率分布)，我们得到以下等式。

注意,像以前一样,我们把一个正态分布的方程,把面前的一切,用另一个点,而不是意味着,占规模的总和除以所有其他点的可能性(不要问我为什么我们摆脱了标准差)。如果我们能使简化特征空间中的点的概率分布近似于原始特征空间中的点，我们就能得到定义良好的聚类。

为了做到这一点，我们使用了一种叫做KL散度的东西。KL散度是衡量一个概率分布与另一个概率分布的差异。

KL散度值越低，两个分布越接近。KL散度为0意味着这两个分布是相同的。

这应该有望带来大量的想法。回想一下，在线性回归的情况下，我们是如何通过使用梯度下降最小化成本函数(即均方误差)来确定最佳拟合线的。在t-SNE中，我们使用梯度下降最小化所有数据点上的Kullback-Leiber差异的总和。

我们对成本函数对每一点求偏导数，以便得到每次更新的方向。

Python代码

很多时候，我们在使用一些库时，并没有真正理解其中的含义。在这一节中，我将尝试以Python代码的形式实现算法和相关的数学方程。为了帮助完成这个过程，我从scikiti-learn库中提取了一些TSNE的源代码。

首先，我们将导入以下库并设置一些属性，这些属性将在我们绘制数据时发挥作用。

import numpy as np
from sklearn.datasets import load_digits
from scipy.spatial.distance import pdist
from sklearn.manifold.t_sne import _joint_probabilities
from scipy import linalg
from sklearn.metrics import pairwise_distances
from scipy.spatial.distance import squareform
from sklearn.manifold import TSNE
from matplotlib import pyplot as plt
import seaborn as sns
sns.set(rc={'figure.figsize':(11.7,8.27)})
palette = sns.color_palette("bright", 10)

在这个例子中，我们将使用手绘数字。scikiti -learn库提供了一种将它们导入程序的方法。

X, y = load_digits(return_X_y=True)

我们要选择2或3作为分量的数量，因为t-SNE是严格用于可视化的我们只能看到三维的东西。另一方面，复杂度与算法中使用的最近邻的数量有关。不同的复杂度可能导致最终结果的剧烈变化。在本例中，我们将其设置为t-SNE(30)的scitkit-learn实现的默认值。根据numpy文档，机器的epsilon是最小的可表示的正数，因此1.0 + eps != 1.0。换句话说，任何低于机器epsilon的数字都不能被计算机操作，因为它缺少必要的比特。我们会看到，代码贡献者使用的是np.maximum来检查矩阵中的值是否小于机器的值，并在它们小于时替换它们。

接下来，我们定义（fit)拟合函数。我们将在转换数据时调用fit函数。

def fit(X):n_samples = X.shape[0]# Compute euclidean distancedistances = pairwise_distances(X, metric='euclidean', squared=True)# Compute joint probabilities p_ij from distances.P = _joint_probabilities(distances=distances, desired_perplexity=perplexity, verbose=False)# The embedding is initialized with iid samples from Gaussians with standard deviation 1e-4.X_embedded = 1e-4 * np.random.mtrand._rand.randn(n_samples, n_components).astype(np.float32)# degrees_of_freedom = n_components - 1 comes from# "Learning a Parametric Embedding by Preserving Local Structure"# Laurens van der Maaten, 2009.degrees_of_freedom = max(n_components - 1, 1)return _tsne(P, degrees_of_freedom, n_samples, X_embedded=X_embedded)

这个函数中有很多东西我们一步一步地把它分解。

我们将样本数量存储在一个变量中，以备将来使用。
我们计算每个数据点之间的欧式距离。这个对应于前一个方程中的。

我们将前一步计算的欧氏距离作为参数传递给_join_probability函数，然后该函数计算并返回一个p_ji值矩阵(使用相同的方程)。
我们使用从标准偏差1e-4的高斯分布中随机选取值来创建减少的特征空间。
我们定义自由度。源代码中有一个注释，告诉您查看这篇解释其推理的论文。基本上，经验表明，当我们使用组件的数量减去1时，我们会得到更好的结果(粗体部分)。

最后调用tsne函数，其实现如下。

def _tsne(P, degrees_of_freedom, n_samples, X_embedded):params = X_embedded.ravel()obj_func = _kl_divergenceparams = _gradient_descent(obj_func, params, [P, degrees_of_freedom, n_samples, n_components])X_embedded = params.reshape(n_samples, n_components)
return X_embedded

这个函数并没有太多的变化。首先，我们使用np.ravel拉平我们的向量到一个一维数组。

然后利用梯度下降最小化KL散度。一旦它完成，我们把嵌入改回一个2D数组并返回它。

接下来，让我们来看看有更多符合它的东西。下面的代码块负责计算KL散度和梯度形式的误差。

def _kl_divergence(params, P, degrees_of_freedom, n_samples, n_components):X_embedded = params.reshape(n_samples, n_components)dist = pdist(X_embedded, "sqeuclidean")dist /= degrees_of_freedomdist += 1.dist **= (degrees_of_freedom + 1.0) / -2.0Q = np.maximum(dist / (2.0 * np.sum(dist)), MACHINE_EPSILON)# Kullback-Leibler divergence of P and Qkl_divergence = 2.0 * np.dot(P, np.log(np.maximum(P, MACHINE_EPSILON) / Q))# Gradient: dC/dYgrad = np.ndarray((n_samples, n_components), dtype=params.dtype)PQd = squareform((P - Q) * dist)for i in range(n_samples):grad[i] = np.dot(np.ravel(PQd[i], order='K'),X_embedded[i] - X_embedded)grad = grad.ravel()c = 2.0 * (degrees_of_freedom + 1.0) / degrees_of_freedomgrad *= c
return kl_divergence, grad

同样，让我们一步一步地遍历代码。

第一部分计算低维图中点的概率分布。

作者实际上使用了上述方程的一个变量，其中包括自由度。

代表学生-t分布的自由度

我们计算KL散度。

我们计算梯度(偏导).dist实际上是yi - yj in:

同样地，他们使用了上面的方程的变化来表示自由度。

其中代表学生-t分布的自由度

梯度下降函数通过最小化KL散度来更新嵌入中的值。当梯度范数低于阈值，或者当我们达到最大迭代次数而没有取得任何进展时，我们会提前停止。

def _gradient_descent(obj_func, p0, args, it=0, n_iter=1000,n_iter_check=1, n_iter_without_progress=300,momentum=0.8, learning_rate=200.0, min_gain=0.01,min_grad_norm=1e-7):p = p0.copy().ravel()update = np.zeros_like(p)gains = np.ones_like(p)error = np.finfo(np.float).maxbest_error = np.finfo(np.float).maxbest_iter = i = itfor i in range(it, n_iter):
error, grad = obj_func(p, *args)
grad_norm = linalg.norm(grad)
inc = update * grad < 0.0dec = np.invert(inc)gains[inc] += 0.2gains[dec] *= 0.8np.clip(gains, min_gain, np.inf, out=gains)grad *= gainsupdate = momentum * update - learning_rate * gradp += update
print("[t-SNE] Iteration %d: error = %.7f,"" gradient norm = %.7f"% (i + 1, error, grad_norm))if error < best_error:best_error = errorbest_iter = ielif i - best_iter > n_iter_without_progress:breakif grad_norm <= min_grad_norm:break
return p

接下来，用数据来调用fit函数。

X_embedded = fit(X)

正如我们所看到的，该模型根据像素的位置将不同的数字进行了相当不错的分离。

sns.scatterplot(X_embedded[:,0], X_embedded[:,1], hue=y, legend='full', palette=palette)

让我们用scikit-learn实现t-SNE做同样的事情。

tsne = TSNE()
X_embedded = tsne.fit_transform(X)

正如我们所看到的，该模型成功地把一个64维数据集投影到一个2维空间中，从而使类似的样本聚类在一起。

sns.scatterplot(X_embedded[:,0], X_embedded[:,1], hue=y, legend='full', palette=palette)

结语

t-SNE是目前来说效果最好的数据降维与可视化方法，但是它的缺点也很明显，比如：

占内存大，运行时间长。
专用于可视化，即嵌入空间只能是2维或3维。
需要尝试不同的初始化点，以防止局部次优解的影响。

但是，当我们想要对高维数据进行分类，又不清楚这个数据集有没有很好的可分性（即同类之间间隔小，异类之间间隔大），可以通过t-SNE投影到2维或者3维的空间中观察一下。如果在低维空间中具有可分性，则数据是可分的；如果在高维空间中不具有可分性，可能是数据不可分，也可能仅仅是因为不能投影到低维空间。参考文献：

[1] https://towardsdatascience.com/t-sne-python-example-1ded9953f26

[2] Visualizing Data using t-SNE