论文笔记：A tractable latent variable model for nonlinear dimensionality reduction

一、论文基本信息

• 发表期刊：《Proceedings of the National Academy of Sciences of the United States of America》
• 期刊最新影响因子:9.412
• 中科院分区：综合性期刊 1 区
• 论文作者及单位：加州大学圣地亚哥分校计算机科学与工程系劳伦斯·k·索尔（Lawrence K. Saul）
• 出版日期:2020年6月24日
• 原文链接：https://www.pnas.org/content/117/27/15403.short
• MATLAB代码链接： https://osf.io/wy793/
• MNIST书写数据集链接： http://yann.lecun.com/exdb/mnist/
• Matlab word2vec 数据集链接：https://github.com/chrisjmccormick/word2vec matlab

二、摘要

潜在变量模型（latent variable model）其他算法的共同点：

• 将高维的数据降低到二维或者三位实现可视化
• 拉近同一族之间的数据点，使得不同族数据之间的距离更加远

潜在变量模型（latent variable model）算法的创新点:

• 推导了一个期望最大化过程，该过程带有封闭形式更新，单调地提高模型的似然性:在这个过程中，通过求解由离散图Laplacian产生的稀疏的、对角占优的线性方程组，迭代地适应嵌入。
• 开发了一个近似的粗粒化过程，以避免对图中不相邻节点进行负抽样。在图像和文本数据集上验证了该模型的有效性。

三、常用非线性降维算法

t分布随机邻域嵌入(t-SNE) ， LargeVis，和一致流形逼近与投影(UMAP)产生了令人惊叹的高维数据集可视化。都采用相似性和差异性的概率概念。尽管有这些共同的基础，然而，这些方法与其他典型的非监督学习概率模型几乎没有相似之处。但是存在一个问题，只有附近的输入被映射到附近的输出，即这些算法在实现降维的过程中存在“拥挤问题”，本文提出的 LVM可以有效的解决这个问题。

四、LVM算法的理解

基本原理：
这些输入随机地映射到一个低维的潜在空间中，然后调整映射的参数，使得在高概率情况下，只有附近输入的图像仍然在附近。，模型的参数是通过最大化一个目标函数来估计的，该目标函数由两项组成:第一项度量附近输入保持在附近的可能性，而第二项度量远处输入保持在远处的可能性。在模型中，这些可能性是通过对低维潜在空间中(随机分配的)输入位置进行积分来计算的。
LVM算法的意义：
典型的LVMs包括因子分析(factor analysis)和高斯混合模型(Gaussian mixed model)，前者解释了大量观察到的变量与少量未观察到的变量之间的相关性，后者揭示了由潜在的多模态分布产生的数据簇。在本文中，我们描述了一个概念简单但同样有效的用于非线性降维(NLDR)的LVM，其目标是发现高维数据的忠实的、保持邻域的嵌入。NLDR工具可以帮助所有科学和工程领域的研究人员更好地理解和可视化他们的数据。我们的方法将NLDR提升到一个问题的家族，可以研究特别易处理LVMs。

五、 LVM算法的基本步骤

（一）构造输入和输出之间的多元高斯分布

本节中的符号及意义：

• $x_i$和$x_j$：表示一对特定的输入对
• $u_i$和$u_j$：表示对应的输出
• $P（h/x_i）$和$P（h/x_j）$：表示多元的高斯分布
• ${\sigma }_i^2$和${\sigma }_j^2$：表示方差
• $h$：表示一个潜在随机变量
• $h^{/}$：表示从多元高斯分布中的样本
• ${\Delta }_{ij}$：表示ij不相似的距离
• ${\delta }_{ij}$：表示ij相似的距离
• $c_{\Delta }$和$c_{\delta }$：表示伯努力随机变量
• $\zeta$：整体log(条件)可能性
• $S_{ij}$和$D_{ij}$非负权值附加在我们认为在一个长度尺度上相似和/或在另一个尺度上不同的输入对上

输入和输出之间的高斯分布 高斯分布如下：

$c_{\Delta }$和$c_{\delta }$为伯努力随机变量可以得到一个条件分布和推导过程如下：

