机器学习报告——非线性降维算法t-SNE

1. 摘要

为了了解在低维、非线性流型上表示高维数据的算法，我选择了t-SNE算法进行学习。本次学习中，我阅读了论文Visualizing Data using t-SNE(2008)，了解算法的详细内容，阅读并分析了sklearn源码中的TSNE类，之后自己进行了手写数字可视化的算法实践。在本篇报告中，我将详细介绍我理解的t-SNE算法和它要求的预备知识，并展示自己的算法实践，最后对本次学习进行总结并为自己列出了几个可继续学习的方向。

2. 选题背景

高维数据的可视化在许多领域都有着特别重要的地位。在机器学习的课程中我们已经学习了一种传统的降维算法：主成分分析（PCA）。通过PCA，我们可以很好地将一些多维数据转化到低维空间上。但PCA是一种线性降维算法，它不能解释特征之间的复杂多项式关系。同时，线性降维算法更关注于将不相似的数据点以远距离的形式在低维度中表示出来。换句话来说，PCA更关注全局结构。但对于在低维、非线性流型上表示高维数据，通常更为重要的是将非常相似的数据点的低维表示保持在一起，也就是局部结构。但这并不是线性降维算法所能做的，因此我进一步选择了能够很好地捕获高维数据的大部分局部结构,且能揭示高维数据全局结构的非线性降维算法t-SNE进行学习。

3. 算法介绍

3.1 Stochastic Neighbor Embedding (SNE)

t-SNE是在SNE的基础上进行的改进。SNE算法是基于“高维空间中相似的数据点，映射到低维空间中的距离也是相似的。”的想法而提出的。SNE通过将数据点之间的欧式距离转化为条件概率来表示数据点之间的相似性。

3.1.1 高维空间

数据点$x_{i}$xi​和$x_{j}$xj​之间的相似度用条件概率$p_{j|i}$pj∣i​来表示，即：$x_{i}$xi​以条件概率$p_{j|i}$pj∣i​选择$x_{j}$xj​作为它的邻近点。考虑以$x_{i}$xi​为中心点的高斯分布，若$x_{j}$xj​越靠近$x_{i}$xi​，则$p_{j|i}$pj∣i​越大。反之，若两者相距较远，则$p_{j|i}$pj∣i​极小。因此，$p_{j|i}$pj∣i​的定义如下：$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})}$pj∣i​=∑k̸​=i​exp(−∥xi​−xk​∥2/2σi2​)exp(−∥xi​−xj​∥2/2σi2​)​其中$\sigma_{i}$σi​表示以$x_{i}$xi​为中心点的高斯分布的方差，$p_{i|i}=0$pi∣i​=0。
关于高斯分布的方差$\sigma_{i}$σi​的选择
由于数据密度在区域之间会有不同，所以通常没有一个单一的$\sigma_{i}$σi​能对数据集中的所有数据点达到最优状态。每一个$\sigma_{i}$σi​会在所有数据点上产生一个概率分布$P_{i}$Pi​，该分布的熵随$\sigma_{i}$σi​的增加而增加。SNE对$\sigma_{i}$σi​的值进行二分检索，由该$\sigma_{i}$σi​产生一个具有用户指定困惑度的分布$P_{i}$Pi​。论文中将困惑度定义为：$Prep(P_{i})=2^{H(P_{i})}$Prep(Pi​)=2H(Pi​)$H(P_{i})=-\sum _{j}p_{j|i}log_{2}p_{j|i}$H(Pi​)=−j∑​pj∣i​log2​pj∣i​

3.1.2 低维空间

低维空间中的数据表示同样用条件概率的形式描述：假设高维数据点$x_{i}$xi​和$x_{j}$xj​在低维空间的映射点分别为$y_{i}$yi​和$y_{j}$yj​，则低维空间中的条件概率用$q_{j|i}$qj∣i​表示，并将所有高斯分布的方差均设定为$\frac{1}{\sqrt{2}}$2​1​，则有：$q_{j|i}= \frac{exp(-\left \| y_{i}- y_{j}\right \|^{2})}{\sum _{k\neq i}exp(-\left \| y_{i}- y_{k}\right \|^{2})}$qj∣i​=∑k̸​=i​exp(−∥yi​−yk​∥2)exp(−∥yi​−yj​∥2)​同理，$q_{i|i}=0$qi∣i​=0。

