遇到了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维空间上随机撒点，那么这些点距离原点的距离分布的问题。

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维能显示更多的信息，但是显示效果容易受视角的影响）。

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

代码

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()

【数据可视化】SNE

SNE

t-SNE

词向量降维可视化

代码

参考