降维打击-观察高维世界

我们知道1维的世界最多有2个点互相距离相同，2维的世界最多有3个点，3维有4个。那么4维呢？数学上的直观可以告诉我们是5个，并且可以证明出来。但是，你能想象这种情况吗？–不能，因为我们天生如此。

你可能觉得这个没有什么意思–数学做的足够好了。但是，对于一些复杂而又不平凡的情况，不是用简单的数学可以解决的。比如如果我们把手写数字的图片想象成一个向量–向量每一项代表图像中每一点的像素值，那么即使对于mnist这种最简单手写数字图片的数据集之一，向量的维度也可以达到 28*28 = 784。

这是MNIST数据集的几个例子：

我们这样表示向量：

先把它表示为一个矩阵（应该28*28的，没有画全）

然后摊平即可

整个数据集构成了784维空间中向量（或者点）的集合，并且难以用一个数学公式描述出它的特性。

这个时候，为了直观表示出高维中数据的“样子”，我们需要降维–即把高维空间投射到3维以内。

降维就其本身来说是非常平凡的东西。最简单的方法是“投影“，就是沿着一个轴把一个方向的数据的变化消除（比如把3维的一个方向上的东西压平到平面上成为2维），或者说将向量的某个分量全部置0。这样每次投影可以减少一个纬度。但是对于784维，至少需要781次投影，这样几乎确信自己很难发现什么（一般来说4维以上的进行投影就基本无效了）。更麻烦的问题是，我们不知道沿着哪个方向投影更好。

我们需要更加漂亮的方式来降低维度，同时维持高维度的”感觉“。

首先我们有一种数学方法。这种方式叫PCA（主成分分析）。PCA的数学定义是：一个正交化线性变换，把数据变换到一个新的坐标系统中，使得这一数据的任何投影的第一大方差在第一个坐标（称为第一主成分）上，第二大方差在第二个坐标（第二主成分）上，依次类推。从信息论角度考虑主成分分析在降维上是最优的（证明略，主要通过信息熵），它解决了投影问题。但是实际过程中我们可以发现它做的并不总是令人满意，比如下图是它对MNIST的784维空间的一个投影(选了一部分数据，更多数据就糟到看不清了)，不同颜色代表不同数字。我们可以看到很多数字“纠缠”到一起（除了左下角的1和右下角的0），和我们对数字容易区分的认知相差很远。

你可能感到不理解：为什么信息论上面可以证明最优的在实际中并不是最优？这是因为，“最优”的前提是对所有情况发生概率一致的基础上得出的，对有先验知识的东西并不成立。这个再次凸显了 “A more general theory is more useless” 的道理。

比如，PCA 处理下面噪声图片的方式和处理上面的没有什么不同。它不关心数据本身，而仅仅是一个方向上面的方差，所以是一个“过大”的理论（但对于方差显著的数据很有效）。

那么我们有没有其他方法呢？

注：

1. PCA 常被称为线性降维方法，因为它对数据进行了线性变换。另外一种常用的线性降维方法是 LDA (Linear Discriminant Analysis)，下面正文要谈的都是非线性的。
2. Laplacian Eigenmaps（由求图的 Laplace 矩阵的特征值而得名）也是一种线性降维方法。
3. LLE (Locally Linear Embedding)， 局部线性嵌入算法，是一种非线性算法。LLE算法认为每一个数据点都可以由其近邻点的线性加权组合构造得到。
4. PCA, LDA, Laplacian Eigenmaps, LLE 都属于谱（spectral）方法（涉及矩阵特征向量的方法）。
5. 线性方法计算效率是非常高的，非线性算法对于大规模数据需要考虑随机游走(random walk)等优化，否则计算代价非常高。

距离表示

仔细想想我们会发现，表示出数据的空间结构并不是一件完全没有规律可遵循的事情。首先，旋转、平移和翻转不会改变数据本身的结构，对于高维空间也是如此。并且结构还明显受点的个数的影响–3个点对于高于2维的空间是没有特殊意义的，我们需要一种东西可以保证不受维数增加的影响。

因此，我们开始考虑用点之间的距离（这里使用欧氏距离）代替点的坐标 – 距离不受旋转和平移的影响，并且只和点数有关。

