【机器学习】Tensorflow：tSNE数据非线性降维

深度学习巨头之一的Hinton大神在数据降维领域有一篇经典论文Visualizing Data using t-SNE。该方法是流形（非线性）数据降维的经典，从发表至今鲜有新的降维方法能全面超越。该方法相比PCA等线性方法能有效将数据投影到低维空间并保持严格的分割界面；缺点是计算复杂度大，一般推荐先线性降维然后再用tSNE降维。python sklearn有相应的实现。我现在用Tensorflow实现这个算法。
tSNE的核心思想是：把高维数据看做高维空间中的点$x_i$，然后用流形方法将其映射到低维空间中的$y_i$，要求保持其空间距离。即高维空间相距较近/远的点，映射到低维空间中仍然较近/远。为了让挨得近的点靠的更近，距离用高斯函数来度量：

$$D_{ij} = e^{-\frac{|x_i-x_j|^2}{2\sigma ^2}}/ \sum_{i,j}{e^{-\frac{|x_i-x_j|^2}{2\sigma ^2}}}$$ 为了均衡点之间的斥力，我们要让 $$p_{ij}=p_{ji}=\frac{D_{ij}+D_{ji}}{2}$$
为了避免噪声、离群点对低维空间映射的干扰，Hinton在低维空间中类比统计学中t分布的公式定义了一个新的距离: $$q_{ij}=\frac{(1+|y_i-y_j|^2)^{-1}}{\sum_{i,j}{(1+|y_i-y_j|^2)^{-1}}}$$
为了保持映射中，远的还是远的，近的还是近的，我们要用KL散度来定义损失，KL散度可以反映两个分布的正关系，如果为0说明两个分布一致： $$loss = \sum{p_{ij} \log {\frac{p_{ij}}{q_{ij}}}}$$
最优化这个损失函数就得到低维空间中的映射$y_i$。这就是tSNE。求解这个损失函数，可以用梯度下降法，对损失函数求导，原文用了相当长的篇幅来解出求导公式： $$\frac{\partial loss}{\partial y_i}=4\sum_j{(p_{ij}-q_{ij})(1+|y_i-y_j|^2)^{-1}(y_i-y_j)}$$
现在我不管这个求导公式，直接上Tensorflow。我用sklearn中的iris数据集来做测试，这是一个[150,4]的数据，其中0-49为1类，50-99为2类，100-149为3类。为了方便计算，我对该数据进行了归一化。$p_{ij}$可以预先用numpy求出来，因为它与优化过程无关。
一般来说直接用tensorflow构建loss函数即可，然而tensorflow求解两两之间的距离是个问题。tensorflow不支持单独操作tensor的某个元素（也许未来会支持?），因此没有办法提取每个样本来循环求解所有的距离。这里就需要一个矩阵技巧。
我知道$|x_i-x_j|^2=x_i^2-2x_ix_j+x_j^2$，我要构建的距离矩阵是[150,150]的对称矩阵pairD，样本矩阵是[150,4]的矩阵X，由上述公式知道 $$pairD = sum(X*X,axis=1)-2X^TX+sum(X^T*X^T,axis=1)$$ 其中*表示元素乘法，做加法的时候默认使用broadcast广播。
另外我发现在构建loss函数时，有$p_{ij,i=j}=0$，为了方便计算没有排除这个元素，导致出现log函数自变量为0的情况，于是无法求解。为了解决这个问题，我强制$p_{ij}=max(p_{ij},0.000001)$，当然也可以设置其他比较小的值。这也是频域逆滤波的常用方法。
然后就是构建tensorflow的损失函数了：

with tf.device(‘/cpu:0‘):
    X = tf.placeholder(‘float‘,(150,150))
    initial = tf.random_normal([150,2]) * 0.0001#映射到二维空间
    Y = tf.Variable(initial)
    A = tf.reduce_sum(Y*Y, axis=1)
    A = tf.reshape(r, [-1, 1])
    #pair wise distance
    pairD = A - 2*tf.matmul(Y, tf.transpose(Y)) + tf.transpose(A) + 1.
    qij = 1./pairD
    sumq = tf.reduce_sum(qij,axis=1)
    qij /= sumq
    loss = tf.reduce_sum( X*tf.log(X / qij) )
    global_step = tf.Variable(0, name = ‘global_step‘,trainable=False)
    starter_learning_rate = 0.1
    learning_rate = tf.train.exponential_decay(starter_learning_rate, global_step,20, 0.95, staircase=True)
    train_op = tf.train.AdamOptimizer(learning_rate=learning_rate).minimize(loss=loss,global_step = global_step) 

直接使用原始数据绘图第0、1维：

求解之后，绘制图像得到

经过tSNE降维映射，不同的数据在低维空间中依然能够被分开。整个最优化过程空间复杂度较高；此外不同的参数会导致最终的降维图不一样，但是降维图中的3类数据全部都被分隔开了。这个降维方法用于可视化居多。全程由Tensorflow自动最优化得到，如果使用GPU的话更加快，tSNE计算慢的问题可以解决了，空间复杂度有点头疼。