理解TSNE算法

结合论文公式与几个python实现理解t-SNE算法。

t-SNE 是一种数据可视化工具,它可以将高维数据降维到2-3维以用于画图,局部相似性被这种embedding所保留。

t-SNE把原空间的数据点之间的距离转换为高斯分布概率,如果两点在高维空间距离越近,那么这个概率值越大。注意到高斯分布的这个标准差$\sigma_i$ 对每个点都是不同的,这也是算法的创新点之一,因为理论上空间不同位置的点的密度是不同的,条件概率如此计算,

$$p_{j|i} = \frac{\exp{(-d(\boldsymbol{x}_i, \boldsymbol{x}_j) / (2 \sigma_i^2)})}{\sum_{i \neq k} \exp{(-d(\boldsymbol{x}_i, \boldsymbol{x}_k) / (2 \sigma_i^2)})}, \quad p_{i|i} = 0,$$
Screen Shot 2018-01-30 at 15.18.45

图中公式是理论方式,实际是先计算条件概率再用下面公式来产生联合分布,

$$p_{ij} = \frac{p_{j|i} + p_{i|j}}{2N}.$$

其中 $\sigma_i$ 将自动确定。这个过程可以通过设置算法的困惑性来影响。

用一个长尾分布(Student-t Distribution,简称为t分布)来表示 embed空间的相似性
$$q_{ij} = \frac{(1 + ||\boldsymbol{y}_i - \boldsymbol{y}_j)||^2)^{-1}}{\sum_{k \neq l} (1 + ||\boldsymbol{y}_k - \boldsymbol{y}_l)||^2)^{-1}},$$
Screen Shot 2018-01-30 at 15.28.56

损失函数是两个分布之间的 Kullback-Leibler divergence(KL散度)

$$KL(P|Q) = \sum_{i \neq j} p_{ij} \log \frac{p_{ij}}{q_{ij}}$$

而为什么说tsne保留的是局部相似性呢?我们从KL散度的公式出发来解释,
Screen Shot 2018-01-30 at 15.33.19
可以看到,当$p_{ij}$很大而$q_{ij}$很小(高维空间距离近,低维空间距离远)惩罚很大,反之惩罚小(高维空间距离远,低维空间距离近)。

而为什么高维空间用高斯分布,低维空间用Student-t Distribution呢?

Screen Shot 2018-01-30 at 15.41.32Screen Shot 2018-01-30 at 15.41.43
原因就是因为降维是必然要带来信息损失,我们要保存局部信息那么必然要损失全局信息,比如我们要把上面的这个2维空间的三个成直角边的点降维到1维,那么把它们放平就保存了局部信息(左中和中右之间的距离保持不变),但是牺牲了全局信息(左右之间的距离变大了)。而Student-t Distribution就能放大这种密度,如下图(tsne默认t分布自由度为1),t分布相比高斯分布更加长尾。
Screen Shot 2018-01-30 at 15.48.47
梯度计算时有优化技巧,如果按下图中的原公式计算,复杂度为$O(N^2)$ Barnes-Hut 树方法就可以优化到$ O(NlogN)$
Screen Shot 2018-01-30 at 15.59.02
Screen Shot 2018-01-30 at 16.01.17
原理类似于用上图中ABC三点中心的距离乘以三来代替计算三者各自的距离。
那么把用barnes树结构来进行深度优先搜索,分别判断其距离是否大于阈值,分块计算距离,这样复杂度就降低了。
Screen Shot 2018-01-30 at 16.01.38

以下是计算损失KL散度的公式,

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
def _kl_divergence(params, P, degrees_of_freedom, n_samples, n_components,
skip_num_points=0):
"""t-SNE objective function: gradient of the KL divergence
of p_ijs and q_ijs and the absolute error."""
X_embedded = params.reshape(n_samples, n_components)
# Q is a heavy-tailed distribution: Student's t-distribution
n = pdist(X_embedded, "sqeuclidean")
n += 1.
n /= degrees_of_freedom
n **= (degrees_of_freedom + 1.0) / -2.0
Q = np.maximum(n / (2.0 * np.sum(n)), MACHINE_EPSILON)
# Optimization trick below: np.dot(x, y) is faster than
# np.sum(x * y) because it calls BLAS
# Objective: C (Kullback-Leibler divergence of P and Q)
kl_divergence = 2.0 * np.dot(P, np.log(P / Q))
# Gradient: dC/dY
grad = np.ndarray((n_samples, n_components))
PQd = squareform((P - Q) * n)
for i in range(skip_num_points, n_samples):
np.dot(_ravel(PQd[i]), X_embedded[i] - X_embedded, out=grad[i])
grad = grad.ravel()
c = 2.0 * (degrees_of_freedom + 1.0) / degrees_of_freedom
grad *= c

用梯度下降(和一些tricks)优化,值得注意的是损失函数非对称,并且不同的训练会导致结果的不同。

sklearn里对于binary search计算 联合分布下面的(_utils._binary_search_perplexity)和Barnes-Hut 树计算梯度(_barnes_hut_tsne.gradient)都是C实现,有空再来研究。

1
2
3
4
5
6
7
8
9
10
11
12
def _joint_probabilities(distances, desired_perplexity, verbose):
"""Compute joint probabilities p_ij from distances."""
# Compute conditional probabilities such that they approximately match
# the desired perplexity
distances = astype(distances, np.float32, copy=False)
conditional_P = _utils._binary_search_perplexity(
distances, None, desired_perplexity, verbose)
P = conditional_P + conditional_P.T
sum_P = np.maximum(np.sum(P), MACHINE_EPSILON)
P = np.maximum(squareform(P) / sum_P, MACHINE_EPSILON)
return P
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
def _kl_divergence_bh(params, P, neighbors, degrees_of_freedom, n_samples,
n_components, angle=0.5, skip_num_points=0,
verbose=False):
"""t-SNE objective function: KL divergence of p_ijs and q_ijs."""
params = astype(params, np.float32, copy=False)
X_embedded = params.reshape(n_samples, n_components)
neighbors = astype(neighbors, np.int64, copy=False)
if len(P.shape) == 1:
sP = squareform(P).astype(np.float32)
else:
sP = P.astype(np.float32)
grad = np.zeros(X_embedded.shape, dtype=np.float32)
error = _barnes_hut_tsne.gradient(sP, X_embedded, neighbors,
grad, angle, n_components, verbose,
dof=degrees_of_freedom)
c = 2.0 * (degrees_of_freedom + 1.0) / degrees_of_freedom
grad = grad.ravel()
grad *= c
return error, grad

t-SNE – Laurens van der Maaten这个链接是作者收集的各种tsne变种及相关实现。