3.1.3 cost function

条件概率分布$P_{i}$Pi​表示：$x_{i}$xi​与其他所有点之间的条件概率。同理在低维空间存在一个条件概率分布$Q_{i}$Qi​且应该与$P_{i}$Pi​一致，用KL距离表示两个分布间的相似性。SNE的目标就是对所有数据点最小化这个KL距离，我们可以使用梯度下降算法最小化如下代价函数：$C=\sum_{i}^{ }KL\left ( P_{i}|| Q_{i}\right )=\sum_{i}^{ }\sum_{j}^{ }p_{j|i}log\frac{p_{j|i}}{q_{j|i}}$C=i∑​KL(Pi​∣∣Qi​)=i∑​j∑​pj∣i​logqj∣i​pj∣i​​SNE代价函数对$y_{i}$yi​求梯度后的形式如下：$\frac{\delta C}{\delta y_{i}}=2\sum_{j}^{ }\left ( p_{j|i}- q_{j|i}+p_{i|j}-q_{i|j}\right )\left ( y_{i}- y_{j} \right )$δyi​δC​=2j∑​(pj∣i​−qj∣i​+pi∣j​−qi∣j​)(yi​−yj​)梯度更新$\gamma ^{\left ( t \right )}= \gamma ^{\left ( t-1 \right )}+\eta \frac{\delta C}{\delta \gamma }+\alpha \left ( t \right )\left (\gamma ^{\left ( t-1 \right )}-\gamma ^{\left ( t-2 \right )} \right )$γ(t)=γ(t−1)+ηδγδC​+α(t)(γ(t−1)−γ(t−2))

3.1.4 SNE的限制

• 距离不对称：$p_{i|j}$pi∣j​与$p_{j|i}$pj∣i​是不相等的，低维空间中$q_{i|j}$qi∣j​与$q_{j|i}$qj∣i​也是不相等的
• SNE更关注局部结构，忽视了全局结构：理论上说，高维空间中两个数据点距离较近但映射到低维空间后距离较远，和高维空间中两个数据点距离较远但映射到低维空间距离较近，都应该得到一个很高的惩罚值。但我们分析SNE的代价函数可以发现：当$p_{j|i}$pj∣i​较大，$q_{j|i}$qj∣i​较小时，代价$C$C较高；而$p_{j|i}$pj∣i​较小，$q_{j|i}$qj∣i​较大时，代价$C$C较低，也就是说SNE对于将距离远的点映射成距离近时的惩罚不够。

3.2 t-Distributed Stochastic Neighbor Embedding (t-SNE)

3.2.1 对称SNE

为了能在高维和低维空间中能分别构造更加通用、合理的联合概率分布P和Q，可以得出要求：对任意$i$i，$j$j，均有$p_{i|j}=p_{j|i}$pi∣j​=pj∣i​,$q_{i|j}=q_{j|i}$qi∣j​=qj∣i​。
所以在低维空间中，$q_{ij}$qij​定义如下：$q_{ij}= \frac{exp(-\left \| y_{i}- y_{j}\right \|^{2})}{\sum _{k\neq l}exp(-\left \| y_{k}- y_{l}\right \|^{2})}$qij​=∑k̸​=l​exp(−∥yk​−yl​∥2)exp(−∥yi​−yj​∥2)​
而在高维空间中，我们这样定义$p_{ij}$pij​：$p_{ij}= \frac{p_{j|i}+p_{i|j}}{2n}$pij​=2npj∣i​+pi∣j​​其中$n$n为数据点的总数。之所以不像低维空间中那样定义，是为了保证该距离既能满足对称性，又能保证某一离群点$x_{i}$xi​的惩罚值不会过小。若有一个离群点$x_{i}$xi​，它与所有其他数据点之间的距离都较大，那么对所有$j$j，$p_{i|j}$pi∣j​的值均较小。则无论该离群点在低维空间中的映射点$y_{i}$yi​处在什么位置，惩罚值都不会太高，这是不符合理论分析的。

