理解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
28
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变种及相关实现。