t-SNE 原理及Python实例
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图，请注意蓝色的点团比绿色的点团更分散。如果我们不解决比例上的差异，绿色点的可能性将大于蓝色点的可能性。为了解释这一事实，我们除以概率的总和。
因此，尽管两点之间的绝对距离不同，但它们被认为是相似的。
让我们试着把这些概念与基本的理论联系起来。数学上，正态分布的方程如下。
P(x)=\frac{1}{\sigma \sqrt{2 \pi}} e^{-(x-\mu)^{2} /\left(2 \sigma^{2}\right)} \\
如果我们把指数前面的所有东西都放下，用另一个点代替均值，同时解决前面讨论的规模问题，我们得到了方程：
Visualizing Data using t-SNE
p_{j | i}=\frac{\exp \left(-\left\|x_{i}-x_{j}\right\|^{2} / 2 \sigma_{i}^{2}\right)}{\sum_{k \neq i} \exp \left(-\left\|x_{i}-x_{k}\right\|^{2} / 2 \sigma_{i}^{2}\right)} \\
接下来，让我们讨论一下如何使用缩小的特征空间。首先，我们创建一个n_samples x n_components矩阵(在本例中为9x1)，并用随机值(即位置)填充它。
如果我们采用与上面相似的方法(测量点之间的距离并将它们映射到一个概率分布)，我们得到以下等式。
q_{j | i}=\frac{\exp \left(-\left\|y_{i}-y_{j}\right\|^{2}\right)}{\sum_{k \neq i} \exp \left(-\left\|y_{i}-y_{k}\right\|^{2}\right)} \\
注意,像以前一样,我们把一个正态分布的方程,把面前的一切,用另一个点,而不是意味着,占规模的总和除以所有其他点的可能性(不要问我为什么我们摆脱了标准差)。 如果我们能使简化特征空间中的点的概率分布近似于原始特征空间中的点，我们就能得到定义良好的聚类。
为了做到这一点，我们使用了一种叫做KL散度的东西。KL散度是衡量一个概率分布与另一个概率分布的差异。
KL散度值越低，两个分布越接近。KL散度为0意味着这两个分布是相同的。
这应该有望带来大量的想法。回想一下，在线性回归的情况下，我们是如何通过使用梯度下降最小化成本函数(即均方误差)来确定最佳拟合线的。在t-SNE中，我们使用梯度下降最小化所有数据点上的Kullback-Leiber差异的总和。
C=\sum_{i} K L\left(P_{i} \| Q_{i}\right)=\sum_{i} \sum_{j} p_{j | i} \log \frac{p_{j | i}}{q_{j | i}} \\
我们对成本函数对每一点求偏导数，以便得到每次更新的方向。
\frac{\delta C}{\delta y_{i}}=2 \sum_{j}\left(p_{j | i}-q_{j | i}+p_{i | j}-q_{i | j}\right)\left(y_{i}-y_{j}\right) \\
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 distance
distances = 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)这个函数中有很多东西我们一步一步地把它分解。
- 我们将样本数量存储在一个变量中，以备将来使用。
- 我们计算每个数据点之间的欧式距离。这个对应于前一个方程中的||xi - xj||^2。
p_{j | i}=\frac{\exp \left(-\left\|\mathbf{x}_{i}-\mathbf{x}_{j}\right\|^{2} / 2 \sigma_{i}^{2}\right)}{\sum_{k \neq i} \exp \left(-\left\|\mathbf{x}_{i}-\mathbf{x}_{k}\right\|^{2} / 2 \sigma_{i}^{2}\right)} \\
- 我们将前一步计算的欧氏距离作为参数传递给
_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_divergence
params = _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_freedom
dist += 1.
dist **= (degrees_of_freedom + 1.0) / -2.0
Q = np.maximum(dist / (2.0 * np.sum(dist)), MACHINE_EPSILON)
# Kullback-Leibler divergence of P and Q
kl_divergence = 2.0 * np.dot(P, np.log(np.maximum(P, MACHINE_EPSILON) / Q))
# Gradient: dC/dY
grad = 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_freedom
grad *= c
return kl_divergence, grad同样，让我们一步一步地遍历代码。
- 第一部分计算低维图中点的概率分布。
q_{i j}=\frac{\left(1+\left\|\mathbf{y}_{i}-\mathbf{y}_{j}\right\|^{2}\right)^{-1}}{\sum_{k \neq l}\left(1+\left\|\mathbf{y}_{k}-\mathbf{y}_{l}\right\|^{2}\right)^{-1}} \\
作者实际上使用了上述方程的一个变量，其中包括自由度。
代表学生-t分布的自由度
- 我们计算KL散度。
K L(P \| Q)=\sum_{i, j} p_{i j} \log \frac{p_{i j}}{q_{i j}} \\
- 我们计算梯度(偏导).dist实际上是yi - yj in:
\frac{\delta C}{\delta y_{i}}=2 \sum_{j}\left(p_{j | i}-q_{j | i}+p_{i | j}-q_{i | j}\right)\left(y_{i}-y_{j}\right) \\
同样地，他们使用了上面的方程的变化来表示自由度。
\frac{\delta C}{\delta f\left(x_{i} | W\right)}=\frac{2 \alpha+2}{\alpha} \sum_{j}\left(p_{i j}-q_{i j}\right)\left(f\left(x_{i} | W\right)-f\left(x_{j} | W\right)\right) \\
其中\alpha代表学生-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).max
best_error = np.finfo(np.float).max
best_iter = i = it
for i in range(it, n_iter):
error, grad = obj_func(p, *args)
grad_norm = linalg.norm(grad)
inc = update * grad < 0.0
dec = np.invert(inc)
gains[inc] += 0.2
gains[dec] *= 0.8
np.clip(gains, min_gain, np.inf, out=gains)
grad *= gains
update = momentum * update - learning_rate * grad
p += update
print("[t-SNE] Iteration %d: error = %.7f,"
" gradient norm = %.7f"
% (i + 1, error, grad_norm))
if error < best_error:
best_error = error
best_iter = i
elif i - best_iter > n_iter_without_progress:
break
if 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
