手撕t-SNE算法

t-SNE（t-Distributed Stochastic Neighbor Embedding）即t分布随机邻近嵌入，读作Tee-Snee，Laurens van der Maaten等人在2008年发表于JMLR上的论文《Visualizing Data using t-SNE》。类似于PCA，是一种数据降维方法，与PCA不同的是，它是非线性的。

1、SNE

讲t-SNE之前，得先讲SNE，Geoffrey Hinton等人2002年发表在NIPS上的《Stochastic Neighbor Embedding》，其实t-SNE和SNE论文的共同作者都是Hinton大佬。

假设我们的原始数据矩阵为 X ∈ R n × m \mathbf{X} \in \mathbb{R}^{n \times m} X∈Rn×m，降维后的矩阵为 Y ∈ R n × 2 \mathbf{Y} \in \mathbb{R}^{n \times 2} Y∈Rn×2，也就是说从 m m m维降为2维。

对于 X \mathbf{X} X：

数据点之间的条件概率为：

p j ∣ i = e − ∥ x i − x j ∥ 2 2 σ i 2 ∑ k ≠ i e − ∥ x i − x k ∥ 2 2 σ i 2 p_{j \mid i}=\frac{e^{- \frac{\left \| x_i - x_j \right \|^2}{2 \sigma_i^2}}}{\sum_{k \ne i} e^{- \frac{\left \| x_i - x_k \right \|^2}{2 \sigma_i^2}}} pj∣i​=∑k​=i​e−2σi2​∥xi​−xk​∥2​e−2σi2​∥xi​−xj​∥2​​

∥ x i − x j ∥ 2 \left \| x_i - x_j \right \|^2 ∥xi​−xj​∥2是也就是欧几里得距离的平方，然后再套一个Softmax使它变成有效的概率。

n n n个数据点，两两分别按上式计算，得到概率矩阵 P n × n P_{n \times n} Pn×n​， i i i就是第几行， j j j就是第几列：

( p 1 ∣ 1 p 2 ∣ 1 ⋯ p n ∣ 1 p 1 ∣ 2 p 2 ∣ 2 ⋯ p n ∣ 2 ⋮ ⋮ ⋮ ⋮ p 1 ∣ n p 2 ∣ n ⋯ p n ∣ n ) \begin{pmatrix} p_{1|1} & p_{2|1} & \cdots & p_{n|1} \\ p_{1|2} & p_{2|2} & \cdots & p_{n|2} \\ \vdots & \vdots & \vdots & \vdots\\ p_{1|n} & p_{2|n} & \cdots & p_{n|n} \end{pmatrix} ⎝⎜⎜⎜⎛​p1∣1​p1∣2​⋮p1∣n​​p2∣1​p2∣2​⋮p2∣n​​⋯⋯⋮⋯​pn∣1​pn∣2​⋮pn∣n​​⎠⎟⎟⎟⎞​

其中 p i ∣ i = 0 p_{i \mid i}=0 pi∣i​=0，也就是说 P P P是一个对角线都是0的矩阵。

那 σ ∈ R n \mathbf{\sigma} \in \mathbb{R}^n σ∈Rn是啥？

用过t-SNE进行可视化就知道有个参数叫困惑度（Perplexity），困惑度 K K K：

K = 2 H ( P i ) K=2^{H(P_i)} K=2H(Pi​)

其中 H H H是香农熵，即：

H ( P i ) = − ∑ j p j ∣ i l o g ( p j ∣ i ) H(P_i)=-\sum_j p_{j|i}log(p_{j|i}) H(Pi​)=−j∑​pj∣i​log(pj∣i​)

困惑度 K K K和 σ i \sigma_i σi​的关系是，当我们设定好 K K K之后， σ i \sigma_i σi​的值也为之确定，二者是单调递增关系，所以可以通过二分搜索来确定 σ i \sigma_i σi​。

例如我们设置K=32，一开始假设范围是 0.001 < σ i < 1000 0.001<\sigma_i<1000 0.001<σi​<1000，猜测 σ i = 1000 + 0.001 2 = 500.0005 \sigma_i=\frac{1000+0.001}{2}=500.0005 σi​=21000+0.001​=500.0005，然后求 H ( P i ) H(P_i) H(Pi​)，再求得 K = 10 K=10 K=10，K小了，所以 σ i \sigma_i σi​在 500.0005 < σ i < 1000 500.0005<\sigma_i<1000 500.0005<σi​<1000范围内，再猜测 σ i = 1000 + 500.0005 2 = 750.00025 \sigma_i=\frac{1000+500.0005}{2}=750.00025 σi​=21000+500.0005​=750.00025，以此类推，直到迭代次数满足或者算得的K比真实的K=32相差不多时停止。

注意 σ \mathbf{\sigma} σ是一个维度为 n n n的向量，也就是说有多少个点，就有多少个 σ i \sigma_i σi​。困惑度K只有一个值。

σ \mathbf{\sigma} σ越大， P P P中每一行概率单纯形就越平均（趋向 1 n \frac{1}{n} n1​），熵就越大，困惑度 K K K也就越大。

对于 Y \mathbf{Y} Y：

数据点之间的条件概率为：

q j ∣ i = e − ∥ x i − x j ∥ 2 ∑ k ≠ i e − ∥ x i − x k ∥ 2 q_{j \mid i}=\frac{e^{- \left \| x_i - x_j \right \|^2}}{\sum_{k \ne i} e^{- \left \| x_i - x_k \right \|^2} } qj∣i​=∑k​=i​e−∥xi​−xk​∥2e−∥xi​−xj​∥2​

也就是没有了 σ i \sigma_i σi​。

同理也可以得到概率矩阵 Q n × n Q_{n \times n} Qn×n​。

损失函数

要想使得两个分布差异尽可能小，可以用KL散度度量两个概率分布的差异，最小化差异从而得到损失函数：

C = min ⁡ Y D K L ( P ∥ Q ) = min ⁡ Y ∑ i P i l o g P i Q i \begin{aligned} C & =\underset{\mathbf{Y}}{\min} D_{KL}(P \| Q) \\ & =\underset{\mathbf{Y}}{\min} \sum_i P_i log \frac{P_i}{Q_i}\\ \end{aligned} C​=Ymin​DKL​(P∥Q)=Ymin​i∑​Pi​logQi​Pi​​​

使用梯度下降优化就可以得到 Y \mathbf{Y} Y啦。

2、t-SNE

t-SNE相比SNE有两个区别：

1、对称SNE；

2、对Q的计算进行了修改，使其服从t分布，从而缓解“拥堵”问题。

参考文献

[1] In Raw Numpy: t-SNE