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

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