但是，距离表示会不会丢失某些结构信息呢？不会，我们以2维为例构造式的说明一下：

只有一个距离的时候当然随便什么角度和位置排放点都不会影响结构：

假设不是退化到一维的情况（即存在不共线的点），它们三个点的三个边一定构成唯一形状的三角形

下面任何一个点将和这三个点构成3条边，也就是3个方程，足以确定2维坐标的两个分量。否则那个点不在二维平面上。

其他维度类似。所以各点间的距离确实完整保留了空间结构，同时又去除了不必要的转动和平移翻转。同时，距离本身不直接包含维度信息，不会受到冗余维度的影响。

但是，正是因为距离本身不直接包含维度，所以我们似乎没有什么特别明确的算法。因此我们退一步，把问题转化为优化问题的形式，以备我们有更灵活的选择：

$$\min_{D_Y} E(D_X, D_Y)$$

其中 $D_X, D_Y$ 分别是原空间和当前空间中的边。E（能量函数，损失函数）是一种评估它们“结构相似”的函数。

然后我们试图寻找一个策略，使得它们“最为相似”。这样比直接求 $D_Y$ 有更多选择的余地。

但困难之处还在于，你必须保持低维度空间的距离是合法的（比如三维中一共5点10条边，长度都相等是无解的）。所以我们需要间接的方法：

设 $X$ 是原空间的点（向量），$Y$ 是低维空间的点，$D_X, D_Y$ 分别是它们对应的边。我们寻求一个 $Y$ 最小化 $E$：

$$\min_Y E(D_X, D_Y)$$

这样一来在于，我们将难以直接求解 $Y$，因为 $Y$ 的距离是关于Y的相当复杂的函数（带 $n*dim$ 个自变量，$n*(n-1)/2$ 个因变量，还有平方，根号），再加上 $E$ 的与 $D_X$ 相关的复杂运算，使得我们基本不可能有个可以接受的公式直接根据 $E$ 求得 $Y$ 。

不过不用害怕。这类问题的解法早已熟知，那就是采用各种数值优化算法，对于多数使用梯度下降就足够了，这里暂时不做更多介绍。

所以下面的关建是如何选择 $E$ （能量函数，损失函数）

基于距离相似的 E

MDS (Multidimensional scaling)

首先我们想到的是，显然一般情况下我们没有办法让 $D_Y$和 $D_X$ 完全相同，否则低维空间无法容纳这样的 $D_Y$。这让我们想到了直线拟合：对于大部分实际数据都无法真正过一条直线，因此我们有最小二乘法，即使得直线和数据点之间的平方误差最小（这是少数我们能够有解析解的情况，对于不能转化为直线的一些复杂的情况我们只能采用优化的方式）。按照最小二乘法的思想，我们很容易设计以下的 $E$（按照求和约定省略了部分标记，但保留了求和符号。下同）：

$$E_{MDS} = \frac 1{\sum D_{X_i}^2}\sum (D_{X_i}-D_{Y_i})^2$$

前面的系数是正则化系数，目的是为了保证 $E$ 不受距离加倍的影响（同时把所有距离乘上一个正因子不应该改变结果）

这种方法属于所谓 Multidimensional scaling

结果如下：

可以看到一些团簇，但是仍然不能让人满意。

我们注意到，问题在于一些理应分开的点变的像一团粥一样，Multidimensional scaling 同等的对待所有的距离，这是和PCA一样的问题。因此我们希望能够更加强调一些局部结构，让它们更加紧密一些。

Sammon Mapping

我们稍稍修改一下上面的E，把一部分 Dx 移入到后面参加运算。

$$E_{sammon} = \frac 1{\sum D_{X_i}}\sum \frac{(D_{X_i}-D_{Y_i})^2}{D_{X_i}}$$

这就是所谓 Sammon Mapping。它的效果是强调更短的边（在分母上，越小就意味着上面的平方项会带来更大的损失，从而约束上面的平方项更加小）

然而它的效果似乎并不是那么明显：

然而大量的事实确实证实 Sammon Mapping 在比较简单的情况下很有效。为什么这里就不行了呢？

因为，这是784维的真实数据，这么高的维数几乎是种诅咒。

Isomap

之前我们采用距离是欧几里德距离，它是平坦空间的最短距离。