可以通过在多个输入对上累积这些分数来衡量LVM嵌入的总体质量。正如我们将要看到的，在这个过程中投资一些比其他的作用更大的对投入将是方便的。考虑到分数的加权和，得到了一个整体的log(条件)可能性

（二）、期望最大化算法（EM Algorithm）

当应用到最大效果时，EM算法产生封闭形式的参数更新，单调地提高这些模型的可能性。EM算法被广泛应用于高斯混合模型、因子分析和隐马尔科夫模型的ML估计，其中后验统计量可以有效计算。
补充材料中给出了期望步长(e步长)和最大化步长(m步长)的完整推导,这篇文章提出的LVM算法主要集中以下在两个方面：

• 算法潜在变量是连续的，但后向分布并不总是多元高斯分布，精确的推断仍然是可控的。
• 与大多数用于无监督学习的LVMs不同，该算法不学习数据生成模型;相反，我们估计了LVM 的参数$(u_i,u_j{)}_{i=1}^n$来最大化Eq. 6中的log条件似然，正是这些参数提供了数据的低维表示.

（三）、E-step的推导

简单地用LVM连续潜变量的后验统计量来表示EM算法的参数更新。对于Eq. 6中出现的每对输入(xi, xj)，其后验分布由贝叶斯规则给出：

将Eqs. 1–5 代入到 Eq. 7中， Eq. 7是一个多元的高斯分布，其均值和方差与先验分布不同。Eq. 8不是多元高斯分布，但统计量仍然可以有效计算，定义了以下的优势比：

计算潜变量最重要的统计量是它们的后验均值。推导过程得到下列表达式：

SI附录给出了$E[h^{/}|c_{\sigma }=1,x_i,x_j]$和$E[h^/|c_\sigma=0,x_i,x_j]的对应表达式。注意，在每个表达式中，右边的第一项是先验平均值，E[h|xi] = i，而其余的一项给出了修正。这些修正的不同标志直观地表明，当xi和xj被建模为附近的输入时(c的值为1)，h和h0的后验平均值被拉在一起，而当它们不被建模时(c = 0)，它们被分开。EM算法需要的其他后验统计数据是潜在变量与其先验均值之间的欧式距离的平方期望值。
这些二阶统计量也可以有效地计算出来：

（四）、学习M-step

LVMs的EM算法基于对数似然值(15)的下界最大化。有以下两种方法为LVM派生EM更新如下：

• 第一种(标准)方法在Eq. 6上使用一个单一的下界来派生模型参数 { u i , σ i 2 } i = 1 2 \lbrace u_i,\sigma_i^2 \rbrace_{i=1}^2 {ui​,σi2​}i=12​的联合更新。
• 第二种方法使用单独的(和更紧密的)下界来派生这些参数的交替更新。这里我们展示了第二种方法的结果，它在实践中收敛速度明显更快。
  首先给出输出 { u i } i = 1 n \lbrace u_i \rbrace_{i=1}^n {ui​}i=1n​的更新。由相似度权重$s_{ij}$定义矩阵如下：
  
  eq.18是稀疏对称的。形成对称对角占优(SDD)矩阵
  
  矩阵L由两项的和组成:第一个是通过Eq. 18包含相似权值Sij的加权图Laplacian，第二个是非负的对角矩阵包含不同权重$D_{ij}$。通过求解线性系统，重新估计输出
  
  为了实现这个更新，我们可以利用SDD线性系统的极其高效的求解器,在实践中，更新利用一些二阶结构的目标函数，同时避免了复杂的牛顿拉夫森更新。它还例证了所谓的Laplacian范式(18)，在这种范式中，SDD系统的快速求解器被用作涉及大型图的问题的算法.
  对方差进行更新，采用了一种特别简单的形式，这种更新和Eq. 20中的更新都保证增加Eq. 6中的对数条件似然性(除了在静止点)。:
  

（五）、长度尺度和成对相似性（Length Scales and Pairwise Similarities）

LVM的目标是学习将附近输入映射到附近输出和将远处输入映射到远处输出的嵌入方法。
然而，必须将我们认为邻近的输入对形式化。通常，我们用稀疏图来编码这些邻域关系。然后利用此图推导出Eq. 3、4中的长度比例尺(等式${\delta }_{ij},{\Delta }_{ij}$)和Eq. 6中的权重($S_{ij},D_{ij}$)

