【数据可视化】SNE

发表于 更新于

遇到了t分布。

这学期正巧选了一门数据可视化的课程，不过似乎老师并没有打算讲数学原理方面的东西，而纯粹是以工程的角度秀自己之前的成果……因此也不期望从课上学到什么。但还是在老师的PPT上看到了一些酷炫的数据可视化的图片。一直好奇一些论文中漂亮的散点图是怎么绘制的，比如像下面这种将Fashion MNSIT种的10种衣物一齐展示在一张图上，并且物以类聚，简洁明了。

之前猜测这种图可能是用到了PCA或者Kmeans这些方式，比如先用某个“黑盒子”提取出几种特征，然后再使用PCA筛选最终展示效果规定的维度（比如上面这个例子就是2维），如果使用神经网络的话，似乎直接用全连接层控制输出的维度就好了。不过关于这种猜测并没有真正实践（狗头）。后来在看CS224N的时候，Christopher Manning老师正好提到了t-SNE，于是就开始了这趟递归学习……

SNE

SNE的全称是stochastic neighbor embedding，而且是Hinton提出的。看到“stochastic”和“embedding”的第一感觉就是这个算法一定是某种迭代优化的算法，而不是一种直接映射的方式。

算法的目的是将高维数据分布的特点迁移到低维数据上，如果两个点在高维空间靠近的话，那么在映射到低维空间时，这两个点的距离也应当是靠近的。

先是如何衡量高维上的靠近。最简单的是算出每一对点的距离，不过Hinton并没有这么做，而是使用另一种衡量方式。假设高维上有两个点$x_i$和点$x_j$，定义点i靠近点j的程度为

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

$\sigma_i$表示以$x_i$为中心的高斯分布的方差。
从这个式子可以看出这个式子定义的靠近是概率意义上的，并且这种所谓的靠近并不对称，因为$p_{j|i} \neq p_{i|j}$，就好像KL散度的非对称性一样。实际上之后的损失函数的确用到了KL散度。

然后，高维点集X要被映射到低维上要进行可视化了，假设映射关系为$x_i \rightarrow y_i$。并且定义新点集距离分布为

$$
q_{j|i} = \frac{exp(-\left|y_i-y_j\right|^2)}{\sum_{k \neq i}{exp(-\left|y_i-y_k\right|^2)}}
$$

因为这里的分布是将要计算的，因此可以自定义分布的方差，简单起见就为$1/\sqrt{2}$好了，这样式子就不出现$\sigma$了。

将视角切换到点$x_i$上，点$x_i$到其他点的距离分布最终转化为了概率分布，不妨假设其为$P_i$；同理，低维空间的点$y_i$也拥有一个分布$Q_i$。如何衡量两个分布的相似程度？KL散度出现了。

定义损失函数为

$$
C = \sum_{i}KL(P_i|Q_i) = \sum_{i}{\sum_{j}{p_{j|i}\log{\frac{p_{j|i}}{q_{j|i}}}}}
$$

t-SNE

t-SNE是对SNE的改进。首先是 __对称问题__。因为KL散度是非对称的，KL(P|Q)一般不等于KL(Q|P)，为了做到对称，简单的方法就是两者兼顾，即

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

n表示点的数量。这是高维空间的概率，低维空间的改进是考虑全局距离

$$
q_{ij} = \frac{exp(-\left|y_i-y_j\right|^2)}{\sum_{k \neq l}{exp(-\left|y_k-y_l\right|^2)}}
$$

和SNE的$q_{j|i}$对比一下，可以发现分母计算了每一对点的距离。

损失函数变为

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

其次是 数据拥挤 的问题。考虑在n维空间上随机撒点，那么这些点距离原点的距离分布的问题。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import norm

plt.figure(figsize=(30, 4))

for i, dim in enumerate((1, 2, 4, 8, 16, 32)):
    # 在dim维空间随机撒点
    pts = [np.random.rand(dim) for j in range(1000)]
    dst = list(map(lambda x: norm(x), pts))
    # 绘图
    ax = plt.subplot(1, 6, i + 1)
    ax.set_xlabel('distance')
    if i == 0:
        ax.set_ylabel('count')
    ax.hist(dst, bins=np.linspace(0, np.sqrt(dim), 32))
    ax.set_title('m={0}'.format(str(dim)), loc='left')

plt.show()

可以看出，随机点离原点的距离变化随着维度升高，是趋向于敏感的，高维空间点的分布更加集中。这会给直接降维带来麻烦，若不加处理直接将高维空间距离分布映射到低维空间，那么点会聚成一团；即使使用正则化“拉伸”这些点的距离，效果也不会很明显。

需要注意的是，点是随机撒的，那么无论是高维空间还是低维空间，空间中点的分布都是大致均匀的，并不存在挤成一团的问题。计算距离也算是一种降维处理，这里就能体现出简单降维处理存在的拥挤问题。

比较理想的处理方式是使用一种映射方式，让高维距离分布拥有低维距离分布“矮胖”的特点，这就用到了t-分布。

t-分布的自由度越低，曲线越矮。如果把高维空间的分布视为正态分布，将要映射的分布视为自由度为1的t-分布，就可以缓解数据拥挤的问题！这样，高维数据距离公式不需要变动，只需要将低维距离公式改变，让其满足t-分布的特点即可。

t-分布这种“矮胖”的身材还能很好的消除 小样本 中存在离群点的缺点。

$q_{ij}$的公式贴在下面。之所以使用自由度为1的t-分布，是因为这种形式符合“平方反比定律”，能够在低维空间中很好地表示巨多的点对之间的关系（作者说的:)）。

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