KL代价函数为：$C=\sum_{i}^{ }KL\left ( P|| Q\right )=\sum_{i}^{ }\sum_{j}^{ }p_{ij}log\frac{p_{ij}}{q_{ij}}$C=i∑​KL(P∣∣Q)=i∑​j∑​pij​logqij​pij​​梯度变为：$\frac{\delta C}{\delta y_{i}}=4\sum_{j}^{ }\left ( p_{ij}-q_{ij}\right )\left ( y_{i}- y_{j} \right )$δyi​δC​=4j∑​(pij​−qij​)(yi​−yj​)
这样的改动提高了梯度计算效率，解决了算法中距离不对称的问题。但相对于原SNE的优化效果比较微弱，因为它没有解决“拥挤问题”。

3.2.2 拥挤问题

对“拥挤问题”最直接的理解就是看图。在SNE的可视化结果图中，不同类别的簇挤在一起，不能很好地区分开来，这就是拥挤问题。由于高维空间距离分布和低维空间距离分布的差异，SNE和对称SNE都存在拥挤问题。
高维空间距离分布和低维空间距离分布的差异表现在：在映射后的二维空间中，能容纳中远距离点的面积没有能容纳邻近点的面积大。也就是说：我们为了能精确地对邻近点进行建模，我们只能将在高维中中远距离的点在二维空间中放得更远。这样一来各个簇之间便不能很好的区分开。
我们希望：同一簇内的点(距离较近)聚合的更紧密，不同簇之间的点(距离较远)更加疏远。由于t-分布相对于高斯分布来说，对异常点更加不敏感，能够活得更好的拟合效果，所以t-SNE中采用t-分布来实现这一目标，解决拥挤问题。

上图中有高斯分布和t分布两条曲线，高斯分布对应高维空间，t分布对应低维空间。**那么为了满足$p_{ij}=q_{ij}$pij​=qij​，对于高维空间中相距较近的点，它们在低维空间中的距离需要更小一点；而对于高维空间中相距较远的点，他们在低维空间中的距离需要更远。**由图可知，t-分布能满足这一要求，所以我们使用自由度为1的t分布重新定义q_{ij}：$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}}$qij​=∑k̸​=l​(1+∥yk​−yl​∥2)−1(1+∥yi​−yj​∥2)−1​此时在KL距离下，梯度变为：$\frac{\delta C}{\delta y_{i}}=4\sum_{j}^{ }\left ( p_{ij}-q_{ij}\right )\left ( y_{i}- y_{j} \right )(1+\left \| y_{i}- y_{j}\right \|^{2})^{-1}$δyi​δC​=4j∑​(pij​−qij​)(yi​−yj​)(1+∥yi​−yj​∥2)−1
对于拥挤问题，还可以用梯度的物理意义解释：分子之间的引力和斥力。低维空间中点$y_{i}$yi​的位置是由其他所有点对其作用力的合力所决定的。其中某个点$y_{j}$yj​对其作用力是沿着$y_{i}-y_{j}$yi​−yj​方向的，具体是引力还是斥力占主导就取决于$y_{i}$yi​与$y_{j}$yj​之间的距离了，其实就与$(p_{j|i}-q_{j|i}+p_{i|j}-q_{i|j})$(pj∣i​−qj∣i​+pi∣j​−qi∣j​)这一项有关。所以那些中远距离的点对数据点$i$i是有轻微的力的作用的，而这种力就会阻碍不同簇之间的间隙形成。

3.2.3 t-SNE小结

t-SNE算法如下

t-SNE在原版SNE上解决了距离不对称和拥挤问题两点限制，同时t-SNE存在以下限制：

• 计算成本高：通过之前的分析我们可以知道t-SNE涉及到大量的运算，所以在面对一些很大的数据时，t-SNE的效率不算很突出。
• 算法是随机的：不同的随机种子会产生不同的嵌入空间效果，但我们自然是希望能产生最优的结果。
• 全局结构没有特别明确的保留：t-SNE能够保存一定的全局结构，但我们含可以通过使用PCA来初始化（init = 'PCA'）点来改善这种情况。