Isomap 把距离推广为了“测地距离”。“测地距离”是两点之间最短路径长度（离散空间）或者某种加权的测度（连续的空间）。一般Isomap会先把问题转化为离散的处理，这里不举例子了。

附：维数灾难

我们来做一个有趣的实验：

对于k维，立方体的顶点数是 $2^k$，原因是在单位正交坐标系中，k维立方体的顶点可以用k个0/1构成的坐标描述，比如一个正方形顶点为 $(0,0),(0,1),(1,0),(1,1)$ 。这种广义的正方体有很多奇妙的性质，比如相邻的点仅有一个分量不同，可以从一个点出发不重复沿着边走遍所有点等等。

然后我们计算每个点到其他点的距离的方差（写个程序就行了，注意一些技巧）。结论是方差随着维数增多缓慢减小，并稳定在 0.125 左右。

这意味着什么？我们知道，一个分布的方差是由数据分布的“形状”决定的，越弥散的分布有越大的方差。方差不变意味着分布形状基本不变，但是引入了过多的数据，这样一来一个分布局部区间必然包含了更多的数据，也就是更多距离相似的边，虽然立方体这个数学上的模型没有变化。

比如，对于2维数据，对角线和边长相差 $\sqrt{2}-1\approx 0.414$，但是对于1000维的数据，存在 499,500 个长 $\sqrt{998}$ 和 166,167,000 个长 $\sqrt{997}$ 的边，而它们的长度仅仅相差 0.016，而相对相差只有 0.0005。

这个简单的例子就是所谓“维数灾难”的一种形式，即随着维数增大，边之间的“分辨率”会降低。这可能就是 Sammon Mapping 没有取得很好结果的原因，即它除的分母太相似了，没有过多的区分度。

基于物理体系模拟的 E

之前我们公式中的主项，$\sum (D_{X_i}-D_{Y_i})^2$ 是否有什么物理意义呢？

如果把现有空间中的边看作弹簧，那么上面的公式就是劲度系数为1的弹簧体系的弹性势能。其中对于弹簧 $D_{Y_i}$，$D_{X_i}$ 是它的平衡位置。

我们观察到 Sammon Mapping 试图靠拢相邻点的做法失败了，那么我们是不是可以考虑推远不相邻的点？

做法很简单，就是引入一个势能场，比如让所有点都带上单位同种电荷。

这样一来能量函数（名副其实）就是

$$E_{ForceField} = \sum\frac1{D_{Y_i}} + \sum k(D_{X_i}-D_{Y_i})^2$$

其中k是劲度系数。(注：静电势能选择反比形式，是因为考虑到三维中平方受力的稳定性较好。对于不同维度，维持静电场为有源无旋场根据高斯定理会导致不同的势能形式。)

效果如下：

嗯，好像没有什么改进。原因同 Sammon Mapping，前面的静电势能项因为大家距离都差不多而没有起到实质作用。

那么我们有没有改进方案呢？在我们脑海中，一个弹簧和斥力的体系会自己根据受力平衡“布局”到一个好的位置，这也是布图的一种常见算法。

它的能量函数是 (k是是劲度系数，$D_0$ 是弹簧的平衡位置)：

$$E_{ForceField} = \sum\frac1{D_{Y_i}} + \sum_{Y_i^* \in edges} k(D_{Y_i^*}-D_0)^2$$

你可能感到不可思议的是，$D_X$ 到哪里去了？我们怎么可能不需要原来的距离就得到结果呢？

答案是，$D_X$ 用来确定edges，也就是保留哪些边。

注意，之前我们是把所有的边都纳入到能量体系当中，而现在我们需要根据 $D_X$ 选择一些作为“弹簧”。

选择的方法是使用古老的 KNN (k-nestest neighbors)算法，也就是选择每个点最近的k个边。

这样我们就有

$$E_{KNN-ForceField} = \sum\frac1{D_{Y_i}} + \sum_{Y_i^* \in knnedges (D_X)} k(D_{Y_i^*}-D_0)^2$$

knn 的 k = 4 的时候，

k = 3