衡量两个距离分布相似性还是使用KL散度。计算可以直接使用梯度计算

$$
\frac{\partial C}{\partial y_i} = 4 \sum_{j}{(p_{ij} - q_{ij})(y_i - y_j)(1 + \left|y_i - y_j\right|^2)^{-1}}
$$

最后说说t-SNE算法的缺点吧，其实缺点已经很明显了……和训练词向量一样，当数据量很大的时候，会算得很慢，并且每次训练的结果都是不一样的。为了解决速度的问题，有很多人尝试使用“树”结构来加速（比如kd-tree？），还有人使用“负采样”来加速，说到底也是想用小样本估计总体样本，虽然精度可能会出现偏差，不过反正训练出来的结果是给人看的，又何必在意一两个点的偏移呢？

词向量降维可视化

既然本篇是在学习CS224N的过程中深挖的产物，那么实验环节用词向量降维是再合适不过的了。通常使用的词向量动辄上百维，非常不适合碳基生物的理解，如果需要查看训练之后词向量表达意思的能力，至少需要将其降维到3维空间。

以GLOVE为例，glove.6B.50d为高维数据，然后要将其投射到3维平面。该数据集有40w条预训练的词向量，每条词向量拥有50个维度。为了能比较清楚地看出词和词之间的关系，使用2维空间可视化这些向量（3维能显示更多的信息，但是显示效果容易受视角的影响）。

训练代码附在后面，训练结果如下：

代码

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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
import copy
import math

import matplotlib.pyplot as plt
import numpy as np
from tqdm import tqdm

# 读取GLOVE预训练词向量
filename = r"D:\data\trained_net\glove.6B.50d.txt"
glove = {}

with open(filename, encoding="utf8") as f:
    for line in f.readlines():
        line = line.split()
        glove[line[0]] = list(map(lambda w: float(w), line[1:]))

# ！注意words自己从百度上找，我是用的是四级高频词汇，一共700个！
words = [...]

# 筛选存在GLOVE的词
selected_vecs = {word: glove[word] for word in words if word in glove.keys()}
selected_words = list(selected_vecs.keys())
sigma = [np.var(selected_vecs[w]) for w in selected_words]
n = len(selected_vecs)

# 距离矩阵
_p = np.zeros((n, n))


def vec_distance(word1, word2):
    # 计算两个词向量差的距离
    vec1 = np.array(selected_vecs[word1])
    vec2 = np.array(selected_vecs[word2])
    vec = vec1 - vec2
    return np.dot(vec, vec)


# 计算高维向量分布矩阵
for i in tqdm(range(n)):
    dominator = 0.0  # 计算分母
    for k in range(n):
        if i == k:
            continue
        d = vec_distance(selected_words[i], selected_words[k])
        dominator += np.exp(-d / (2 * sigma[i]**2))
    # 计算p矩阵的项
    for j in range(n):
        if i == j:
            continue
        d = vec_distance(selected_words[i], selected_words[j])
        _p[i][j] = np.exp(-d / (2 * sigma[i]**2)) / dominator
    pass

p = (_p.T + _p[j][i]) / (2 * n)


# 训练低维空间距离
EPOCH = 100
DIM = 2
lr = 1e3
pre_kl = []
window = 5

y = np.random.randn(n, DIM)  # 默认是服从N(0,1)分布


def KL(P, Q, epsilon=1e-10):
    _P = P + epsilon
    return (_P * np.log2(_P / (Q + epsilon))).sum()


for epoch in range(EPOCH):
    # 计算分母
    dy = y.reshape(n, 1, DIM) - y.reshape(1, n, DIM)  # n * n * 3
    dominator = (dy**2).sum(axis=2) + 1         # n * n * 1
    dominator = 1 / dominator                   # n * n * 1
    dominator_sum = dominator.sum()
    # 低维距离矩阵
    q = dominator / dominator_sum

    # 计算loss
    kl = KL(p, q)
    print(f"epoch {epoch + 1}/{EPOCH}\tKL = {kl:.6f}\tLR = {lr}")
    # 动态调整学习率
    if len(pre_kl) > 0:
        kl_mean = sum(pre_kl) / len(pre_kl)
        if kl_mean > kl:
            lr *= 1.2
        else:
            lr *= 0.5
    pre_kl.append(kl)
    if len(pre_kl) > window:
        pre_kl.pop(0)

    d_pq = p - q
    # 计算梯度
    grad = np.zeros((n, DIM))
    for i in range(n):
        _grad = 0
        for j in range(n):
            _grad += d_pq[i][j] * dy[i][j] * dominator[i][j]
        grad[i] = _grad * 4
    # 更新低维向量
    y = y - grad * lr

y -= y.mean(axis=0)

# 词向量展示
SHOW = 150  # 为了防止太过密集，选取部分词向量展示

if DIM == 3:
    from mpl_toolkits.mplot3d import Axes3D # 需要有Axes3D才能绘画
    
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    
    ax.scatter(y[:SHOW, 0], y[:SHOW, 1], y[:SHOW, 2], s=1)
    for i in range(SHOW):
        word = selected_words[i]
        ax.text(y[i, 0], y[i, 1], y[i, 2], word, fontsize=10)
else:
    plt.scatter(y[:SHOW, 0], y[:SHOW, 1], s=1)
    for i in range(SHOW):
        word = selected_words[i]
        plt.text(y[i, 0], y[i, 1], word, fontsize=10)

plt.show()

参考