4. 算法实践

4.1 sklearn源码分析

sklearn中用TSNE类实现了t-SNE算法，完成对高维数据的降维操作。

参数说明

参数	描述
n_components	int (default：2),嵌入空间的维度
perplexity	float (default:30),数据集越大所需参数越大,建议值：5~50
early_exaggeration	float (default: 12),控制原始空间中的自然集群在嵌入空间中的紧密程度和其间的空间
learning_rate :	float (default:200),建议取值为10.0-1000.0
n_iter	int (default:1000), 最大迭代次数
n_iter_without_progress	int (default:300), 另一种方法的最大迭代次数,要求为50的倍数
min_grad_norm	float (default:1e-7), 梯度低于该值时停止算法
metric	string or callable, 精确度的计量方式
init	string or numpy array (default:”random”)
verbose	int (default:0), 训练过程是否可视
random_state	int,RandomState instance or None (default:None)，控制随机数的生成
method	string (default:’barnes_hut’), 对于大数据集用默认值，对于小数据集用’exact’
angle	float (default:0.5), 只有method='barnes_hut’时可用

属性说明

属性	描述
embedding_	嵌入向量
kl_divergence_	最后的KL散度
n_iter_	迭代次数

方法说明

方法	描述
fit	将输入数据投影到嵌入空间
fit_transform	将输入数据投影到嵌入空间，并返回转换后的数据

4.2 手写数字可视化

4.2.1 自己的t-SNE

@PengTingyu_201692392

from time import time
import numpy as np
from sklearn import datasets
import matplotlib.pyplot as plt

def neg_squared_euc_dists(X):
    '''
    # 参数:
        X: [N,D]的矩阵
    # 返回值:
        -D：[N,N]的矩阵, -D[i,j]表示X_i & X_j 之间的欧式距离的平方的负数
    '''
    sum_X = np.sum(np.square(X), 1)
    D = np.add(np.add(-2 * np.dot(X, X.T), sum_X).T, sum_X)
    return -D


def softmax(X, diag_zero=True):
    '''
    softmax 计算 exp()/∑exp()
    '''
    # 归一化
    e_x = np.exp(X - np.max(X, axis=1).reshape([-1, 1]))
    # 设置对角线值为0
    if diag_zero:
        np.fill_diagonal(e_x, 0.)
    e_x = e_x + 1e-8  # 数值稳定
    return e_x / e_x.sum(axis=1).reshape([-1, 1])


def calc_prob_matrix(distances, sigmas=None):
    '''
    将距离矩阵转化为概率矩阵
    # 参数：
        distances：距离矩阵
        sigmas:高斯分布的方差
    '''
    if sigmas is not None:
        two_sig_sq = 2. * np.square(sigmas.reshape((-1, 1)))
        return softmax(distances / two_sig_sq)
    else:
        return softmax(distances)


def p_conditional_to_joint(P):
    '''
    将高维中的距离变为对称的
    # 参数：
        P: 条件概率矩阵
    '''
    return (P + P.T) / (2. * P.shape[0])


def p_joint(X, target_perplexity):
    '''
    根据给的原数据X和困惑度，计算最后高维中的联合概率分布矩阵    
    # 参数        
        X: 原数据矩阵    
    # 返回值:        
        P: 高维中的联合概率分布矩阵.
    '''
    distances = neg_squared_euc_dists(X)
    # 找到最优的高斯分布方差sigma
    sigmas = find_optimal_sigmas(distances, target_perplexity)
    p_conditional = calc_prob_matrix(distances, sigmas)
    P = p_conditional_to_joint(p_conditional)
    return P