（六）、邻域图（neighborhood graph）

通过KNN生成了一个稀疏对称图，其边连接邻近的输入，但在实践中，它经常连接不应该被视为邻近的输入。当数据包含异常值时，当输入没有被密集采样时，或者当潜在密度是可变的时，问题就会出现(19)。许多流形学习的失败模式可以追溯到这样形成的图，KNN算法如下图所示：
图一、KNN基本原理
对于我们的LVMs，我们使用了一个简单的双参数过程来构造稀疏邻域图，其中两个参数(k, s)都是小的正整数。参数k在很大程度上决定了邻域大小，为了可视化，可以将其设置为缺省值s = 1。
计算每个输入的kNN开始，这些邻域关系用二进制矩阵K编码：

$K_{ij}$这个矩阵是不对称的，在实践中，我们发现，即使k值很小，这个图作为应该映射到附近输出的输入对的表示也可能过于宽松，所以我们程序的其余部分被设计为在这个图中提取一个表示更健壮的邻域关系的边子集(19)。其过程主要分为一下三部：

1. 通过识别产生最小连通子图的输入对，计算邻接矩阵$K_{ij}+K_{ji}-K_{ij}{K}_{ji}$i的无向图的最小生成树（如果其中一个输入是另一个输入的kNN，则其边连接$x_i到x_j$，如果图是不连通的，我们只考虑它最大的连通部分），假设(对于一个合适的k值)以包含除少数异常值之外的所有异常值。我们用T表示这个最小加权生成树的邻接矩阵。矩阵T编码了保持输入的大规模连通性的最小边集。
2. 通过确定无疑值得视为附近的输入对
3. 结合这些发现。现在我们将更详细地描述每一个步骤。

在第二步中，我们累加矩阵K的前s次幂;s的值(通常很小)可以看作是一个有向随机游走的步数。

s的值越大，s步随机游走的范围越广。我们特别感兴趣的是输入对(xi, xj)，它们通过随机游走相互可达:这些对$R_{ij}{R}_{ji}$ > 0，它们往往不包括异常值。我们已经观察到，在这个有向随机游走中的互可达性比简单的kNN揭示了更多有意义的邻域.
我们的过程的最后一步结合了前两个结果。我们定义一个有边的邻域图:

通过构造，这个邻域图只连接了这些输入对的一个子集，这些输入是kNN;特别是，它保留了那些需要保持图表连接的对，以及那些相互关系表明了额外程度的稳健性的对。注意，这些邻域关系通常不是对称的，因此我们有时可能认为$x_j$在$x_i$附近(当Eij = 1)，但反之亦然(当Eji = 0时)。

(七)、粗粒化（coarse graining）

下图中出现的边表示应该映射到附近输出的输入对，而不存在的边表示不应该映射到附近输出的输入对，图中出现的边表示应该映射到附近输出的输入对，而不存在的边表示不应该映射到附近输出的输入对，当邻域图是稀疏的时，后一类中的对数将大到足以枚举除最小数据集之外的所有数据集。递归的级别数$(\iota )$可以调整来适应大量的输入中的学习。

图二、粗粒化程序的效果，用红色边缘连接显式建模为不相似的输入(例如手写数字的图像)。(左)粗粒化之前，许多对这样的输入将大量的项贡献到模型s log条件似然中。(右)粗粒化后，只有指定为界标的输入被明确地建模为不相似的，但它们对对数条件似然的贡献由它们所代表的输入的数量加权;此外，由虚线方格内的相同地标所表示的输入被建模为与地标(由新邻域图中的蓝色边表示)微弱相似。

1、递归级数 L=0
选择权重和长度尺度的目的主要有一下两点：

• 缩小邻居之间的距离
• 确保邻居比非邻居更亲密（解决拥挤问题）
  在这个过程中设置权重$S_{ij}=E_{ij}$和$D_{ij}=1-E_{ij}$（对角线上保留非零元素），在进行EM算法之前对权重进行处理，使得${\Sigma }_{ij}{S}_{ij}={\Sigma }_{i}j{D}_{ij}$