嗯，总体上要好很多。主要贡献应该归于用KNN预处理了 $D_X$，KNN仅仅选择最短的几个边，即使各个边长度差距不大，这样大大增加了区分度。之前的方法主要通过把 $D_{X_i}$ 放在分母，区分度是不够的。

但是引入KNN的方法把问题变成了黑盒：k 对数据本身的意义是什么？这个不是很清楚。另外边要不选中，要不不选中，似乎也是一种很极端的做法。

基于概率分布的 E

我们希望有一种算法，能够一定程度上克服使用 KNN 的缺陷。

一种方法是用概率分布取代“是否”这种决定式的方法。这类方法称为 SNE （随机邻居嵌入）。

首先，它将距离转换为概率分布。概率分布中概率比较大的表示更有可能“靠在一起”。

附：将有限个连续数量值转换为概率分布

首先，由于数量值是连续的（比如实数），所以我们需要一个连续的概率分布的概率密度函数(pdf)。

我们直接先把概率密度函数作用上去

$$p_1,p_2,p_3...p_n = f(d_1),f(d_2),f(d_3)...f(d_n)$$

这个时候得到了一些离散的概率值，因为数量值的个数是有限的。

接下来我们归一化这些概率值：

$$p_i^* = \frac{p_i}{\sum p_j}$$

现在 $p^*$ 变成了一个标准的离散的概率分布，它是连续分布的一种近似。

$\square$

然后我们评测原来和现在概率的“相似性”。

对于有重叠的概率分布，一种非常有效而且有很强理论支持的用于表征相似的测度叫 KL divergence （KL散度），一般用记号 ${\rm KL}\left(\cdot|\cdot\right)$ 表示。它也叫相对熵。

所以形式表示为

$$E_{SNE} = {\rm KL}\left(K_X(D_X)\|K_Y(D_Y)\right)$$

符号意义同上文。其中K表示某种映射，将边集合转换为概率分布（将有限个连续数量值转换为概率分布）。

SNE - Part I 将原始高维空间中的每个点相关的距离（邻居距离）转化为概率分布

我们注意到，先前 Sammon Mapping 这些效果不好的原因是对于距离的区分度不是特别好。因此我们考虑通过控制分布的形状来控制区分度。

我们选择的是均值为0的高斯分布（正态分布）。

形式如下：

$$N(0, \sigma) = \frac1{\sqrt{2\pi}\sigma}\exp\left(-\frac{x^2}{2\sigma^2}\right)$$

高斯分布是在方差相等的分布中信息熵最大的，其信息熵为 $\ln\left(\sigma\sqrt{2\pi e}\right)$。在分布未知的情况下根据中心极限定理我们非常有理由选择 Gaussian distribution，因为它是对未知分布的良好近似。

然后我们把这种分布应用到每个点到其他点的距离，得到一组归一化分布 $g^{(i)}_\sigma$，其中 i 代表第 i 个点，$\sigma$ 是所选的 gaussian 分布的方差 。（注意，这里是每个点而言的分布，而不是对所有边的分布）

现在问题是，我们如何选择 $\sigma$ ？

我们用信息熵评价分布的区分度（即 ${\rm H}(g^{(i)}_\sigma)$ ），熵越高则区分度越低，来控制 $\sigma$ 的取值，则可以证明 ${\rm H}(g^{(i)}_\sigma)$ 是关于 $\sigma$ 的单调递增函数。

这个结论并不奇怪。$\sigma$ 代表方差，方差越小则分布越“尖锐”，区分度越高；否则分布越平坦，区分度越低。

但是直接指定信息熵看上去不像是一个良好的做法。因为对于连续分布本身，

$${\rm H}\left(N(0,\sigma)\right) = \ln\left(\sigma\sqrt{2\pi e}\right)$$

即

$$\sigma \propto e^{H(0,\sigma)}$$

这意味着 $\sigma$ 对信息量过于敏感，不便于我们调节。

我们定义 $Perplexity = 2^{H}$，则 $\sigma \propto Perplexity$，相应的离散距离分布也大致满足这个关系。

（这里我们不区分自然对数，2对数和常用对数，因为仅仅影响熵的“单位”而已）

所以，总结一下，就是

$$Find\ \sigma, s.t. {\rm H}(g^{(i)}_\sigma) = \ln (Perplexity)$$