def binary_search(eval_fn, target, tol=1e-10, max_iter=10000, lower=1e-20, upper=1000.):
    '''
    二分查找
    # 参数：
        eval_fn: 优化的函数
        target: 目标输出值
        tol: Float, guess < tol时停止搜索
        max_iter: int, 最大迭代次数
        lower: Float, 搜索下界
        upper: Float, 搜索上界
    # Returns:
        Float, 优化函数的最优输入值
    '''
    for i in range(max_iter):
        guess = (lower + upper) / 2.
        val = eval_fn(guess)
        if val > target:
            upper = guess
        else:
            lower = guess
        if np.abs(val - target) <= tol:
            break
    return guess


def calc_perplexity(prob_matrix):
    '''
    计算困惑度
    # 参数：
        prob_matrix: 困惑度分布
    '''
    entropy = -np.sum(prob_matrix * np.log2(prob_matrix), 1)
    perplexity = 2 ** entropy
    return perplexity


def perplexity(distances, sigmas):
    '''
    根据距离矩阵和高斯方差计算困惑度
    '''
    return calc_perplexity(calc_prob_matrix(distances, sigmas))


def find_optimal_sigmas(distances, target_perplexity):
    '''
    给距离矩阵的每一行找一个能产生最优困惑度分布的sigma[].
    '''
    sigmas = [] 
    # 距离矩阵的每一行代表一个数据点
    for i in range(distances.shape[0]):
        # Make fn that returns perplexity of this row given sigma
        eval_fn = lambda sigma:perplexity(distances[i:i+1, :], np.array(sigma))
        # 执行二分查找
        correct_sigma = binary_search(eval_fn, target_perplexity)
        sigmas.append(correct_sigma)
    return np.array(sigmas)


def estimate_tsne(X, y, P, rng, num_iters, q_fn, grad_fn, learning_rate, momentum):
    '''
    模拟函数
    # 参数：
        X: 原数据矩阵
        y: labels
        P: 高维中的联合分布矩阵
        rng: np.random.RandomState().
        num_iters: 迭代次数
        q_fn: 获得Y和低维空间中的联合概率分布矩阵Q
    # 返回值:
        Y: X在低维空间中的分布.
    '''
    # Y随机初始化
    Y = rng.normal(0., 0.0001, [X.shape[0], 2])

    if momentum:
        Y_m2 = Y.copy()
        Y_m1 = Y.copy()

    # 梯度下降
    for i in range(num_iters):
        Q, distances = q_fn(Y)
        grads = grad_fn(P, Q, Y, distances)
        # 更新Y
        Y = Y - learning_rate * grads
        if momentum:  
            Y += momentum * (Y_m1 - Y_m2)
            # 更新
            Y_m2 = Y_m1.copy()
            Y_m1 = Y.copy()

    return Y


def q_tsne(Y):
    '''
    根据t-SNE低维空间的分布Y，计算低维空间的联合概率分布矩阵
    '''
    distances = neg_squared_euc_dists(Y)
    inv_distances = np.power(1. - distances, -1)
    np.fill_diagonal(inv_distances, 0.)
    return inv_distances / np.sum(inv_distances), inv_distances


def tsne_grad(P, Q, Y, inv_distances):
    '''
    t-SNE的梯度
    '''
    pq_diff = P - Q
    pq_expanded = np.expand_dims(pq_diff, 2)
    y_diffs = np.expand_dims(Y, 1) - np.expand_dims(Y, 0)

    distances_expanded = np.expand_dims(inv_distances, 2)
    y_diffs_wt = y_diffs * distances_expanded
    grad = 4. * (pq_expanded * y_diffs_wt).sum(1)
    return grad


def plot_embedding(data, label, title):
    '''
    可视化嵌入空间
    '''
    x_min, x_max = np.min(data, 0), np.max(data, 0)
    data = (data - x_min) / (x_max - x_min)

    fig = plt.figure()
    ax = plt.subplot(111)
    for i in range(data.shape[0]):
        plt.text(data[i, 0], data[i, 1], str(label[i]),
                 color=plt.cm.Set1(label[i] / 10.),
                 fontdict={'weight': 'bold', 'size': 9})
    plt.xticks([])
    plt.yticks([])
    plt.title(title)
    return fig


# 参数
PERPLEXITY = 20 
SEED = 1  # 随机种子
MOMENTUM = 0.9 
LEARNING_RATE = 10. 
NUM_ITERS = 500  

def main():
    rng = np.random.RandomState(SEED)
    # 加载数据  
    digits = datasets.load_digits(n_class=6) # 只选取数字0~5进行可视化
    X, y = digits.data, digits.target
    P = p_joint(X, PERPLEXITY)

    # t-SNE
    t0 = time()
    Y = estimate_tsne(X, y, P, rng, num_iters=NUM_ITERS, q_fn=q_tsne, grad_fn=tsne_grad,\
            learning_rate=LEARNING_RATE, momentum=MOMENTUM)
    # 可视化
    fig = plot_embedding(Y, y,
                         't-SNE embedding of the digits (time %.2fs)'
                         % (time() - t0))
    plt.show(fig)

if __name__ == '__main__':
    main()

实验结果

4.2.2 sklearn的t-SNE

@PengTingyu_201692392

# coding='utf-8'
"""t-SNE对手写数字进行可视化"""
from time import time
import numpy as np
import matplotlib.pyplot as plt
from sklearn import manifold, datasets

'''数据准备'''
digits = datasets.load_digits(n_class=6) # 只选取数字0~5进行可视化
X, y = digits.data, digits.target
n_samples, n_features = X.shape

'''显示原始数据图片'''
n = 30  # 每行30个数字，每列30个数字
img = np.zeros((10 * n, 10 * n))
for i in range(n):
    ix = 10 * i + 1
    for j in range(n):
        iy = 10 * j + 1
        img[ix:ix + 8, iy:iy + 8] = X[i * n + j].reshape((8, 8))
plt.figure(figsize=(8, 8))
plt.imshow(img, cmap=plt.cm.binary)
plt.xticks([])
plt.yticks([])
plt.show()

'''t-SNE'''
print('Computing t-SNE embedding')
tsne = manifold.TSNE(n_components=2, init='pca', verbose=1, random_state=501)
t0 = time()
X_tsne = tsne.fit_transform(X)
print("Org data dimension is {}. Embedded data dimension is {}".format(X.shape[-1], X_tsne.shape[-1]))

'''可视化'''
x_min, x_max = X_tsne.min(0), X_tsne.max(0)
X_norm = (X_tsne - x_min) / (x_max - x_min)  # 坐标归一化
plt.figure(figsize=(8, 8))
for i in range(X_norm.shape[0]):
    plt.text(X_norm[i, 0], X_norm[i, 1], str(y[i]), color=plt.cm.Set1(y[i]), 
             fontdict={'weight': 'bold', 'size': 9})
plt.xticks([])
plt.yticks([])
# 时间
title = 't-SNE embedding of the digits (time %.2fs)' % (time() - t0) 
plt.title(title)
plt.show()

实验结果
原始数据显示

可视化结果

训练过程显示

总结

1. 通过对比PCA和t-SNE算法，我了解了线性和非线性降维算法各自的“偏好”。对降维有了更深的了解，且意识到了：在对数据进行可视化前，一定要先分析数据形式或特征，以便能选出更适用的算法。
2. SNE是t-SNE的基础，t-SNE通过使用t-分布改进算法，解决了SNE留下的问题，取得了很好的优化效果。这个改动非常巧妙，我们在平常的学习中不能忽视任何一种基础或经典的模型、算法，它们在很多领域上都可以有很巧妙的应用，而这些应用很可能给我们的算法带来很好的优化效果。我还了解到：有人使用多种树的算法和负采样将t-SNE的时间复杂度降至线性级。未来可以将其作为继续学习的方向。
3. 在阅读sklearn源码和实现手写数字可视化项目的过程中，提高了我对算法结构的理解程度，加强了记忆。同时，自己对数据预处理和代码结构的理解都有了一定程度的提高。之后可以试着进行其他实验，如：将t-SNE算法用到到其他项目的高维数据降维中，以降低任务难度，如：人脸识别。

参考资料

1. Van Der Maaten L., Hinton G.: Visualizing Data using t-SNE(2008). J Mach Learn Res 9, 2579-2605(2008), 85.
2. sklearn v0.20.1 官方文档.