2、递归级数 L=1
图1显示了粗粒化背后的基本直觉。假设n太大，无法对$O(n^2)$对不同输入进行建模。相反，我们随机指定$n_l<n$个输入作为地标，并将每个输入赋给它最近的地标;类似的思想(20,23)已经被用于许多NLDR算法中。让我表示分配给landmark $x_i$的输入集合(包括它自身)，让${\Omega }_i$表示这个集合的大小。粗粒化并不改变分配给同一个landmark的输入之间的不同权重。做了三个改变:

• 在地标之间与其他输入之间添加单元相似度权重
• 用不同地标表示的输入之间的不同权重
• 在不同的地标(比如xi和xj)之间添加不同的权重
  
  如果地标被广泛地分开，而它们所代表的输入相对地聚集在一起，那么，通过三角形不等式，不同地标所代表的输入也将被很好地分开。Eq. 26中的长度标度就是基于这一点来选择的，这样，在粗粒化之后的ML估计将把属于不同地标的输入至少与粗粒化之前的量分开。
  2、递归级数 L=>2
  对于较大的n，它可能仍然禁止显式地建模地标对之间的不同和/或分配给同一地标的对输入之间的不同。在这种情况下，可以递归地应用粗粒化过程。

六、结果

图三、7万手写数字的2D可视化(k = 9, s = 1，L = 1)

图4、图三的热图

图5、UMAP和T-SNE的可视化效果

图6、PCA与潜变量建模降维后手写数字9NN分类器的测试错误率(k = 9, s = 2)，当L= 2时，大概是因为额外的粗粒化引起了三角形不等式的松弛。作为比较，对于具有L = 1和L= 2级粗粒化的d = 10嵌入件，400次EM迭代分别需要859和579 s
还对在GoogleNews语料库上训练的word2vec模型的标准化词向量子集进行了实验(7)。原始模型为300万单词的词汇表中的每个单词生成维数为300的向量;该子集包含n = 207,147个单词向量，通过过滤出WordNet(25)中没有列出的单词(如拼写错误和垃圾短语)获得。

图7、对于n = 207,147个归一化词向量的数据集，利用潜在变量建模，通过PCA和NLDR进行线性降维(k = 4, s = 3， L= 2)，得到r处的第一近的召回率。

表1、表1显示了我们的LVM (k = 4, s = 3， L= 2)将单词向量嵌入10维之前和之后的4个最近邻。该嵌入的计算时间主要由初始搜索精确的最近近邻所占时间超过2小时;剩下的计算大约花了30分钟.
实验补充图：

图8、这张图显示了与图2相同的嵌入，但排除了最不确定的输出——那些方差超过95%分位数的输出。结果图的散斑更少，表明非密集簇的输出对应的输出具有更大的不确定性。

图9、上面的图显示了与图2相同的嵌入，但是圈出了方差最大的100个输出。下图显示了与圈出的输出对应的数字图像。

图10、有相同的方差:结果嵌入不分离最相似的数字类;它也有一个相当低的可能性

图11、不同k和s值的MNIST位数的嵌入。也显示了MNIST测试集在低维潜在空间中的1-最近邻误差

图12、左:特写的顶部集群3、5、8。右:特写底部的一组4、7、9。与下一幅图进行比较，看看进一步粗粒化的效果。

图13、左:特写的底部集群2、3、5、8。右:特写顶部的4和9的集群。这个图说明了粗粒化的一个副作用:集群之间的距离被放大了

图14、通过不同级别的粗粒化(’ = 0 vs ’ = 1) LVMs (k = 9, s= 1)对n= 7000 MNIST手写数字进行二维可视化。在前面的图中，我们观察到粗粒化倾向于放大集群之间的距离

图15、使用不同级别的粗粒化(’ = 0 vs ’ = 1) LVMs (k = 9, s= 1)对瑞士滚动数据集(n = 5000)进行二维可视化

图16、完成400次EM迭代的时间(秒)与LVMs嵌入的维数(见图3)。

七、未来的发展方向

• NLDR的一般问题有比可视化更大的目标，换句话说就是可视化大于3维
• 可以加入额外的潜在变量——而可控制的LVMs可以作为不太容易控制的LVMs的近似推断的基础