寻找我们使用二分法即可（由于 Perplexity 随着 $\sigma$ 单调增），因为函数太过复杂不适合使用牛顿迭代法等需要导数的方式。

SNE - Part II 将原始高维空间中邻居距离转化为所有边距离的概率分布

到这里仍然没有结束。因为我们仅仅求得关于每个点的邻居距离的概率分布（对有n个点的数据，就有n个分布），目的是控制距离的显著程度。

我们把这些分布排列成矩阵，其中每行都是一个归一化的概率分布。

$$\left[ \begin{matrix} g^{(1)} \\ g^{(2)} \\ \vdots \\ g^{(n)} \end{matrix} \right] = \left[ \begin{matrix} g^{(1,1)} & g^{(1,2)} & \cdots & g^{(1,n)}\\ g^{(2,1)} & g^{(2,2)} & \cdots & g^{(2,n)}\\ \vdots & \vdots & \ddots & \vdots\\ g^{(n,1)} & g^{(n,2)} & \cdots & g^{(n,n)} \end{matrix} \right] %\tag{3}$$

下面我们需要将这n个独立的分转化为关于所有边的分布。这样一来 $g^{(i,j)}+g^{(j,i)}$ 代表点 $i$ 和 $j$ 的距离相应的概率。

首先我们进行归一化。由于有n个已经归一化的分布，所以为了使得总体分布归一化，我们把矩阵中每一项除以 $n$ 即可。

然后进行对称化。由于 $g^{(i,j)}$, $g^{(j,i)}$ 对应同一个距离，不应该有差异，所以我们让它们各取一半，即

$$p^{(i,j)} = p^{(j,i)} = \frac12\left(g^{(i,j)}+g^{(j,i)}\right)$$

最终我们把距离分布转化为了概率分布，即

$$\left[ \begin{matrix} d^{(1,1)} & d^{(1,2)} & \cdots & d^{(1,n)}\\ d^{(2,1)} & d^{(2,2)} & \cdots & d^{(2,n)}\\ \vdots & \vdots & \ddots & \vdots\\ d^{(n,1)} & d^{(n,2)} & \cdots & d^{(n,n)} \end{matrix} \right] \rightarrow \left[ \begin{matrix} p^{(1,1)} & p^{(1,2)} & \cdots & p^{(1,n)}\\ p^{(2,1)} & p^{(2,2)} & \cdots & p^{(2,n)}\\ \vdots & \vdots & \ddots & \vdots\\ p^{(n,1)} & p^{(n,2)} & \cdots & p^{(n,n)} \end{matrix} \right]$$

我们将最后得到的概率矩阵记为 $P_X$

SNE - Part IV 将低维度维空间中距离转化为概率分布

低维度中我们不应采用先前那样的带参数的概率密度函数，原因是如果使用就代表我们“先验”地定义结果的性质而不是输入的性质，这是一种“自欺欺人”的做法，即我们不应该为了更好的结果而改变结果的形式。

传统的 SNE 使用的是 $N(0, 1)$，即标准正态分布。但是事实证明使用自由度为1的（0均值）t 分布效果更好，对应的方法称为 t-SNE。

t 分布是统计学三大分布（ $\chi^2$ 分布，t 分布，F 分布）之一，又称为学生分布 (和一段署名为“学生”的匿名发表轶事有关）。自由度为1的密度函数为

$$t(X) = \frac1{\pi(1+x^2)}$$

它比标准正态分布更加平缓（红色为t分布）：

为什么 t 分布效果更好？一个可能的原因在于 t 分布本身的性质。正态分布用于估计方差为$\sigma$ 的数据的均值，而 t 分布用于估计方差未知的数据的均值，也就是t分布更加低维数据没有先验方差的情况，而标准正态分布假定了方差为1。

我们使用相同的方法将t分布应用于 $D_Y$ 并且归一化得到所有边的概率分布，记为 $Q_Y$。

SNE - Part V 优化

最终的E为

$$E_{t-SNE} = {\rm KL}(P_X\|Q_Y), with\ certain\ perplexity$$

优化目标为

$$\min_Y {\rm KL}(P_X\|Q_Y), with\ certain\ perplexity$$

t-SNE 的优化目标相对之前介绍的内容复杂很多，而且是非凸的（即存在大量的局部最优值），往往需要加上绝热退火、动量下降等方法。

下面是结果：

可以看出是相当不错的。

下面的图是 t-SNE 论文中的对比结果，t-SNE是碾压的胜利。

t-SNE 也有一些问题，比如速度慢，结果因为有局部最优所以有点随机，对于非常一般的一些数据处理不好等等。但是对于大型的真实数据集和对深度神经网络的分析上面 t-SNE 优势很大，所以常常用于深度学习的结果分析中。

t-SNE 的 perplexity 是必要的参数吗？

相信 t-SNE 让人不满意的地方之一在于有个 perplexity 参数。这个参数是必要的吗？

几乎可以说是必要的，而且很有意义，甚至是 t-SNE 优势之一。

我们考虑下面的情况：

我们把一条直线上面的点绕折构成一个“平面”：

然后把平面像寿司卷那样卷成寿司棒（下面的图切开了 T T）

侧视图

然后我们延长这个“棒子”，把它作为“直线”，继续如上的构造。

问题是，在二维空间它应该呈现出什么样子才能体现它的性质？一条线？一个面？一个圆饼？似乎都有道理，因为在于尺度的大小。上述构造的实质是不断把各个维度的“卷曲”作为信息加入到数据中，这些信息对应了多级结构，我们需要获取这些信息的时候就必须考虑到不同尺度的结构。

t-SNE 的 perplexity 恰恰控制的是尺度，原因如下：

由于采用了 gaussian 分布，所以我们知道大约68%的数据集中在一个标准差 $\sigma$ 内，也就是大约在一个 perplexity 的范围内。

perplexity 非常小的时候，在计算邻居距离的过程中，如下图所示，刚好能够覆盖最小尺度的“邻居”。

在全局的归一化和对称化完成后，实质上形成了这样的分布：

这样 t-SNE 会展现“直线”的姿态（一级结构）。

如果perplexity更大一点就会像下面一样

但是侧面覆盖范围还不够大

所以最终会得到展开的平面（二级结构）

再大就足以覆盖三级结构(棒状)

所以说 perplexity 一定程度上决定了我们看到的结构尺度。

Big Data talks

中心极限定理告诉我们，对于独立同分布的数据，数据样本越多则能越精确地估计数据的均值。统计学三大分布也告诉我们这个结论对估计方差比、均值差等也成立。

我们可以直觉上感到，随着数据增多，对于数据其他特性的估计也会变的更好，包括数据的“形状”。

之前采用的例子都是 1,000 个 MNIST 数字图像。我们试试用 2,000 个：

使用2000个图像的 Sammon mapping，改进有限。

使用2000个图像的 KNN-ForceField（k=3），进步挺大

使用2000个图像的 t-SNE (perplexity = 40)，进步更大:

t-SNE 处理的 MNIST 全体大合影（60,000个），效果非常好

可以看出 t-SNE 甚至保留了很多规律，比如右侧的1的位置和它的形态明显相关

另外图中也暗示了某些数字难以区分。

维度

之前我们都是把高维度投射到2维的空间中，如果是3维那么如何？ 下面我们先用 1000 张图片作测试。

Sammon Mapping，有提升

KNN-ForceField (k=3)

t-SNE 效果也不错

将数目增加到 2000个

KNN-ForceField (k=3) 侧视图

t-SNE 侧视图1

t-SNE 侧视图2

效果是不错的

一般情况下，降维到三维要优于到二维，这个也符合常理。但是三维需要“转着看”，所以论文和著作中更偏向二维的。

尾注

相关的代码参见 https://github.com/suquark/ConvNetJS6 中的 dim_reduction 模块 数据参见 dataset 模块（如果你觉得不错就 star 一个吧 ~~~）

可视化使用了 d3js 和 threejs，其中 d3js 用于控制颜色和选择控件，threejs 是个使用 WebGL 进行渲染的库。另外还有很多辅助的代码，参见 visualize 模块（轮子还没有完全造好）

这个是可以在浏览器里面跑的 demo（动态运行的），等项目整理好了我出demo。

另外参见了 http://colah.github.io/ 和 Andrej Karpathy 的一些blog和项目。很高兴我回过来看的时候发现很多想法是和他们想到一起了。