数据科学：t-SNE 的原理与应用

Posted on Edited on In Waline:

使用 19 世纪文献生成的词嵌入的 T-SNE 可视化

t 分布随机邻域嵌入(t-SNE) 是一种统计方法，通过为每个数据点在二维或三维地图中赋予一个位置来可视化高维数据。它是一种非线性降维技术，用于在二维或三维低维空间中嵌入高维数据以进行可视化。具体而言，它通过二维或三维点对每个高维对象进行建模，这样相似的对象由附近的点建模，而不同的对象则以高概率由远处的点建模。
t-SNE 算法包括两个主要阶段。首先，t-SNE 在高维对象对上构建一个概率分布，这样相似的对象被分配更高的概率，而不同的点被分配更低的概率。其次，t-SNE 在低维图中的点上定义一个相似的概率分布，并最小化两个分布之间关于图中点位置的Kullback-Leibler 散度（KL 散度）。
对于包含n 个元素的数据集，t-SNE 的运行时间为$O(n^2)$，空间为$O( n^2)$。

实现原理

给定一组高维数据对象$x_i$,…,$x_N$,t-SNE计算$x_i$与$x_j$的相似度成正比的概率

$$
p_{i|j} = \frac{\exp\left(-|x_i - x_j|^2 / 2 \sigma_i^2\right)}{\sum_{k \neq i} \exp\left(-|x_i - x_k|^2 / 2 \sigma_i^2\right)}( i\neq j)
$$

其中$\sigma_i$是以数据点$x_i$为中心的高斯方差,$|x_i - x_j|$表示数据$x_i$与$x_j$的欧式距离,且定义$p_{i|i}$=0;$p_{i|j}$表示$x_i$按照其概率密度的比例并以$x_i$为高斯中心选择$x_j$作为邻居的条件概率;
且有定义$p_{ij} = \frac{p_{j|i} + p_{i|j}}{2N}$。
t-SNE旨在学习一个映射,通过该映射,高维数据$x_i$,…,$x_N$将映射为对应的低维数据$y_i$,…,$y_N$(维数一般为2或者3),该低维数据之间计算所得相似度$q_{ij}$尽可能接近$p_{ij}$。$q_{ij}$定义如下:

$$
q_{ij} = \frac{(1+|y_i - y_j|^2 )^{-1}}{\sum_{k}\sum_{l \neq k}(1+|y_k - y_l|^2 )^{-1}}( i\neq j)
$$

为了得到逼近高维数据之间计算得到相似度的低维数据，我们需要通过梯度下降（当前采用的主流方法为Barnes-Hut算法）的方法最小化P分布与Q分布的KL散度：

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

代码实现

案例选取一个长度为3072随机生成的二维信号，通过t-SNE计算并使用matplotlib库进行绘制,使用python进行t-SNE计算的库主要有sklearn.manifold.TSNE、fitsne、tsnecuda三个。
测试算法时笔者的测试平台的相关配置信息如下：

上述命令得到的信息依次是：CPU逻辑核数与型号信息、CPU物理个数、每个物理CPU的核数、总CPU逻辑个数（总CPU逻辑个数 = 物理CPU个数 x 每个物理CPU的核数 x 超线程数）。
GPU信息如下所示：

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
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
from transformers import set_seed
import os
import torch
import random
import numpy as np
from sklearn.manifold import TSNE
from fitsne import FItSNE
from tsnecuda import TSNE as CudaTSNE
import math
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator
import time
from functools import reduce


def quzheng(num:float):
    if num<0:
        return math.floor(num)
    else:
        return math.ceil(num)

if __name__ == "__main__":
    #sure result can repeat
    seed = 1234
    setup_seed(seed)

    #gen random 2d data
    m = 2
    n = 3072 
    data1 = np.random.uniform(-100, 100, size=(n, 1)).astype(np.float32) #random data range (-100,100),type is float32
    data2 = np.random.uniform(-100, 100, size=(n, 1)).astype(np.float32) #random data range (-100,100),type is float32
  
    data = np.column_stack((data1, data2))
    origin_data = data.copy()
  
    data_dict = dict()

    tsne = TSNE(n_components = 2,perplexity=30,verbose=1,max_iter=2000,random_state=seed,init='random')
    start_time = time.perf_counter()
    tsne_data = tsne.fit_transform(data)
    end_time = time.perf_counter()
    execution_time = end_time-start_time
    print("tsne excute cost ",execution_time)
    if not np.array_equal(data,origin_data):
        print("data change!")
        exit(0)
    data_dict['tsne'] = tsne_data


    start_time = time.perf_counter()
    fit_tsne_data = FItSNE(data.astype(np.float64), perplexity=30, max_iter=2000,fft_not_bh=True,rand_seed=seed,initialization='random')
    end_time = time.perf_counter()
    execution_time = end_time-start_time
    print("fit-tsne excute cost ",execution_time)
    if not np.array_equal(data,origin_data):
        print("data change!")
        exit(0)
    data_dict['fit-tsne'] = fit_tsne_data

    cuda_tsne = CudaTSNE(n_iter=2000, verbose=1, perplexity=30,random_seed=seed)
    start_time = time.perf_counter()
    cuda_tsne_data = cuda_tsne.fit_transform(data)
    end_time = time.perf_counter()
    execution_time = end_time-start_time
    print("cuda-tsne excute cost ",execution_time)
    if not np.array_equal(data,origin_data):
        print("data change!")
        exit(0)
    data_dict['cuda-tsne'] = cuda_tsne_data
  

  
    xy_min_max_list =  list(map(quzheng,reduce(lambda x,y:[min(x[0],y[0]),min(x[1],y[1]),max(x[2],y[2]),max(x[3],y[3])],map(lambda x:[np.min(x[:,0]),np.min(x[:,1]),np.max(x[:,0]),np.max(x[:,1])],data_dict.values()))))
    print("get xy_range is ",xy_min_max_list)
    x_min = xy_min_max_list[0]
    y_min = xy_min_max_list[1]
    x_max = xy_min_max_list[2]
    y_max = xy_min_max_list[3]
  
    plt.figure(figsize=(32, 24))
    plt.suptitle(f"Tsne Test",fontsize=40,    # 字体大小
                fontweight='bold',  # 字体粗细
                fontstyle='italic',  # 字体风格（斜体）
                color='blue'   # 字体颜色
            )
    keys_list = list(data_dict.keys())
    for i in  range(0,len(keys_list)):
        sample = data_dict[keys_list[i]]
        ax = plt.subplot(1,len(keys_list),i+1)
        ax.xaxis.set_major_locator(MultipleLocator(25))  # x轴主刻度间隔为25
        ax.yaxis.set_major_locator(MultipleLocator(25))  # y轴主刻度间隔为25
        ax.scatter(sample[:, 0], sample[:, 1], s=5, alpha=0.6)
        ax.set_title(f'{keys_list[i]}')
        ax.set_ylim(y_min,y_max)
        ax.set_xlim(x_min,x_max)
        ax.set_aspect('equal', 'box')
    plt.tight_layout()
    save_file_path = "./tsne"
    plt.savefig(save_file_path,dpi=300)
    plt.close()

运行得到上述三种算法时间的花销如下：

由于tsnecuda底层实现并未采用随机数种子，而是基于当前时间生成的种子,详情参见issue所以对于同一数据重复运行也无法得到相同结果
上述运行结果可视化如下: