t-SNE完整笔记 (附Python代码)

来源:datakit1 11,969

t-SNE(t-distributed stochastic neighbor embedding)是用于降维的一种机器学习算法,是由 Laurens van der Maaten 和 Geoffrey Hinton在08年提出来。此外,t-SNE 是一种非线性降维算法,非常适用于高维数据降维到2维或者3维,进行可视化。

t-SNE是由SNE(Stochastic Neighbor Embedding, SNE; Hinton and Roweis, 2002)发展而来。我们先介绍SNE的基本原理,之后再扩展到t-SNE。最后再看一下t-SNE的实现以及一些优化。

1.SNE

1.1基本原理

SNE是通过仿射(affinitie)变换将数据点映射到概率分布上,主要包括两个步骤:

  • SNE构建一个高维对象之间的概率分布,使得相似的对象有更高的概率被选择,而不相似的对象有较低的概率被选择。
  • SNE在低维空间里在构建这些点的概率分布,使得这两个概率分布之间尽可能的相似。

我们看到t-SNE模型是非监督的降维,他跟kmeans等不同,他不能通过训练得到一些东西之后再用于其它数据(比如kmeans可以通过训练得到k个点,再用于其它数据集,而t-SNE只能单独的对数据做操作,也就是说他只有fit_transform,而没有fit操作)

1.2 SNE原理推导

SNE是先将欧几里得距离转换为条件概率来表达点与点之间的相似度。具体来说,给定一个N个高维的数据 \(x_1, ... , x_N\)(注意N不是维度), t-SNE首先是计算概率\(p_{ij}\),正比于\(x_i\)和\(x_j\)之间的相似度(这种概率是我们自主构建的),即:

\[{p_ {j \mid i} = \frac{\exp(- \mid \mid x_i -x_j \mid \mid ^2 / (2 \sigma^2_i ))} {\sum_{k \neq i} \exp(- \mid \mid x_i - x_k \mid \mid ^2 / (2 \sigma^2_i))}}\]

这里的有一个参数是\(\sigma_i\),对于不同的点\(x_i\)取值不一样,后续会讨论如何设置。此外设置\(p_{x \mid x}=0\),因为我们关注的是两两之间的相似度。

那对于低维度下的\(y_i\),我们可以指定高斯分布为方差为\(\frac{1}{\sqrt{2}}\),因此它们之间的相似度如下:

\[{q_ {j \mid i} = \frac{\exp(- \mid \mid x_i -x_j \mid \mid ^2)} {\sum_{k \neq i} \exp(- \mid \mid x_i - x_k \mid \mid ^2)}}\]

同样,设定\(q_{i \mid i} = 0\).

如果降维的效果比较好,局部特征保留完整,那么 \(p_{i \mid j} = q_{i \mid j}\), 因此我们优化两个分布之间的距离-KL散度(Kullback-Leibler divergences),那么目标函数(cost function)如下:

\[C = \sum_i KL(P_i \mid \mid Q_i) = \sum_i \sum_j p_{j \mid i} \log \frac{p_{j \mid i}}{q_{j \mid i}}\]

这里的\(P_i\)表示了给定点\(x_i\)下,其他所有数据点的条件概率分布。需要注意的是KL散度具有不对称性,在低维映射中不同的距离对应的惩罚权重是不同的,具体来说:距离较远的两个点来表达距离较近的两个点会产生更大的cost,相反,用较近的两个点来表达较远的两个点产生的cost相对较小(注意:类似于回归容易受异常值影响,但效果相反)。即用较小的 \(q_{j \mid i}=0.2\) 来建模较大的 \(p_{j \mid i}=0.8\), cost=\(p \log(\frac{p}{q})\)=1.11,同样用较大的\(q_{j \mid i}=0.8\)来建模较大的\(p_{j \mid i}=0.2\), cost=-0.277, 因此,SNE会倾向于保留数据中的局部特征

思考:了解了基本思路之后,你会怎么选择\(\sigma\),固定初始化?

下面我们开始正式的推导SNE。首先不同的点具有不同的\(\sigma_i\),\(P_i\)的熵(entropy)会随着\(\sigma_i\)的增加而增加。SNE使用困惑度(perplexity)的概念,用二分搜索的方式来寻找一个最佳的\(\sigma\)。其中困惑度指:

\[Perp(P_i) = 2^{H(P_i)}\]

这里的\(H(P_i)\)是\(P_i\)的熵,即:

\[H(P_i) = -\sum_j p_{j \mid i} \log_2 p_{j \mid i}\]

困惑度可以解释为一个点附近的有效近邻点个数。SNE对困惑度的调整比较有鲁棒性,通常选择5-50之间,给定之后,使用二分搜索的方式寻找合适的\(\sigma\)

那么核心问题是如何求解梯度了,目标函数等价于\(\sum \sum - p log(q)\)这个式子与softmax非常的类似,我们知道softmax的目标函数是\(\sum -y \log p\),对应的梯度是\(y - p\)(注:这里的softmax中y表示label,p表示预估值)。 同样我们可以推导SNE的目标函数中的i在j下的条件概率情况的梯度是\(2(p_{i \mid j}-q_{i \mid j})(y_i-y_j)\), 同样j在i下的条件概率的梯度是\(2(p_{j \mid i}-q_{j \mid i})(y_i-y_j)\), 最后得到完整的梯度公式如下:

\[\frac{\delta C}{\delta y_i} = 2 \sum_j (p_{j \mid i} - q_{j \mid i} + p_{i \mid j} - q_{i \mid j})(y_i - y_j)\]

在初始化中,可以用较小的\(\sigma\)下的高斯分布来进行初始化。为了加速优化过程和避免陷入局部最优解,梯度中需要使用一个相对较大的动量(momentum)。即参数更新中除了当前的梯度,还要引入之前的梯度累加的指数衰减项,如下:

\[Y^{(t)} = Y^{(t-1)} + \eta \frac{\delta C}{\delta Y} + \alpha(t)(Y^{(t-1)} - Y^{(t-2)})\]

这里的\(Y^{(t)}\)表示迭代t次的解,\(\eta\)表示学习速率,\(\alpha(t)\)表示迭代t次的动量。

此外,在初始优化的阶段,每次迭代中可以引入一些高斯噪声,之后像模拟退火一样逐渐减小该噪声,可以用来避免陷入局部最优解。因此,SNE在选择高斯噪声,以及学习速率,什么时候开始衰减,动量选择等等超参数上,需要跑多次优化才可以。

思考:SNE有哪些不足? 面对SNE的不足,你会做什么改进?

2.t-SNE

尽管SNE提供了很好的可视化方法,但是他很难优化,而且存在”crowding problem”(拥挤问题)。后续中,Hinton等人又提出了t-SNE的方法。与SNE不同,主要如下:

  • 使用对称版的SNE,简化梯度公式
  • 低维空间下,使用t分布替代高斯分布表达两点之间的相似度

t-SNE在低维空间下使用更重长尾分布的t分布来避免crowding问题和优化问题。在这里,首先介绍一下对称版的SNE,之后介绍crowding问题,之后再介绍t-SNE。

2.1 Symmetric SNE

优化\(p_{i \mid j}\)和\(q_{i \mid j}\)的KL散度的一种替换思路是,使用联合概率分布来替换条件概率分布,即P是高维空间里各个点的联合概率分布,Q是低维空间下的,目标函数为:

\[C = KL(P \mid \mid Q) = \sum_i \sum_j p_{i,j} \log \frac{p_{ij}}{q_{ij}}\]

这里的\(p_{ii}\),\(q_{ii}\)为0,我们将这种SNE称之为symmetric SNE(对称SNE),因为他假设了对于任意i,\(p_{ij} = p_{ji}, q_{ij} = q_{ji}\),因此概率分布可以改写为:

\[p_{ij} = \frac{\exp(- \mid \mid x_i - x_j \mid \mid ^2 / 2\sigma^2)}{\sum_{k \neq l} \exp(- \mid \mid x_k-x_l \mid \mid ^2 / 2\sigma^2)} \ \ \ \ q_{ij} = \frac{\exp(- \mid \mid y_i - y_j \mid \mid ^2)}{\sum_{k \neq l} \exp(- \mid \mid y_k-y_l \mid \mid ^2)}\]

这种表达方式,使得整体简洁了很多。但是会引入异常值的问题。比如\(x_i\)是异常值,那么\(\mid \mid x_i - x_j \mid \mid ^2\)会很大,对应的所有的j, \(p_{ij}\)都会很小(之前是仅在\(x_i\)下很小),导致低维映射下的\(y_i\)对cost影响很小。

思考: 对于异常值,你会做什么改进?\(p_i\)表示什么?

为了解决这个问题,我们将联合概率分布定义修正为: \(p_{ij} = \frac{p_{i \mid j} + p_{j \mid i}}{2}\), 这保证了\(\sum_j p_{ij} \gt \frac{1}{2n}\), 使得每个点对于cost都会有一定的贡献。对称SNE的最大优点是梯度计算变得简单了,如下:

\(\frac{\delta C}{\delta y_i} = 4 \sum_j (p_{ij} - q_{ij})(y_i - y_j)\)

实验中,发现对称SNE能够产生和SNE一样好的结果,有时甚至略好一点。

2.2 Crowding问题

拥挤问题就是说各个簇聚集在一起,无法区分。比如有一种情况,高维度数据在降维到10维下,可以有很好的表达,但是降维到两维后无法得到可信映射,比如降维如10维中有11个点之间两两等距离的,在二维下就无法得到可信的映射结果(最多3个点)。
进一步的说明,假设一个以数据点\(x_i\)为中心,半径为r的m维球(三维空间就是球),其体积是按\(r^m\)增长的,假设数据点是在m维球中均匀分布的,我们来看看其他数据点与\(x_i\)的距离随维度增大而产生的变化。

t-SNE完整笔记 (附Python代码)-图片1

从上图可以看到,随着维度的增大,大部分数据点都聚集在m维球的表面附近,与点\(x_i\)的距离分布极不均衡。如果直接将这种距离关系保留到低维,就会出现拥挤问题。

怎么解决crowding问题呢?

Cook et al.(2007) 提出一种slight repulsion的方式,在基线概率分布(uniform background)中引入一个较小的混合因子\(\rho\),这样\(q_{ij}\)就永远不会小于\(\frac{2 \rho}{n(n-1)}\) (因为一共了n(n-1)个pairs),这样在高维空间中比较远的两个点之间的\(q_{ij}\)总是会比\(p_{ij}\)大一点。这种称之为UNI-SNE,效果通常比标准的SNE要好。优化UNI-SNE的方法是先让\(\rho\)为0,使用标准的SNE优化,之后用模拟退火的方法的时候,再慢慢增加\(\rho\).
直接优化UNI-SNE是不行的(即一开始\(\rho\)不为0),因为距离较远的两个点基本是一样的\(q_{ij}\)(等于基线分布), 即使\(p_{ij}\)很大,一些距离变化很难在\(q_{ij}\)中产生作用。也就是说优化中刚开始距离较远的两个聚类点,后续就无法再把他们拉近了。

2.3 t-SNE

对称SNE实际上在高维度下
另外一种减轻”拥挤问题”的方法:在高维空间下,在高维空间下我们使用高斯分布将距离转换为概率分布,在低维空间下,我们使用更加偏重长尾分布的方式来将距离转换为概率分布,使得高维度下中低等的距离在映射后能够有一个较大的距离。

 

t-SNE完整笔记 (附Python代码)-图片2

我们对比一下高斯分布和t分布(如上图,code见probability/distribution.md), t分布受异常值影响更小,拟合结果更为合理,较好的捕获了数据的整体特征。

使用了t分布之后的q变化,如下:

\[q_{ij} = \frac{(1 + \mid \mid y_i -y_j \mid \mid ^2)^{-1}}{\sum_{k \neq l} (1 + \mid \mid y_i -y_j \mid \mid ^2)^{-1}}\]

此外,t分布是无限多个高斯分布的叠加,计算上不是指数的,会方便很多。优化的梯度如下:

\[\frac{\delta C}{\delta y_i} = 4 \sum_j(p_{ij}-q_{ij})(y_i-y_j)(1+ \mid \mid y_i-y_j \mid \mid ^2)^{-1}\]

t-SNE完整笔记 (附Python代码)-图片3

t-sne的有效性,也可以从上图中看到:横轴表示距离,纵轴表示相似度, 可以看到,对于较大相似度的点,t分布在低维空间中的距离需要稍小一点;而对于低相似度的点,t分布在低维空间中的距离需要更远。这恰好满足了我们的需求,即同一簇内的点(距离较近)聚合的更紧密,不同簇之间的点(距离较远)更加疏远。

总结一下,t-SNE的梯度更新有两大优势:

  • 对于不相似的点,用一个较小的距离会产生较大的梯度来让这些点排斥开来。
  • 这种排斥又不会无限大(梯度中分母),避免不相似的点距离太远。

2.4 算法过程

算法详细过程如下:

  • Data: \(X = {x_1, ... , x_n}\)
  • 计算cost function的参数:困惑度Perp
  • 优化参数: 设置迭代次数T, 学习速率\(\eta\), 动量\(\alpha(t)\)
  • 目标结果是低维数据表示 \(Y^T = {y_1, ... , y_n}\)
  • 开始优化
    • 计算在给定Perp下的条件概率\(p_{j \mid i}\)(参见上面公式)
    • 令 \(p_{ij} = \frac{p_{j \mid i} + p_{i \mid j}}{2n}\)
    • 用 \(N(0, 10^{-4}I)\) 随机初始化 Y
    • 迭代,从 t = 1 到 T, 做如下操作:
      • 计算低维度下的 \(q_{ij}\)(参见上面的公式)
      • 计算梯度(参见上面的公式)
      • 更新 \(Y^{t} = Y^{t-1} + \eta \frac{dC}{dY} + \alpha(t)(Y^{t-1} - Y^{t-2})\)
    • 结束
  • 结束

优化过程中可以尝试的两个trick:

  • 提前压缩(early compression):开始初始化的时候,各个点要离得近一点。这样小的距离,方便各个聚类中心的移动。可以通过引入L2正则项(距离的平方和)来实现。
  • 提前夸大(early exaggeration):在开始优化阶段,\(p_{ij}\)乘以一个大于1的数进行扩大,来避免因为\(q_{ij}\)太小导致优化太慢的问题。比如前50次迭代,\(p_{ij}\)乘以4

优化的过程动态图如下:

t-SNE完整笔记 (附Python代码)-图片4

2.5 不足

主要不足有四个:

  • 主要用于可视化,很难用于其他目的。比如测试集合降维,因为他没有显式的预估部分,不能在测试集合直接降维;比如降维到10维,因为t分布偏重长尾,1个自由度的t分布很难保存好局部特征,可能需要设置成更高的自由度。
  • t-SNE倾向于保存局部特征,对于本征维数(intrinsic dimensionality)本身就很高的数据集,是不可能完整的映射到2-3维的空间
  • t-SNE没有唯一最优解,且没有预估部分。如果想要做预估,可以考虑降维之后,再构建一个回归方程之类的模型去做。但是要注意,t-sne中距离本身是没有意义,都是概率分布问题。
  • 训练太慢。有很多基于树的算法在t-sne上做一些改进

3.变种

后续有机会补充。

  • multiple maps of t-SNE
  • parametric t-SNE
  • Visualizing Large-scale and High-dimensional Data

4.参考文档

  • Maaten, L., & Hinton, G. (2008). Visualizing data using t-SNE. Journal of Machine Learning Research.

5. 代码

文中的插图绘制:

  1. # coding:utf-8
  2.  
  3. import numpy as np
  4. from numpy.linalg import norm
  5. from matplotlib import pyplot as plt
  6. plt.style.use('ggplot')
  7.  
  8. def sne_crowding():
  9. npoints = 1000 # 抽取1000个m维球内均匀分布的点
  10. plt.figure(figsize=(20, 5))
  11. for i, m in enumerate((2, 3, 5, 8)):
  12. # 这里模拟m维球中的均匀分布用到了拒绝采样,
  13. # 即先生成m维立方中的均匀分布,再剔除m维球外部的点
  14. accepts = []
  15. while len(accepts) < 1000:
  16. points = np.random.rand(500, m)
  17. accepts.extend([d for d in norm(points, axis=1)
  18. if d <= 1.0]) # 拒绝采样
  19. accepts = accepts[:npoints]
  20. ax = plt.subplot(1, 4, i+1)
  21. if i == 0:
  22. ax.set_ylabel('count')
  23. if i == 2:
  24.  
  25. ax.set_xlabel('distance')
  26. ax.hist(accepts, bins=np.linspace(0., 1., 50))
  27. ax.set_title('m=%s' %m)
  28. plt.savefig("./images/sne_crowding.png")
  29.  
  30. x = np.linspace(0, 4, 100)
  31. ta = 1 / (1 + np.square(x))
  32. tb = np.sum(ta) - 1
  33. qa = np.exp(-np.square(x))
  34. qb = np.sum(qa) - 1
  35.  
  36. def sne_norm_t_dist_cost():
  37. plt.figure(figsize=(8, 5))
  38. plt.plot(qa/qb, c="b", label="normal-dist")
  39. plt.plot(ta/tb, c="g", label="t-dist")
  40. plt.plot((0, 20), (0.025, 0.025), 'r--')
  41. plt.text(10, 0.022, r'$q_{ij}$')
  42. plt.text(20, 0.026, r'$p_{ij}$')
  43.  
  44. plt.plot((0, 55), (0.005, 0.005), 'r--')
  45. plt.text(36, 0.003, r'$q_{ij}$')
  46. plt.text(55, 0.007, r'$p_{ij}$')
  47.  
  48. plt.title("probability of distance")
  49. plt.xlabel("distance")
  50. plt.ylabel("probability")
  51. plt.legend()
  52. plt.savefig("./images/sne_norm_t_dist_cost.png")
  53.  
  54. if __name__ == '__main__':
  55. sne_crowding()
  56. sne_norm_t_dist_cost()

附录一下t-sne的完整代码实现:

  1. # coding utf-8
  2. '''
  3. 代码参考了作者Laurens van der Maaten的开放出的t-sne代码, 并没有用类进行实现,主要是优化了计算的实现
  4. '''
  5. import numpy as np
  6.  
  7.  
  8. def cal_pairwise_dist(x):
  9. '''计算pairwise 距离, x是matrix
  10. (a-b)^2 = a^w + b^2 - 2*a*b
  11. '''
  12. sum_x = np.sum(np.square(x), 1)
  13. dist = np.add(np.add(-2 * np.dot(x, x.T), sum_x).T, sum_x)
  14. return dist
  15.  
  16.  
  17. def cal_perplexity(dist, idx=0, beta=1.0):
  18. '''计算perplexity, D是距离向量,
  19. idx指dist中自己与自己距离的位置,beta是高斯分布参数
  20. 这里的perp仅计算了熵,方便计算
  21. '''
  22. prob = np.exp(-dist * beta)
  23. # 设置自身prob为0
  24. prob[idx] = 0
  25. sum_prob = np.sum(prob)
  26. perp = np.log(sum_prob) + beta * np.sum(dist * prob) / sum_prob
  27. prob /= sum_prob
  28. return perp, prob
  29.  
  30.  
  31. def seach_prob(x, tol=1e-5, perplexity=30.0):
  32. '''二分搜索寻找beta,并计算pairwise的prob
  33. '''
  34.  
  35. # 初始化参数
  36. print("Computing pairwise distances...")
  37. (n, d) = x.shape
  38. dist = cal_pairwise_dist(x)
  39. pair_prob = np.zeros((n, n))
  40. beta = np.ones((n, 1))
  41. # 取log,方便后续计算
  42. base_perp = np.log(perplexity)
  43.  
  44. for i in range(n):
  45. if i % 500 == 0:
  46. print("Computing pair_prob for point %s of %s ..." %(i,n))
  47.  
  48. betamin = -np.inf
  49. betamax = np.inf
  50. perp, this_prob = cal_perplexity(dist[i], i, beta[i])
  51.  
  52. # 二分搜索,寻找最佳sigma下的prob
  53. perp_diff = perp - base_perp
  54. tries = 0
  55. while np.abs(perp_diff) > tol and tries < 50:
  56. if perp_diff > 0:
  57. betamin = beta[i].copy()
  58. if betamax == np.inf or betamax == -np.inf:
  59. beta[i] = beta[i] * 2
  60. else:
  61. beta[i] = (beta[i] + betamax) / 2
  62. else:
  63. betamax = beta[i].copy()
  64. if betamin == np.inf or betamin == -np.inf:
  65. beta[i] = beta[i] / 2
  66. else:
  67. beta[i] = (beta[i] + betamin) / 2
  68.  
  69. # 更新perb,prob值
  70. perp, this_prob = cal_perplexity(dist[i], i, beta[i])
  71. perp_diff = perp - base_perp
  72. tries = tries + 1
  73. # 记录prob值
  74. pair_prob[i,] = this_prob
  75. print("Mean value of sigma: ", np.mean(np.sqrt(1 / beta)))
  76. return pair_prob
  77.  
  78.  
  79. def pca(x, no_dims = 50):
  80. ''' PCA算法
  81. 使用PCA先进行预降维
  82. '''
  83. print("Preprocessing the data using PCA...")
  84. (n, d) = x.shape
  85. x = x - np.tile(np.mean(x, 0), (n, 1))
  86. l, M = np.linalg.eig(np.dot(x.T, x))
  87. y = np.dot(x, M[:,0:no_dims])
  88. return y
  89.  
  90.  
  91. def tsne(x, no_dims=2, initial_dims=50, perplexity=30.0, max_iter=1000):
  92. """Runs t-SNE on the dataset in the NxD array x
  93. to reduce its dimensionality to no_dims dimensions.
  94. The syntaxis of the function is Y = tsne.tsne(x, no_dims, perplexity),
  95. where x is an NxD NumPy array.
  96. """
  97.  
  98. # Check inputs
  99. if isinstance(no_dims, float):
  100. print("Error: array x should have type float.")
  101. return -1
  102. if round(no_dims) != no_dims:
  103. print("Error: number of dimensions should be an integer.")
  104. return -1
  105.  
  106. # 初始化参数和变量
  107. x = pca(x, initial_dims).real
  108. (n, d) = x.shape
  109. initial_momentum = 0.5
  110. final_momentum = 0.8
  111. eta = 500
  112. min_gain = 0.01
  113. y = np.random.randn(n, no_dims)
  114. dy = np.zeros((n, no_dims))
  115. iy = np.zeros((n, no_dims))
  116. gains = np.ones((n, no_dims))
  117.  
  118. # 对称化
  119. P = seach_prob(x, 1e-5, perplexity)
  120. P = P + np.transpose(P)
  121. P = P / np.sum(P)
  122. # early exaggeration
  123. P = P * 4
  124. P = np.maximum(P, 1e-12)
  125.  
  126. # Run iterations
  127. for iter in range(max_iter):
  128. # Compute pairwise affinities
  129. sum_y = np.sum(np.square(y), 1)
  130. num = 1 / (1 + np.add(np.add(-2 * np.dot(y, y.T), sum_y).T, sum_y))
  131. num[range(n), range(n)] = 0
  132. Q = num / np.sum(num)
  133. Q = np.maximum(Q, 1e-12)
  134.  
  135. # Compute gradient
  136. PQ = P - Q
  137. for i in range(n):
  138. dy[i,:] = np.sum(np.tile(PQ[:,i] * num[:,i], (no_dims, 1)).T * (y[i,:] - y), 0)
  139.  
  140. # Perform the update
  141. if iter < 20:
  142. momentum = initial_momentum
  143. else:
  144. momentum = final_momentum
  145. gains = (gains + 0.2) * ((dy > 0) != (iy > 0)) + (gains * 0.8) * ((dy > 0) == (iy > 0))
  146. gains[gains < min_gain] = min_gain
  147. iy = momentum * iy - eta * (gains * dy)
  148. y = y + iy
  149. y = y - np.tile(np.mean(y, 0), (n, 1))
  150. # Compute current value of cost function
  151. if (iter + 1) % 100 == 0:
  152. if iter > 100:
  153. C = np.sum(P * np.log(P / Q))
  154. else:
  155. C = np.sum( P/4 * np.log( P/4 / Q))
  156. print("Iteration ", (iter + 1), ": error is ", C)
  157. # Stop lying about P-values
  158. if iter == 100:
  159. P = P / 4
  160. print("finished training!")
  161. return y
  162.  
  163.  
  164. if __name__ == "__main__":
  165. # Run Y = tsne.tsne(X, no_dims, perplexity) to perform t-SNE on your dataset.
  166. X = np.loadtxt("mnist2500_X.txt")
  167. labels = np.loadtxt("mnist2500_labels.txt")
  168. Y = tsne(X, 2, 50, 20.0)
  169. from matplotlib import pyplot as plt
  170. plt.scatter(Y[:,0], Y[:,1], 20, labels)
  171. plt.show()

 

 

    • 探索者 0

      t-sne怎么可视化自己的图像数据啊

    发表评论

    匿名网友

    拖动滑块以完成验证
    加载失败