系统学习机器学习之维度归约(完整篇)

这里,我们讨论特征选择和特征提取,前者选取重要的特征子集,后者由原始输入形成较少的新特征,理想情况下,无论是分类还是回归,我们不应该将特征选择或特征提取作为一个单独的进程,分类或者回归方法应该能够利用任何必要的特征,而丢弃不相关的特征。但是,考虑到算法存储量和时间的复杂度,或者输入不必要的特征等原因,还是需要降维。较简单的模型在小数据上更为鲁棒,有小方差,模型的变化更依赖于样本的特殊性,包括噪声,离群点等。同时,低纬度描述数据,方便我们队数据绘图,可视化分析数据结构和离群点。

降低方法,一般为特征选择,特征提取,按监督和非监督分:
非监督:PCA、FA ,LLE
监督:LDA,MDS
其中,PCA,FA,MDS,LDA,都是线性投影方法,非线性维度规约有等距特征映射,局部线性嵌入LLE。

维度归约使用数据编码或变换,以便得到原数据的归约或“压缩”表示。如果原数据可以由压缩数据重新构造而不丢失任何信息,则该数据归约是无损的。如果我们只能重新构造原数据的近似表示,则该数据归约是有损的。有一些很好的串压缩算法。尽管它们通常是无损的,但是只允许有限的数据操作。

一.先说特征选择,即子集选择。

在子集选择中,我们选择最佳子集,其含的维度最少,但对正确率的贡献最大。在维度较大时,采用启发式方法,在合理的时间内得到一个合理解(但不是最优解)。维度较小时,对所有子集做检验。
有两种方法,向前选择,即从空集开始逐渐增加特征,每次添加一个降低误差最多的变量,直到进一步添加不会降低误差或者降低很少。同时,可以用浮动搜索,每一步可以改变增加和去掉的特征数量,以此来加速。
向后选择中,从所有变量开始,逐个排除他们,每次排除一个降低误差最多的变量,直到进一步的排除会显著提高误差。如果我们预料有许多无用特征时,向前选择更可取。
在两种情况下,误差检测都应在不同于训练集的验证集上做,因为我们想要检验泛化准确率。使用更多的特征,我们一般会有更低的训练误差,但不一定有更低的验证误差。
像人脸识别这样的应用中,特征选择不是很好的降维方法,因为个体像素本身并不携带很多识别信息,携带脸部识别信息的是许多像素值的组合。这可以用特征提取来归约。

二. 主成分分析(PCA)

1. 问题

     真实的训练数据总是存在各种各样的问题:

1、 比如拿到一个汽车的样本,里面既有以“千米/每小时”度量的最大速度特征,也有“英里/小时”的最大速度特征,显然这两个特征有一个多余。

2、 拿到一个数学系的本科生期末考试成绩单,里面有三列,一列是对数学的兴趣程度,一列是复习时间,还有一列是考试成绩。我们知道要学好数学,需要有浓厚的兴趣,所以第二项与第一项强相关,第三项和第二项也是强相关。那是不是可以合并第一项和第二项呢?

3、 拿到一个样本,特征非常多,而样例特别少,这样用回归去直接拟合非常困难,容易过度拟合。比如北京的房价:假设房子的特征是(大小、位置、朝向、是否学区房、建造年代、是否二手、层数、所在层数),搞了这么多特征,结果只有不到十个房子的样例。要拟合房子特征->房价的这么多特征,就会造成过度拟合。

4、 这个与第二个有点类似,假设在IR中我们建立的文档-词项矩阵中,有两个词项为“learn”和“study”,在传统的向量空间模型中,认为两者独立。然而从语义的角度来讲,两者是相似的,而且两者出现频率也类似,是不是可以合成为一个特征呢?

5、 在信号传输过程中,由于信道不是理想的,信道另一端收到的信号会有噪音扰动,那么怎么滤去这些噪音呢?

     回顾我们之前介绍的《模型选择和规则化》,里面谈到的特征选择的问题。但在那篇中要剔除的特征主要是和类标签无关的特征。比如“学生的名字”就和他的“成绩”无关,使用的是互信息的方法。

     而这里的特征很多是和类标签有关的,但里面存在噪声或者冗余。在这种情况下,需要一种特征降维的方法来减少特征数,减少噪音和冗余,减少过度拟合的可能性。

     下面探讨一种称作主成分分析(PCA)的方法来解决部分上述问题。PCA的思想是将n维特征映射到k维上(k<n),这k维是全新的正交特征。这k维特征称为主元,是重新构造出来的k维特征,而不是简单地从n维特征中去除其余n-k维特征。

2. PCA计算过程

     首先介绍PCA的计算过程:

     假设我们得到的2维数据如下:

     clip_image001[4]

     行代表了样例,列代表特征,这里有10个样例,每个样例两个特征。可以这样认为,有10篇文档,x是10篇文档中“learn”出现的TF-IDF,y是10篇文档中“study”出现的TF-IDF。也可以认为有10辆汽车,x是千米/小时的速度,y是英里/小时的速度,等等。

     第一步分别求x和y的平均值,然后对于所有的样例,都减去对应的均值。这里x的均值是1.81,y的均值是1.91,那么一个样例减去均值后即为(0.69,0.49),得到

     clip_image002[4]

     第二步,求特征协方差矩阵,如果数据是3维,那么协方差矩阵是

     clip_image003[4]

     这里只有x和y,求解得

     clip_image004[4]

     对角线上分别是x和y的方差,非对角线上是协方差。协方差大于0表示x和y若有一个增,另一个也增;小于0表示一个增,一个减;协方差为0时,两者独立。协方差绝对值越大,两者对彼此的影响越大,反之越小。

     第三步,求协方差的特征值和特征向量,得到

     clip_image005[4]

     上面是两个特征值,下面是对应的特征向量,特征值0.0490833989对应特征向量为clip_image007[4],这里的特征向量都归一化为单位向量。

    第四步,将特征值按照从大到小的顺序排序,选择其中最大的k个,然后将其对应的k个特征向量分别作为列向量组成特征向量矩阵。

     这里特征值只有两个,我们选择其中最大的那个,这里是1.28402771,对应的特征向量是clip_image009[6]

     第五步,将样本点投影到选取的特征向量上。假设样例数为m,特征数为n,减去均值后的样本矩阵为DataAdjust(m*n),协方差矩阵是n*n,选取的k个特征向量组成的矩阵为EigenVectors(n*k)。那么投影后的数据FinalData为

     clip_image011[4]

     这里是

     FinalData(10*1) = DataAdjust(10*2矩阵)×特征向量clip_image009[7]

     得到结果是

     clip_image012[4]

     这样,就将原始样例的n维特征变成了k维,这k维就是原始特征在k维上的投影。

     上面的数据可以认为是learn和study特征融合为一个新的特征叫做LS特征,该特征基本上代表了这两个特征。

     上述过程有个图描述:

     clip_image013[4]

     正号表示预处理后的样本点,斜着的两条线就分别是正交的特征向量(由于协方差矩阵是对称的,因此其特征向量正交),最后一步的矩阵乘法就是将原始样本点分别往特征向量对应的轴上做投影。

     如果取的k=2,那么结果是

     clip_image014[4]

     这就是经过PCA处理后的样本数据,水平轴(上面举例为LS特征)基本上可以代表全部样本点。整个过程看起来就像将坐标系做了旋转,当然二维可以图形化表示,高维就不行了。上面的如果k=1,那么只会留下这里的水平轴,轴上是所有点在该轴的投影。

     这样PCA的过程基本结束。在第一步减均值之后,其实应该还有一步对特征做方差归一化。比如一个特征是汽车速度(0到100),一个是汽车的座位数(2到6),显然第二个的方差比第一个小。因此,如果样本特征中存在这种情况,那么在第一步之后,求每个特征的标准差clip_image016[6],然后对每个样例在该特征下的数据除以clip_image016[7]

     归纳一下,使用我们之前熟悉的表示方法,在求协方差之前的步骤是:

     clip_image017[4]

     其中clip_image019[6]是样例,共m个,每个样例n个特征,也就是说clip_image019[7]是n维向量。clip_image021[4]是第i个样例的第j个特征。clip_image023[4]是样例均值。clip_image025[4]是第j个特征的标准差。

     整个PCA过程貌似及其简单,就是求协方差的特征值和特征向量,然后做数据转换。但是有没有觉得很神奇,为什么求协方差的特征向量就是最理想的k维向量?其背后隐藏的意义是什么?整个PCA的意义是什么?

3. PCA理论基础

     要解释为什么协方差矩阵的特征向量就是k维理想特征,我看到的有三个理论:分别是最大方差理论、最小错误理论和坐标轴相关度理论。这里简单探讨前两种,最后一种在讨论PCA意义时简单概述。

3.1 最大方差理论

     在信号处理中认为信号具有较大的方差,噪声有较小的方差,信噪比就是信号与噪声的方差比,越大越好。如前面的图,样本在横轴上的投影方差较大,在纵轴上的投影方差较小,那么认为纵轴上的投影是由噪声引起的。

因此我们认为,最好的k维特征是将n维样本点转换为k维后,每一维上的样本方差都很大。

     比如下图有5个样本点:(已经做过预处理,均值为0,特征方差归一)

     clip_image026[4]

     下面将样本投影到某一维上,这里用一条过原点的直线表示(前处理的过程实质是将原点移到样本点的中心点)。

     clip_image028[4]

     假设我们选择两条不同的直线做投影,那么左右两条中哪个好呢?根据我们之前的方差最大化理论,左边的好,因为投影后的样本点之间方差最大。

     这里先解释一下投影的概念:

     QQ截图未命名

     红色点表示样例clip_image037[14],蓝色点表示clip_image037[15]在u上的投影,u是直线的斜率也是直线的方向向量,而且是单位向量。蓝色点是clip_image037[16]在u上的投影点,离原点的距离是clip_image039[4](即clip_image030[4]或者clip_image041[4])由于这些样本点(样例)的每一维特征均值都为0,因此投影到u上的样本点(只有一个到原点的距离值)的均值仍然是0。

     回到上面左右图中的左图,我们要求的是最佳的u,使得投影后的样本点方差最大。

     由于投影后均值为0,因此方差为:

     clip_image042[4]

     中间那部分很熟悉啊,不就是样本特征的协方差矩阵么(clip_image037[17]的均值为0,一般协方差矩阵都除以m-1,这里用m)。

     用clip_image044[10]来表示clip_image046[4]clip_image048[6]表示clip_image050[4],那么上式写作

     clip_image052[4] 

     由于u是单位向量,即clip_image054[4],上式两边都左乘u得,clip_image056[4]

     即clip_image058[4]

     We got it!clip_image044[11]就是clip_image048[7]的特征值,u是特征向量。最佳的投影直线是特征值clip_image044[12]最大时对应的特征向量,其次是clip_image044[13]第二大对应的特征向量,依次类推。

     因此,我们只需要对协方差矩阵进行特征值分解,得到的前k大特征值对应的特征向量就是最佳的k维新特征,而且这k维新特征是正交的。得到前k个u以后,样例clip_image037[18]通过以下变换可以得到新的样本。

     clip_image059[4]

     其中的第j维就是clip_image037[19]clip_image061[4]上的投影。

     通过选取最大的k个u,使得方差较小的特征(如噪声)被丢弃。

     这是其中一种对PCA的解释

3.2 最小平方误差理论

     clip_image001

     假设有这样的二维样本点(红色点),回顾我们前面探讨的是求一条直线,使得样本点投影到直线上的点的方差最大。本质是求直线,那么度量直线求的好不好,不仅仅只有方差最大化的方法。再回想我们最开始学习的线性回归等,目的也是求一个线性函数使得直线能够最佳拟合样本点,那么我们能不能认为最佳的直线就是回归后的直线呢?回归时我们的最小二乘法度量的是样本点到直线的坐标轴距离。比如这个问题中,特征是x,类标签是y。回归时最小二乘法度量的是距离d。如果使用回归方法来度量最佳直线,那么就是直接在原始样本上做回归了,跟特征选择就没什么关系了。

     因此,我们打算选用另外一种评价直线好坏的方法,使用点到直线的距离d’来度量。

     现在有n个样本点clip_image003,每个样本点为m维(这节内容中使用的符号与上面的不太一致,需要重新理解符号的意义)。将样本点clip_image005在直线上的投影记为clip_image007,那么我们就是要最小化

     clip_image009

     这个公式称作最小平方误差(Least Squared Error)。

     而确定一条直线,一般只需要确定一个点,并且确定方向即可。

     第一步确定点:

     假设要在空间中找一点clip_image011来代表这n个样本点,“代表”这个词不是量化的,因此要量化的话,我们就是要找一个m维的点clip_image011[1],使得

     clip_image012

     最小。其中clip_image014是平方错误评价函数(squared-error criterion function),假设m为n个样本点的均值:

     clip_image015

     那么平方错误可以写作:

     clip_image017

     后项与clip_image019无关,看做常量,而clip_image021,因此最小化clip_image014[1]时,

     clip_image023 

     clip_image019[1]是样本点均值。

     第二步确定方向:

     我们从clip_image019[2]拉出要求的直线(这条直线要过点m),假设直线的方向是单位向量e。那么直线上任意一点,比如clip_image007[1]就可以用点m和e来表示

     clip_image025 

     其中clip_image027clip_image029到点m的距离。

     我们重新定义最小平方误差:

     clip_image030

     这里的k只是相当于i。clip_image032就是最小平方误差函数,其中的未知参数是clip_image034和e。

     实际上是求clip_image032[1]的最小值。首先将上式展开:

     clip_image036

     我们首先固定e,将其看做是常量,clip_image038,然后对clip_image027[1]进行求导,得

     clip_image039

     这个结果意思是说,如果知道了e,那么将clip_image041与e做内积,就可以知道了clip_image043在e上的投影离m的长度距离,不过这个结果不用求都知道。

     然后是固定clip_image027[2],对e求偏导数,我们先将公式(8)代入clip_image032[2],得 

     clip_image044

     其中clip_image045 与协方差矩阵类似,只是缺少个分母n-1,我们称之为散列矩阵(scatter matrix)。

     然后可以对e求偏导数,但是e需要首先满足clip_image038[1],引入拉格朗日乘子clip_image047,来使clip_image049最大(clip_image032[3]最小),令

     clip_image050

     求偏导

     clip_image051

     这里存在对向量求导数的技巧,方法这里不多做介绍。可以去看一些关于矩阵微积分的资料,这里求导时可以将clip_image049[1]看作是clip_image053,将clip_image055看做是clip_image057

     导数等于0时,得

     clip_image058

     两边除以n-1就变成了,对协方差矩阵求特征值向量了。

     从不同的思路出发,最后得到同一个结果,对协方差矩阵求特征向量,求得后特征向量上就成为了新的坐标,如下图:

     clip_image059

     这时候点都聚集在新的坐标轴周围,因为我们使用的最小平方误差的意义就在此。

4. PCA理论意义

     PCA将n个特征降维到k个,可以用来进行数据压缩,如果100维的向量最后可以用10维来表示,那么压缩率为90%。同样图像处理领域的KL变换使用PCA做图像压缩。但PCA要保证降维后,还要保证数据的特性损失最小。再看回顾一下PCA的效果。经过PCA处理后,二维数据投影到一维上可以有以下几种情况:

     clip_image060

     我们认为左图好,一方面是投影后方差最大,一方面是点到直线的距离平方和最小,而且直线过样本点的中心点。为什么右边的投影效果比较差?直觉是因为坐标轴之间相关,以至于去掉一个坐标轴,就会使得坐标点无法被单独一个坐标轴确定。

     PCA得到的k个坐标轴实际上是k个特征向量,由于协方差矩阵对称,因此k个特征向量正交。看下面的计算过程。

     假设我们还是用clip_image062来表示样例,m个样例,n个特征。特征向量为e,clip_image064表示第i个特征向量的第1维。那么原始样本特征方程可以用下面式子来表示:

     前面两个矩阵乘积就是协方差矩阵clip_image066(除以m后),原始的样本矩阵A是第二个矩阵m*n。

     clip_image068

     上式可以简写为clip_image070

     我们最后得到的投影结果是clip_image072,E是k个特征向量组成的矩阵,展开如下:

     clip_image074

     得到的新的样例矩阵就是m个样例到k个特征向量的投影,也是这k个特征向量的线性组合。e之间是正交的。从矩阵乘法中可以看出,PCA所做的变换是将原始样本点(n维),投影到k个正交的坐标系中去,丢弃其他维度的信息。举个例子,假设宇宙是n维的(霍金说是11维的),我们得到银河系中每个星星的坐标(相对于银河系中心的n维向量),然而我们想用二维坐标去逼近这些样本点,假设算出来的协方差矩阵的特征向量分别是图中的水平和竖直方向,那么我们建议以银河系中心为原点的x和y坐标轴,所有的星星都投影到x和y上,得到下面的图片。然而我们丢弃了每个星星离我们的远近距离等信息。

     clip_image075

5. 总结与讨论

     这一部分来自http://www.cad.zju.edu.cn/home/chenlu/pca.htm

     PCA技术的一大好处是对数据进行降维的处理。我们可以对新求出的“主元”向量的重要性进行排序,根据需要取前面最重要的部分,将后面的维数省去,可以达到降维从而简化模型或是对数据进行压缩的效果。同时最大程度的保持了原有数据的信息。

     PCA技术的一个很大的优点是,它是完全无参数限制的。在PCA的计算过程中完全不需要人为的设定参数或是根据任何经验模型对计算进行干预,最后的结果只与数据相关,与用户是独立的。
但是,这一点同时也可以看作是缺点。如果用户对观测对象有一定的先验知识,掌握了数据的一些特征,却无法通过参数化等方法对处理过程进行干预,可能会得不到预期的效果,效率也不高。

     clip_image076

     图表 4:黑色点表示采样数据,排列成转盘的形状。
     容易想象,该数据的主元是clip_image077或是旋转角clip_image078

     如图表 4中的例子,PCA找出的主元将是clip_image077[1]。但是这显然不是最优和最简化的主元。clip_image077[2]之间存在着非线性的关系。根据先验的知识可知旋转角clip_image078[1]是最优的主元(类比极坐标)。则在这种情况下,PCA就会失效。但是,如果加入先验的知识,对数据进行某种划归,就可以将数据转化为以clip_image078[2]为线性的空间中。这类根据先验知识对数据预先进行非线性转换的方法就成为kernel-PCA,它扩展了PCA能够处理的问题的范围,又可以结合一些先验约束,是比较流行的方法。

     有时数据的分布并不是满足高斯分布。如图表 5所示,在非高斯分布的情况下,PCA方法得出的主元可能并不是最优的。在寻找主元时不能将方差作为衡量重要性的标准。要根据数据的分布情况选择合适的描述完全分布的变量,然后根据概率分布式

     clip_image079

     来计算两个向量上数据分布的相关性。等价的,保持主元间的正交假设,寻找的主元同样要使clip_image080。这一类方法被称为独立主元分解(ICA)。

     clip_image081

     图表 5:数据的分布并不满足高斯分布,呈明显的十字星状。
     这种情况下,方差最大的方向并不是最优主元方向。

     另外PCA还可以用于预测矩阵中缺失的元素。

6. 其他参考文献

     A tutorial on Principal Components Analysis LI Smith – 2002

     A Tutorial on Principal Component Analysis J Shlens

     http://www.cmlab.csie.ntu.edu.tw/~cyy/learning/tutorials/PCAMissingData.pdf

     http://www.cad.zju.edu.cn/home/chenlu/pca.htm

所谓主成分,就是在某方向上样本投影之后,被广泛散布,使得样本之间的差别变得最明显。为了得到唯一解并且使得该方向成为最重要因素,要求该向量长度为1,利用拉格朗日,计算出在方差最大情况下,我们选择具有最大特征值的特征向量为该投影向量。后面次方差最大。。。如此类推。由于协方差矩阵是对称的,因此,对于两个不同的特征值,特征向量是正交的。
通常情况下,我们考虑贡献90%以上方差的前K个分量,或者利用斜坡图,目视化分析拐点k,也可以忽略特征值小于平均输入方差的特征向量。
当然,如果原维之间不相关,则PCA就没有收益。在许多图像和语音处理任务中,邻近的输入时高度相关的。
如果原维的方差变化显著,则他们对主成分方向的影响比相关性大。因此,一般在PCA之前预处理下,使得每个维都具有0均值和单位方差。或者,为了使协方差而不是个体方差起作用,我们用协相关矩阵R而不是协方差矩阵S的本征向量。
离群点对方差有很大影响,从而影响特征向量,鲁棒的估计方法是允许计算离群点存在时的参数,如计算数据的马氏距离,丢弃哪些远离的孤立点数据。
当原维度很大时,我们可以直接从数据计算特征向量,特征值,不必计算协方差矩阵。
在所有正交线性投影中,PCA最小化重构误差。
这里,作为一种维度归约方法,我们直观地介绍主成分分析。
假定待归约的数据由n个属性或维描述的元组或数据向量组成。主成分分析(principal components analysis)或PCA(又称Karhunen-Loeve或K-L方法)搜索k个最能代表数据的n维正交向量,其中k≤n。这样,原来的数据投影到一个小得多的空间,导致维度归约。不像属性子集选择通过保留原属性集的一个子集来减少属性集的大小,PCA通过创建一个替换的、更小的变量集“组合”属性的基本要素。原数据可以投影到该较小的集合中。PCA常常揭示先前未曾察觉的联系,并因此允许解释不寻常的结果。

基本过程如下:

(1)对输入数据规范化,使得每个属性都落入相同的区间。此步有助于确保具有较大定义域的属性不会支配具有较小定义域的属性。
(2)PCA计算k个标准正交向量,作为规范化输入数据的基。这些是单位向量,每一个方向都垂直于另一个。这些向量称为主成分。输入数据是主成分的线性组合。

(3)对主成分按“重要性”或强度降序排列。主成分基本上充当数据的新坐标轴,提供关于方差的重要信息。也就是说,对坐标轴进行排序,使得第一个坐标轴显示数据的最大方差,第二个显示次大方差,如此下去。例如,图2-17显示原来映射到轴X1和X2的给定数据集的前两个主成分Y1和Y2。这一信息帮助识别数据中的分组或模式。

(4)既然主成分根据“重要性”降序排列,就可以通过去掉较弱的成分(即方差较小)来归约数据的规模。使用最强的主成分,应当能够重构原数据的很好的近似。

图2-17 主成分分析。Y1和Y2是给定数据的前两个主成分

PCA计算开销低,可以用于有序和无序的属性,并且可以处理稀疏和倾斜数据。多于2维的多维数据可以通过将问题归约为2维问题来处理。主成分可以用作多元回归和聚类分析的输入。与小波变换相比,PCA能够更好地处理稀疏数据,而小波变换更适合高维数据。


三.FA


FA也是非监督的,但与PCA相反,其旨在找到较少数量的因子,刻画观测变量之间的依赖性。但是,对于维度归约,除了因子可解释性,允许识别公共原因,简单解释,知识提取外,FA与PCA相比,并无优势。

1 问题

     之前我们考虑的训练数据中样例clip_image002的个数m都远远大于其特征个数n,这样不管是进行回归、聚类等都没有太大的问题。然而当训练样例个数m太小,甚至m<<n的时候,使用梯度下降法进行回归时,如果初值不同,得到的参数结果会有很大偏差(因为方程数小于参数个数)。另外,如果使用多元高斯分布(Multivariate Gaussian distribution)对数据进行拟合时,也会有问题。让我们来演算一下,看看会有什么问题:

多元高斯分布的参数估计公式如下:

     clip_image004

     clip_image006

     分别是求mean和协方差的公式,clip_image002[1]表示样例,共有m个,每个样例n个特征,因此clip_image008是n维向量,clip_image010是n*n协方差矩阵。

     当m<<n时,我们会发现clip_image010[1]是奇异阵(clip_image012),也就是说clip_image014不存在,没办法拟合出多元高斯分布了,确切的说是我们估计不出来clip_image010[2]

     如果我们仍然想用多元高斯分布来估计样本,那怎么办呢?

2 限制协方差矩阵

     当没有足够的数据去估计clip_image010[3]时,那么只能对模型参数进行一定假设,之前我们想估计出完全的clip_image010[4](矩阵中的全部元素),现在我们假设clip_image010[5]就是对角阵(各特征间相互独立),那么我们只需要计算每个特征的方差即可,最后的clip_image010[6]只有对角线上的元素不为0

     clip_image016

     回想我们之前讨论过的二维多元高斯分布的几何特性,在平面上的投影是个椭圆,中心点由clip_image008[1]决定,椭圆的形状由clip_image010[7]决定。clip_image010[8]如果变成对角阵,就意味着椭圆的两个轴都和坐标轴平行了。

     clip_image018

     如果我们想对clip_image010[9]进一步限制的话,可以假设对角线上的元素都是等值的。

     clip_image020

     其中

     clip_image022

     也就是上一步对角线上元素的均值,反映到二维高斯分布图上就是椭圆变成圆。

     当我们要估计出完整的clip_image010[10]时,我们需要m>=n+1才能保证在最大似然估计下得出的clip_image010[11]是非奇异的。然而在上面的任何一种假设限定条件下,只要m>=2都可以估计出限定的clip_image010[12]

     这样做的缺点也是显然易见的,我们认为特征间独立,这个假设太强。接下来,我们给出一种称为因子分析的方法,使用更多的参数来分析特征间的关系,并且不需要计算一个完整的clip_image010[13]

3 边缘和条件高斯分布

     在讨论因子分析之前,先看看多元高斯分布中,条件和边缘高斯分布的求法。这个在后面因子分析的EM推导中有用。

     假设x是有两个随机向量组成(可以看作是将之前的clip_image024分成了两部分)

     clip_image026

     其中clip_image028clip_image030,那么clip_image032。假设x服从多元高斯分布clip_image034,其中

     clip_image035

     其中clip_image037clip_image039,那么clip_image041clip_image043,由于协方差矩阵是对称阵,因此clip_image045

     整体看来clip_image047clip_image049联合分布符合多元高斯分布。

     那么只知道联合分布的情况下,如何求得clip_image047[1]的边缘分布呢?从上面的clip_image008[2]clip_image010[14]可以看出,

     clip_image051,clip_image053,下面我们验证第二个结果

     clip_image054

     由此可见,多元高斯分布的边缘分布仍然是多元高斯分布。也就是说clip_image056

     上面Cov(x)里面有趣的是clip_image058,这个与之前计算协方差的效果不同。之前的协方差矩阵都是针对一个随机变量(多维向量)来说的,而clip_image058[1]评价的是两个随机向量之间的关系。比如clip_image060={身高,体重},clip_image049[1]={性别,收入},那么clip_image062求的是身高与身高,身高与体重,体重与体重的协方差。而clip_image058[2]求的是身高与性别,身高与收入,体重与性别,体重与收入的协方差,看起来与之前的大不一样,比较诡异的求法。

     上面求的是边缘分布,让我们考虑一下条件分布的问题,也就是clip_image064的问题。根据多元高斯分布的定义,clip_image066

     且

     clip_image067

     这是我们接下来计算时需要的公式,这两个公式直接给出,没有推导过程。如果想了解具体的推导过程,可以参见Chuong B. Do写的《Gaussian processes》。

4 因子分析例子

     下面通过一个简单例子,来引出因子分析背后的思想。

     因子分析的实质是认为m个n维特征的训练样例clip_image069的产生过程如下:

     1、 首先在一个k维的空间中按照多元高斯分布生成m个clip_image071(k维向量),即

     clip_image073

     2、 然后存在一个变换矩阵clip_image075,将clip_image071[1]映射到n维空间中,即

     clip_image077

     因为clip_image071[2]的均值是0,映射后仍然是0。

     3、 然后将clip_image079加上一个均值clip_image008[3](n维),即

     clip_image081

     对应的意义是将变换后的clip_image079[1](n维向量)移动到样本clip_image002[2]的中心点clip_image008[4]

     4、 由于真实样例clip_image002[3]与上述模型生成的有误差,因此我们继续加上误差clip_image083(n维向量),

     而且clip_image083[1]符合多元高斯分布,即

     clip_image085

     clip_image087

     5、 最后的结果认为是真实的训练样例clip_image002[4]的生成公式

     clip_image089

     让我们使用一种直观方法来解释上述过程:

     假设我们有m=5个2维的样本点clip_image002[5](两个特征),如下:

clip_image090

     那么按照因子分析的理解,样本点的生成过程如下:

     1、 我们首先认为在1维空间(这里k=1),存在着按正态分布生成的m个点clip_image092,如下

­clip_image095

     均值为0,方差为1。

     2、 然后使用某个clip_image097将一维的z映射到2维,图形表示如下:

clip_image100

     3、 之后加上clip_image102,即将所有点的横坐标移动clip_image104,纵坐标移动clip_image106,将直线移到一个位置,使得直线过点clip_image008[5],原始左边轴的原点现在为clip_image008[6](红色点)。

clip_image111

     然而,样本点不可能这么规则,在模型上会有一定偏差,因此我们需要将上步生成的点做一些扰动(误差),扰动clip_image113

     4、 加入扰动后,我们得到黑色样本clip_image002[6]如下:

clip_image118

     5、 其中由于z和clip_image119的均值都为0,因此clip_image008[7]也是原始样本点(黑色点)的均值。

     由以上的直观分析,我们知道了因子分析其实就是认为高维样本点实际上是由低维样本点经过高斯分布、线性变换、误差扰动生成的,因此高维数据可以使用低维来表示。

5 因子分析模型

     上面的过程是从隐含随机变量z经过变换和误差扰动来得到观测到的样本点。其中z被称为因子,是低维的。

     我们将式子再列一遍如下:

     clip_image121

     clip_image085[1]

     clip_image123

     其中误差clip_image119[1]和z是独立的。

     下面使用的因子分析表示方法是矩阵表示法,在参考资料中给出了一些其他的表示方法,如果不明白矩阵表示法,可以参考其他资料。

     矩阵表示法认为z和x联合符合多元高斯分布,如下

     clip_image125

     求clip_image127之前需要求E[x]

     clip_image129

           clip_image131

           clip_image133

     我们已知E[z]=0,因此

     clip_image135

     下一步是计算clip_image010[15]

     其中clip_image137

     接着求clip_image139

     clip_image140

     这个过程中利用了z和clip_image119[2]独立假设(clip_image142)。并将clip_image144看作已知变量。

     接着求clip_image146

     clip_image147

     然后得出联合分布的最终形式

     clip_image148

     从上式中可以看出x的边缘分布clip_image150

     那么对样本clip_image152进行最大似然估计

     clip_image154

     然后对各个参数求偏导数不就得到各个参数的值了么?

     可惜我们得不到closed-form。想想也是,如果能得到,还干嘛将z和x放在一起求联合分布呢。根据之前对参数估计的理解,在有隐含变量z时,我们可以考虑使用EM来进行估计。

6 因子分析的EM估计

     我们先来明确一下各个参数,z是隐含变量,clip_image156是待估参数。

     回想EM两个步骤:

循环重复直到收敛 {

(E步)对于每一个i,计算

clip_image158

(M步)计算

clip_image159

     我们套用一下:

     (E步):

     clip_image161

     根据第3节的条件分布讨论,

     clip_image163

     因此

     clip_image164

     那么根据多元高斯分布公式,得到

     clip_image166

     (M步):

     直接写要最大化的目标是

     clip_image167

     其中待估参数是clip_image169

     下面我们重点求clip_image171的估计公式

     首先将上式简化为:

     clip_image173

     这里clip_image175表示clip_image177服从clip_image179分布。然后去掉与clip_image144[1]不相关的项(后两项),得

     clip_image181

     去掉不相关的前两项后,对clip_image144[2]进行导,

     clip_image182

     第一步到第二步利用了tr a = a(a是实数时)和tr AB = tr BA。最后一步利用了    

     tr就是求一个矩阵对角线上元素和。

     最后让其值为0,并且化简得

     clip_image184

     然后得到

     clip_image185

     到这里我们发现,这个公式有点眼熟,与之前回归中的最小二乘法矩阵形式类似

     clip_image186

     这里解释一下两者的相似性,我们这里的x是z的线性函数(包含了一定的噪声)。在E步得到z的估计后,我们找寻的clip_image144[3]实际上是x和z的线性关系。而最小二乘法也是去找特征和结果直接的线性关系。

到这还没完,我们需要求得括号里面的值

     根据我们之前对z|x的定义,我们知道

     clip_image187

     第一步根据z的条件分布得到,第二步根据clip_image189得到

     将上面的结果代入(7)中得到

     clip_image190

     至此,我们得到了clip_image144[4],注意一点是E[z]和clip_image192的不同,后者需要求z的协方差。

     其他参数的迭代公式如下:

     clip_image193

     均值clip_image008[8]在迭代过程中值不变。

     clip_image195

     然后将clip_image197上的对角线上元素抽取出来放到对应的clip_image199中,就得到了clip_image199[1]

7 总结

     根据上面的EM的过程,要对样本X进行因子分析,只需知道要分解的因子数(z的维度)即可。通过EM,我们能够得到转换矩阵clip_image144[5]和误差协方差clip_image199[2]

     因子分析实际上是降维,在得到各个参数后,可以求得z。但是z的各个参数含义需要自己去琢磨。

     下面从一个ppt中摘抄几段话来进一步解释因子分析。

     因子分析(factor analysis)是一种数据简化的技术。它通过研究众多变量之间的内部依赖关系,探求观测数据中的基本结构,并用少数几个假想变量来表示其基本的数据结构。这几个假想变量能够反映原来众多变量的主要信息。原始的变量是可观测的显在变量,而假想变量是不可观测的潜在变量,称为因子。

     例如,在企业形象或品牌形象的研究中,消费者可以通过一个有24个指标构成的评价体系,评价百货商场的24个方面的优劣。

     但消费者主要关心的是三个方面,即商店的环境、商店的服务和商品的价格。因子分析方法可以通过24个变量,找出反映商店环境、商店服务水平和商品价格的三个潜在的因子,对商店进行综合评价。而这三个公共因子可以表示为:

     clip_image201

     这里的clip_image203就是样例x的第i个分量,clip_image207就是clip_image008[9]的第i个分量,clip_image209就是clip_image144[6]的第i行第j列元素,clip_image211是z的第i个分量,clip_image213clip_image215

clip_image211[1]是不可观测的潜在因子。24个变量共享这三个因子,但是每个变量又有自己的个性,不被包含的部分clip_image213[1],称为特殊因子。

     注:

     因子分析与回归分析不同,因子分析中的因子是一个比较抽象的概念,而回归因子有非常明确的实际意义;

     主成分分析分析与因子分析也有不同,主成分分析仅仅是变量变换,而因子分析需要构造因子模型。

     主成分分析:原始变量的线性组合表示新的综合变量,即主成分;

     因子分析:潜在的假想变量和随机影响变量的线性组合表示原始变量。

     PPT地址

     http://www.math.zju.edu.cn/webpagenew/uploadfiles/attachfiles/2008123195228555.ppt

     其他值得参考的文献

     An Introduction to Probabilistic Graphical Models by Jordan Chapter 14

     主成分分析和因子分析的区别http://cos.name/old/view.php?tid=10&id=82

四.ICA

1. 问题:

     1、PCA是一种数据降维的方法,但是只对符合高斯分布的样本点比较有效,那么对于其他分布的样本,有没有主元分解的方法呢?

     2、经典的鸡尾酒宴会问题(cocktail party problem)。假设在party中有n个人,他们可以同时说话,我们也在房间中一些角落里共放置了n个声音接收器(Microphone)用来记录声音。宴会过后,我们从n个麦克风中得到了一组数据clip_image002,i表示采样的时间顺序,也就是说共得到了m组采样,每一组采样都是n维的。我们的目标是单单从这m组采样数据中分辨出每个人说话的信号。

     将第二个问题细化一下,有n个信号源clip_image004clip_image006,每一维都是一个人的声音信号,每个人发出的声音信号独立。A是一个未知的混合矩阵(mixing matrix),用来组合叠加信号s,那么

     clip_image008

     x的意义在上文解释过,这里的x不是一个向量,是一个矩阵。其中每个列向量是clip_image010clip_image012

     表示成图就是

     clip_image014

     这张图来自

     http://amouraux.webnode.com/research-interests/research-interests-erp-analysis/blind-source-separation-bss-of-erps-using-independent-component-analysis-ica/

     clip_image033

     clip_image035的每个分量都由clip_image037的分量线性表示。A和s都是未知的,x是已知的,我们要想办法根据x来推出s。这个过程也称作为盲信号分离。

     令clip_image039,那么clip_image041

     将W表示成

     clip_image042

     其中clip_image044,其实就是将clip_image046写成行向量形式。那么得到:

     clip_image048

2. ICA的不确定性(ICA ambiguities)

     由于w和s都不确定,那么在没有先验知识的情况下,无法同时确定这两个相关参数。比如上面的公式s=wx。当w扩大两倍时,s只需要同时扩大两倍即可,等式仍然满足,因此无法得到唯一的s。同时如果将人的编号打乱,变成另外一个顺序,如上图的蓝色节点的编号变为3,2,1,那么只需要调换A的列向量顺序即可,因此也无法单独确定s。这两种情况称为原信号不确定。

     还有一种ICA不适用的情况,那就是信号不能是高斯分布的。假设只有两个人发出的声音信号符合多值正态分布,clip_image050,I是2*2的单位矩阵,s的概率密度函数就不用说了吧,以均值0为中心,投影面是椭圆的山峰状(参见多值高斯分布)。因为clip_image052,因此,x也是高斯分布的,均值为0,协方差为clip_image054

     令R是正交阵clip_image056clip_image058。如果将A替换成A’。那么clip_image060。s分布没变,因此x’仍然是均值为0,协方差clip_image062

     因此,不管混合矩阵是A还是A’,x的分布情况是一样的,那么就无法确定混合矩阵,也就无法确定原信号。

3. 密度函数和线性变换

     在讨论ICA具体算法之前,我们先来回顾一下概率和线性代数里的知识。

     假设我们的随机变量s有概率密度函数clip_image064(连续值是概率密度函数,离散值是概率)。为了简单,我们再假设s是实数,还有一个随机变量x=As,A和x都是实数。令clip_image066是x的概率密度,那么怎么求clip_image066[1]

     令clip_image039[1],首先将式子变换成clip_image068,然后得到clip_image070,求解完毕。可惜这种方法是错误的。比如s符合均匀分布的话(clip_image072),那么s的概率密度是clip_image074,现在令A=2,即x=2s,也就是说x在[0,2]上均匀分布,可知clip_image076。然而,前面的推导会得到clip_image078。正确的公式应该是

     clip_image080

     推导方法

     clip_image082

     clip_image084

     更一般地,如果s是向量,A可逆的方阵,那么上式子仍然成立。

4. ICA算法

     ICA算法归功于Bell和Sejnowski,这里使用最大似然估计来解释算法,原始的论文中使用的是一个复杂的方法Infomax principal。

     我们假定每个clip_image086有概率密度clip_image088,那么给定时刻原信号的联合分布就是

     clip_image090

     这个公式代表一个假设前提:每个人发出的声音信号各自独立。有了p(s),我们可以求得p(x)

     clip_image092

     左边是每个采样信号x(n维向量)的概率,右边是每个原信号概率的乘积的|W|倍。

     前面提到过,如果没有先验知识,我们无法求得W和s。因此我们需要知道clip_image094,我们打算选取一个概率密度函数赋给s,但是我们不能选取高斯分布的密度函数。在概率论里我们知道密度函数p(x)由累计分布函数(cdf)F(x)求导得到。F(x)要满足两个性质是:单调递增和在[0,1]。我们发现sigmoid函数很适合,定义域负无穷到正无穷,值域0到1,缓慢递增。我们假定s的累积分布函数符合sigmoid函数

     clip_image096

     求导后

     clip_image098

     这就是s的密度函数。这里s是实数。

     如果我们预先知道s的分布函数,那就不用假设了,但是在缺失的情况下,sigmoid函数能够在大多数问题上取得不错的效果。由于上式中clip_image100是个对称函数,因此E[s]=0(s的均值为0),那么E[x]=E[As]=0,x的均值也是0。

     知道了clip_image100[1],就剩下W了。给定采样后的训练样本clip_image002[1],样本对数似然估计如下:

     使用前面得到的x的概率密度函数,得

     clip_image101

     大括号里面是clip_image103

     接下来就是对W求导了,这里牵涉一个问题是对行列式|W|进行求导的方法,属于矩阵微积分。这里先给出结果,在文章最后再给出推导公式。

     clip_image105

     最终得到的求导后公式如下,clip_image107的导数为clip_image109(可以自己验证):

     clip_image110

     其中clip_image112是梯度上升速率,人为指定。

     当迭代求出W后,便可得到clip_image114来还原出原始信号。

     注意:我们计算最大似然估计时,假设了clip_image116clip_image118之间是独立的,然而对于语音信号或者其他具有时间连续依赖特性(比如温度)上,这个假设不能成立。但是在数据足够多时,假设独立对效果影响不大,同时如果事先打乱样例,并运行随机梯度上升算法,那么能够加快收敛速度。

     回顾一下鸡尾酒宴会问题,s是人发出的信号,是连续值,不同时间点的s不同,每个人发出的信号之间独立(clip_image086[1]clip_image120之间独立)。s的累计概率分布函数是sigmoid函数,但是所有人发出声音信号都符合这个分布。A(W的逆阵)代表了s相对于x的位置变化,x是s和A变化后的结果。

5. 实例

     clip_image122

     s=2时的原始信号

     clip_image124

     观察到的x信号

     clip_image126

     使用ICA还原后的s信号

6. 行列式的梯度

     对行列式求导,设矩阵A是n×n的,我们知道行列式与代数余子式有关,

     clip_image127

     clip_image129是去掉第i行第j列后的余子式,那么对clip_image131求导得

     clip_image132

     adj(A)跟我们线性代数中学的clip_image134是一个意思,因此

 clip_image135

7. ICA算法扩展描述

     上面介绍的内容基本上是讲义上的,与我看的另一篇《Independent Component Analysis:

Algorithms and Applications》(Aapo Hyvärinen and Erkki Oja)有点出入。下面总结一下这篇文章里提到的一些内容(有些我也没看明白)。

     首先里面提到了一个与“独立”相似的概念“不相关(uncorrelated)”。Uncorrelated属于部分独立,而不是完全独立,怎么刻画呢?

     如果随机变量clip_image002clip_image004是独立的,当且仅当clip_image006

     如果随机变量clip_image002[1]clip_image004[1]是不相关的,当且仅当clip_image008

     第二个不相关的条件要比第一个独立的条件“松”一些。因为独立能推出不相关,不相关推不出独立。

     证明如下:

     clip_image009

     clip_image010

     clip_image012

     反过来不能推出。

     比如,clip_image002[2]clip_image004[2]的联合分布如下(0,1),(0,-1),(1,0),(-1,0)。

     clip_image014

     因此clip_image002[3]clip_image004[3]不相关,但是

     clip_image016

     因此clip_image002[4]clip_image004[4]不满足上面的积分公式,clip_image002[5]clip_image004[5]不是独立的。

     上面提到过,如果clip_image018是高斯分布的,A是正交的,那么clip_image020也是高斯分布的,且clip_image020[1]clip_image022之间是独立的。那么无法确定A,因为任何正交变换都可以让clip_image020[2]达到同分布的效果。但是如果clip_image018[1]中只有一个分量是高斯分布的,仍然可以使用ICA。

     那么ICA要解决的问题变为:如何从x中推出s,使得s最不可能满足高斯分布?

     中心极限定理告诉我们:大量独立同分布随机变量之和满足高斯分布。

     clip_image023

     我们一直假设的是clip_image020[3]是由独立同分布的主元clip_image018[2]经过混合矩阵A生成。那么为了求clip_image018[3],我们需要计算clip_image018[4]的每个分量clip_image025。定义clip_image027,那么clip_image029,之所以这么麻烦再定义z是想说明一个关系,我们想通过整出一个clip_image031来对clip_image020[4]进行线性组合,得出y。而我们不知道得出的y是否是真正的s的分量,但我们知道y是s的真正分量的线性组合。由于我们不能使s的分量成为高斯分布,因此我们的目标求是让y(也就是clip_image033)最不可能是高斯分布时的w。

     那么问题递归到如何度量y是否是高斯分布的了。

     一种度量方法是kurtosis方法,公式如下:

     clip_image034

     如果y是高斯分布,那么该函数值为0,否则绝大多数情况下值不为0。

     但这种度量方法不怎么好,有很多问题。看下一种方法:

     负熵(Negentropy)度量方法。

     我们在信息论里面知道对于离散的随机变量Y,其熵是

     clip_image035

     连续值时是

     clip_image036

     在信息论里有一个强有力的结论是:高斯分布的随机变量是同方差分布中熵最大的。也就是说对于一个随机变量来说,满足高斯分布时,最随机。

     定义负熵的计算公式如下:

     clip_image037

     也就是随机变量y相对于高斯分布时的熵差,这个公式的问题就是直接计算时较为复杂,一般采用逼近策略。

     clip_image038

     这种逼近策略不够好,作者提出了基于最大熵的更优的公式:

     clip_image039

     之后的FastICA就基于这个公式。

     另外一种度量方法是最小互信息方法:

     clip_image040

     这个公式可以这样解释,前一个H是clip_image042的编码长度(以信息编码的方式理解),第二个H是y成为随机变量时的平均编码长度。之后的内容包括FastICA就不再介绍了,我也没看懂。

8. ICA的投影追踪解释(Projection Pursuit)

     投影追踪在统计学中的意思是去寻找多维数据的“interesting”投影。这些投影可用在数据可视化、密度估计和回归中。比如在一维的投影追踪中,我们寻找一条直线,使得所有的数据点投影到直线上后,能够反映出数据的分布。然而我们最不想要的是高斯分布,最不像高斯分布的数据点最interesting。这个与我们的ICA思想是一直的,寻找独立的最不可能是高斯分布的s。

     在下图中,主元是纵轴,拥有最大的方差,但最interesting的是横轴,因为它可以将两个类分开(信号分离)。

     clip_image043

9. ICA算法的前处理步骤

     1、中心化:也就是求x均值,然后让所有x减去均值,这一步与PCA一致。

     2、漂白:目的是将x乘以一个矩阵变成clip_image045,使得clip_image045[1]的协方差矩阵是clip_image047。解释一下吧,原始的向量是x。转换后的是clip_image045[2]

     clip_image045[3]的协方差矩阵是clip_image047[1],即

     clip_image048

     我们只需用下面的变换,就可以从x得到想要的clip_image045[4]

     clip_image049

     其中使用特征值分解来得到E(特征向量矩阵)和D(特征值对角矩阵),计算公式为

     clip_image050

     下面用个图来直观描述一下:

     假设信号源s1和s2是独立的,比如下图横轴是s1,纵轴是s2,根据s1得不到s2。

     clip_image051

     我们只知道他们合成后的信号x,如下

     image

     此时x1和x2不是独立的(比如看最上面的尖角,知道了x1就知道了x2)。那么直接代入我们之前的极大似然概率估计会有问题,因为我们假定x是独立的。

     因此,漂白这一步为了让x独立。漂白结果如下:

     clip_image053

     可以看到数据变成了方阵,在clip_image045[5]的维度上已经达到了独立。

     然而这时x分布很好的情况下能够这样转换,当有噪音时怎么办呢?可以先使用前面提到的PCA方法来对数据进行降维,滤去噪声信号,得到k维的正交向量,然后再使用ICA。

10. 小结

     ICA的盲信号分析领域的一个强有力方法,也是求非高斯分布数据隐含因子的方法。从之前我们熟悉的样本-特征角度看,我们使用ICA的前提条件是,认为样本数据由独立非高斯分布的隐含因子产生,隐含因子个数等于特征数,我们要求的是隐含因子。

     而PCA认为特征是由k个正交的特征(也可看作是隐含因子)生成的,我们要求的是数据在新特征上的投影。同是因子分析,一个用来更适合用来还原信号(因为信号比较有规律,经常不是高斯分布的),一个更适合用来降维(用那么多特征干嘛,k个正交的即可)。有时候也需要组合两者一起使用。这段是我的个人理解,仅供参考。

五.MDS
假设给定N个点,切每对点之间的距离dij.但我们不知道这些点的确切坐标,也不知道他们的维度,以及该距离是如何计算的。多维定标(MDS)是把这些点映射到低维空间(如2维)的方法,使他们在低维空间的欧式距离尽可能的接近原空间给定的距离dij.实际上,在相关矩阵上而不是协方差矩阵上做PCA等价于标准化的欧式距离做MDS,其中每个变量有单位方差。这里也不做详细解释。


六.线性判别分析(LDA)

1. 问题

     之前我们讨论的PCA、ICA也好,对样本数据来言,可以是没有类别标签y的。回想我们做回归时,如果特征太多,那么会产生不相关特征引入、过度拟合等问题。我们可以使用PCA来降维,但PCA没有将类别标签考虑进去,属于无监督的。

     比如回到上次提出的文档中含有“learn”和“study”的问题,使用PCA后,也许可以将这两个特征合并为一个,降了维度。但假设我们的类别标签y是判断这篇文章的topic是不是有关学习方面的。那么这两个特征对y几乎没什么影响,完全可以去除。

     再举一个例子,假设我们对一张100*100像素的图片做人脸识别,每个像素是一个特征,那么会有10000个特征,而对应的类别标签y仅仅是0/1值,1代表是人脸。这么多特征不仅训练复杂,而且不必要特征对结果会带来不可预知的影响,但我们想得到降维后的一些最佳特征(与y关系最密切的),怎么办呢?

2. 线性判别分析(二类情况)

     回顾我们之前的logistic回归方法,给定m个n维特征的训练样例clip_image002(i从1到m),每个clip_image004对应一个类标签clip_image006。我们就是要学习出参数clip_image008,使得clip_image010(g是sigmoid函数)。

     现在只考虑二值分类情况,也就是y=1或者y=0。

     为了方便表示,我们先换符号重新定义问题,给定特征为d维的N个样例,clip_image012,其中有clip_image014个样例属于类别clip_image016,另外clip_image018个样例属于类别clip_image020

     现在我们觉得原始特征数太多,想将d维特征降到只有一维,而又要保证类别能够“清晰”地反映在低维数据上,也就是这一维就能决定每个样例的类别。

     我们将这个最佳的向量称为w(d维),那么样例x(d维)到w上的投影可以用下式来计算

     clip_image022

     这里得到的y值不是0/1值,而是x投影到直线上的点到原点的距离。

     当x是二维的,我们就是要找一条直线(方向为w)来做投影,然后寻找最能使样本点分离的直线。如下图:

     clip_image024

     从直观上来看,右图比较好,可以很好地将不同类别的样本点分离。

     接下来我们从定量的角度来找到这个最佳的w。

     首先我们寻找每类样例的均值(中心点),这里i只有两个

     clip_image026

     由于x到w投影后的样本点均值为

     clip_image028

     由此可知,投影后的的均值也就是样本中心点的投影。

     什么是最佳的直线(w)呢?我们首先发现,能够使投影后的两类样本中心点尽量分离的直线是好的直线,定量表示就是:

     clip_image030

     J(w)越大越好。

     但是只考虑J(w)行不行呢?不行,看下图

     clip_image031

     样本点均匀分布在椭圆里,投影到横轴x1上时能够获得更大的中心点间距J(w),但是由于有重叠,x1不能分离样本点。投影到纵轴x2上,虽然J(w)较小,但是能够分离样本点。因此我们还需要考虑样本点之间的方差,方差越大,样本点越难以分离。

     我们使用另外一个度量值,称作散列值(scatter),对投影后的类求散列值,如下

     clip_image033

     从公式中可以看出,只是少除以样本数量的方差值,散列值的几何意义是样本点的密集程度,值越大,越分散,反之,越集中。

     而我们想要的投影后的样本点的样子是:不同类别的样本点越分开越好,同类的越聚集越好,也就是均值差越大越好,散列值越小越好。正好,我们可以使用J(w)和S来度量,最终的度量公式是

     clip_image035

     接下来的事就比较明显了,我们只需寻找使J(w)最大的w即可。

     先把散列值公式展开

     clip_image037

     我们定义上式中中间那部分

     clip_image039

     这个公式的样子不就是少除以样例数的协方差矩阵么,称为散列矩阵(scatter matrices)

     我们继续定义

     clip_image041

     clip_image043称为Within-class scatter matrix。

     那么回到上面clip_image045的公式,使用clip_image047替换中间部分,得

     clip_image049

     clip_image051

     然后,我们展开分子

     clip_image052

     clip_image054称为Between-class scatter,是两个向量的外积,虽然是个矩阵,但秩为1。

     那么J(w)最终可以表示为

     clip_image056

     在我们求导之前,需要对分母进行归一化,因为不做归一的话,w扩大任何倍,都成立,我们就无法确定w。因此我们打算令clip_image058,那么加入拉格朗日乘子后,求导

     clip_image059

     其中用到了矩阵微积分,求导时可以简单地把clip_image061当做clip_image063看待。

     如果clip_image043[1]可逆,那么将求导后的结果两边都乘以clip_image065,得

     clip_image066

     这个可喜的结果就是w就是矩阵clip_image068的特征向量了。

     这个公式称为Fisher linear discrimination。

     等等,让我们再观察一下,发现前面clip_image070的公式

     clip_image072

     那么

     clip_image074

     代入最后的特征值公式得

     clip_image076

     由于对w扩大缩小任何倍不影响结果,因此可以约去两边的未知常数clip_image078clip_image080,得到

     clip_image082

     至此,我们只需要求出原始样本的均值和方差就可以求出最佳的方向w,这就是Fisher于1936年提出的线性判别分析。

     看上面二维样本的投影结果图:

     clip_image083

3. 线性判别分析(多类情况)

     前面是针对只有两个类的情况,假设类别变成多个了,那么要怎么改变,才能保证投影后类别能够分离呢?

我们之前讨论的是如何将d维降到一维,现在类别多了,一维可能已经不能满足要求。假设我们有C个类别,需要K维向量(或者叫做基向量)来做投影。

     将这K维向量表示为clip_image085

     我们将样本点在这K维向量投影后结果表示为clip_image087,有以下公式成立

     clip_image089

     clip_image091

     为了像上节一样度量J(w),我们打算仍然从类间散列度和类内散列度来考虑。

     当样本是二维时,我们从几何意义上考虑:

     clip_image092

     其中clip_image094clip_image043[2]与上节的意义一样,clip_image096是类别1里的样本点相对于该类中心点clip_image098的散列程度。clip_image100变成类别1中心点相对于样本中心点clip_image102的协方差矩阵,即类1相对于clip_image102[1]的散列程度。

     clip_image043[3]

     clip_image104

     clip_image106的计算公式不变,仍然类似于类内部样本点的协方差矩阵

     clip_image108

     clip_image054[1]需要变,原来度量的是两个均值点的散列情况,现在度量的是每类均值点相对于样本中心的散列情况。类似于将clip_image094[1]看作样本点,clip_image102[2]是均值的协方差矩阵,如果某类里面的样本点较多,那么其权重稍大,权重用Ni/N表示,但由于J(w)对倍数不敏感,因此使用Ni。

     clip_image110

     其中

     clip_image112

     clip_image102[3]是所有样本的均值。

     上面讨论的都是在投影前的公式变化,但真正的J(w)的分子分母都是在投影后计算的。下面我们看样本点投影后的公式改变:

     这两个是第i类样本点在某基向量上投影后的均值计算公式。

     clip_image114

     clip_image116

     下面两个是在某基向量上投影后的clip_image043[4]clip_image070[1]

     clip_image118

     clip_image120

     其实就是将clip_image102[4]换成了clip_image122

     综合各个投影向量(w)上的clip_image124clip_image126,更新这两个参数,得到

     clip_image128

     clip_image130

     W是基向量矩阵,clip_image124[1]是投影后的各个类内部的散列矩阵之和,clip_image126[1]是投影后各个类中心相对于全样本中心投影的散列矩阵之和。

     回想我们上节的公式J(w),分子是两类中心距,分母是每个类自己的散列度。现在投影方向是多维了(好几条直线),分子需要做一些改变,我们不是求两两样本中心距之和(这个对描述类别间的分散程度没有用),而是求每类中心相对于全样本中心的散列度之和。

     然而,最后的J(w)的形式是

     clip_image132

     由于我们得到的分子分母都是散列矩阵,要将矩阵变成实数,需要取行列式。又因为行列式的值实际上是矩阵特征值的积,一个特征值可以表示在该特征向量上的发散程度。因此我们使用行列式来计算(此处我感觉有点牵强,道理不是那么有说服力)。

     整个问题又回归为求J(w)的最大值了,我们固定分母为1,然后求导,得出最后结果(我翻查了很多讲义和文章,没有找到求导的过程)

     clip_image134

     与上节得出的结论一样

     clip_image136

     最后还归结到了求矩阵的特征值上来了。首先求出clip_image138的特征值,然后取前K个特征向量组成W矩阵即可。

     注意:由于clip_image070[2]中的clip_image140 秩为1,因此clip_image070[3]的秩至多为C(矩阵的秩小于等于各个相加矩阵的秩的和)。由于知道了前C-1个clip_image094[2]后,最后一个clip_image142可以有前面的clip_image094[3]来线性表示,因此clip_image070[4]的秩至多为C-1。那么K最大为C-1,即特征向量最多有C-1个。特征值大的对应的特征向量分割性能最好。

     由于clip_image138[1]不一定是对称阵,因此得到的K个特征向量不一定正交,这也是与PCA不同的地方。

4. 实例

      将3维空间上的球体样本点投影到二维上,W1相比W2能够获得更好的分离效果。

      clip_image002

      PCA与LDA的降维对比:

      clip_image004

      PCA选择样本点投影具有最大方差的方向,LDA选择分类性能最好的方向。

      LDA既然叫做线性判别分析,应该具有一定的预测功能,比如新来一个样例x,如何确定其类别?

      拿二值分来来说,我们可以将其投影到直线上,得到y,然后看看y是否在超过某个阈值y0,超过是某一类,否则是另一类。而怎么寻找这个y0呢?

      看

      clip_image006

      根据中心极限定理,独立同分布的随机变量和符合高斯分布,然后利用极大似然估计求

      clip_image008

      然后用决策理论里的公式来寻找最佳的y0,详情请参阅PRML。

      这是一种可行但比较繁琐的选取方法,可以看第7节(一些问题)来得到简单的答案。

5. 使用LDA的一些限制

      1、 LDA至多可生成C-1维子空间

      LDA降维后的维度区间在[1,C-1],与原始特征数n无关,对于二值分类,最多投影到1维。

      2、 LDA不适合对非高斯分布样本进行降维。

      clip_image010

      上图中红色区域表示一类样本,蓝色区域表示另一类,由于是2类,所以最多投影到1维上。不管在直线上怎么投影,都难使红色点和蓝色点内部凝聚,类间分离。

      3、 LDA在样本分类信息依赖方差而不是均值时,效果不好。

      clip_image011

      上图中,样本点依靠方差信息进行分类,而不是均值信息。LDA不能够进行有效分类,因为LDA过度依靠均值信息。

      4、 LDA可能过度拟合数据。

6. LDA的一些变种

1、 非参数LDA

      非参数LDA使用本地信息和K临近样本点来计算clip_image013,使得clip_image013[1]是全秩的,这样我们可以抽取多余C-1个特征向量。而且投影后分离效果更好。

2、 正交LDA

      先找到最佳的特征向量,然后找与这个特征向量正交且最大化fisher条件的向量。这种方法也能摆脱C-1的限制。

3、 一般化LDA

      引入了贝叶斯风险等理论

4、 核函数LDA

      将特征clip_image015,使用核函数来计算。

7. 一些问题

      上面在多值分类中使用的

      clip_image017

      是带权重的各类样本中心到全样本中心的散列矩阵。如果C=2(也就是二值分类时)套用这个公式,不能够得出在二值分类中使用的clip_image013[2]

      clip_image019

      因此二值分类和多值分类时求得的clip_image013[3]会不同,而clip_image021意义是一致的。

      对于二值分类问题,令人惊奇的是最小二乘法和Fisher线性判别分析是一致的。

      下面我们证明这个结论,并且给出第4节提出的y0值得选取问题。

      回顾之前的线性回归,给定N个d维特征的训练样例clip_image023(i从1到N),每个clip_image025对应一个类标签clip_image027。我们之前令y=0表示一类,y=1表示另一类,现在我们为了证明最小二乘法和LDA的关系,我们需要做一些改变

      clip_image029

      就是将0/1做了值替换。

      我们列出最小二乘法公式

      clip_image031

      w和clip_image033是拟合权重参数。

      分别对clip_image033[1]和w求导得

      clip_image035

      clip_image037

      从第一个式子展开可以得到

      clip_image039

      消元后,得

      clip_image041

      clip_image043

      可以证明第二个式子展开后和下面的公式等价

      clip_image045

      其中clip_image047clip_image049与二值分类中的公式一样。

      由于clip_image051

      因此,最后结果仍然是

      clip_image053

      这个过程从几何意义上去理解也就是变形后的线性回归(将类标签重新定义),线性回归后的直线方向就是二值分类中LDA求得的直线方向w。

      好了,我们从改变后的y的定义可以看出y>0属于类clip_image055,y<0属于类clip_image057。因此我们可以选取y0=0,即如果clip_image059,就是类clip_image055[1],否则是类clip_image057[1]


这里主要介绍线性判别式分析(LDA),主要基于Fisher Discriminant Analysis with Kernals[1]和Fisher Linear Discriminant Analysis[2]两篇文献。

  LDA与PCA的一大不同点在于,LDA是有监督的算法,而PCA是无监督的,因为PCA算法没有考虑数据的标签(类别),只是把原数据映射到一些方差比较大的方向(基)上去而已。而LDA算法则考虑了数据的标签。文献[2]中举了一个非常形象的例子,说明了在有些情况下,PCA算法的性能很差,如下图:                                   

 我们用不同的颜色标注C1,C2两个不同类别的数据。根据PCA算法,数据应该映射到方差最大的那个方向,亦即Y轴方向,但是如果映射到Y轴方向,C1,C2两个不同类别的数据将完全混合在一起,很难区分开,所以使用PCA算法进行降维后再进行分类的效果会非常差。但是使用LDA算法,数据会映射到X轴方向。

  LDA算法会考虑到数据的类别属性,给定两个类别C1、C2,我们希望找到一个向量ω,当数据映射到ω的方向上时,来自两个类的数据尽可能的分开,同一个类内的数据尽可能的紧凑。数据的映射公式为:z=ωTx,  其中z是数据x到ω上的投影,因而也是一个d维到1维的维度归约。

  令 m 1和m 1分别表示C1类数据投影之前个投影之后的均值,易知m 1=ωTm1,同理m2=ωTm2

  令s 1 2和s22分别表示C1和C2类数据在投影之后的散布(scatter),亦即s12=∑(ωTxt-m1)2rt,s22=∑(ωTxt-m2)2(1-rt)其中如果xt∈C1,则rt=1,否则rt=0。

  我们希望|m 1-m 2|尽可能的大,而s12+s22尽可能的小,Fisher线性判别式就是最大化下面式子的ω:

  J(ω)=(m 1-m 2)2/(s12+s22)     式子-1

  改写式子-1中的分子:  (m 1-m 2)2=  (ωTm1-ωTm2)2=ωT(m1-m2)(m1-m2)Tω=ωTSBω

  其中 S B=( m 1- m 2)(m1-m2)T    式子-2

  是 类间散布矩阵(between class scatter matrix)。

  改写式子-1中的分母:

  ∑(ω Tx t-m1)2rt=∑ωT(xt-m1)(xt-m1)Tωrt=ωTS1ω, 其中S 1=∑rt(xt-m1)(xt-m1)T是C1的类内散布矩阵(within class scatter matrix)。

  令 S W= S 1+ S 2,是 类内散布的总和,则s 12+s22=ωTSWω。

  所以式子-1可以改写为:

  J(ω)=(ω T S Bω)/(ω TSWω)    式子-3

  我们只需要使式子-3对于ω求导,然后使导数等于0,便可以求出ω的值:ω=c S W-1(m1-m2),其中c是一个参数,我们只对ω的方向感兴趣,所以c可以取值为1.

  另外,最后求得的 J(ω)的值等于λ k,λ k是SW-1SB的最大的特征值,而ω则是SW-1SB的最大特征值所对应的特征向量。

  最后有一些关于LDA算法的讨论,出自文献[1]:

  1. Fisher LDA对数据的分布做了一些很强的假设,比如每个类的数据都是高斯分布,各个类的协方差相等。虽然这些强假设很可能在实际数据中并不满足,但是Fisher LDA已经被证明是非常有效地降维算法,其中的原因是线性模型对于噪音的鲁棒性比较好,不容易过拟合。

  2. 准确的估计数据的散布矩阵是非常重要的,很可能会有较大的偏置。用式子-2进行估计在样本数据比较少(相对于维数来说)时会产生较大的变异性。

  参考文献:

  [1] Fisher Discriminant Analysis with Kernals. Sebastian Mika, Gunnar Ratsch, Jason Weston, Bernhadr Scholkopf, Klaus-Robert Muller.

  [2] Fisher Linear Discriminant Analysis. Max Welling.

  [3] 机器学习导论。 Ethem Alpaydin

七.局部线性嵌入LLE

LLE(Locally Linear Embedding)算法,即局部线性嵌入算法。该算法是针对非线性信号特征矢量维数的优化方法,这种维数优化并不是仅仅在数量上简单的约简,而是在保持原始数据性质不变的情况下,将高维空间的信号映射到低维空间上,即特征值的二次提取[1] 

理论的提出

LLE 是 S.T.Roweis 等人提出了一种针对非线性数据的无监督降维方法,它是流行学习算法中的一种用局部线性反映全局的非线性的算法,并能够使降维的数据保持原有数据的拓扑结构。LLE 是一种用局部线性的非线性维数约简方法,该算法在一定程度上扩大了对维数约简的认识。[2] 

自从在 2000 年 Roweis 提出了局部线性嵌入(LLE)算法后,使得流行学习在最近几年里得到迅速发展,流行学习与传统的线性降维方法相比,能够非常有效地发现非线性高维空间的维数,有利于对高维数据的维数约简和数据分析,在机器学习和认知科学领域受到广大研究者的重视。[3] 
基本原理
LLE 算法假设在局部领域内数据点是线性的,所以邻域内任意一点,都可用局部近邻点来线性表示。LLE 算法是由重构成本函数最小化求出最优权值,各点的局部邻域权值能够在多尺度变换下仍保持不变。LLE 算法无迭代计算过程,可使计算复杂度大幅度的减小。
假设样本集由 N 个 D 维矢量Xi 组成。每一样本点都可以用它的近邻域点带权值线性组合表示,其中权值能反映出局部邻域的信息。根据这些信息,可使这种低维空间中仍保留原高维空间中的几何性质,通过重叠的局部邻域,得到整体的信息。这种方法的关键就在于:用局部线性结构来反映全局的非线性结构。在每一小块局部上,流行可近似的看成是平坦的,故这种算法不会产生较大的错误
将 N 个 D 维矢量组成的矩阵作为输入,输出矢量为一个 N 个 d 维矢量(d<D)组成的矩阵。矩阵 Y 的第 k 列对应是矩阵 X 中的第 k 列。
优缺点
优点
(1)LLE 算法能够突破主元分析法在非线性数据的局限,可以处理、分析非线性信号。
(2)该算法可以很好表达数据的内在流行结构,能够保留数据的本质特征,这样可以很好的保留原有数据特征,这点在故障诊断中有重要的意义。
(3)该算法本身参数的选择很少,故能更好的进行特征参数优化,这为故障检测和故障诊断打下坚实的基础。
缺点
(1)LLE 算法需要进行稠密采样;
(2)该算法的局部邻域参数 k、嵌入维数 d 和信号中的噪声,会影响高维空间的降维效果;
(3)LLE 无法处理等距流行等。[4] 
主要应用
LLE 主要应用在人脸识别图像处理等领域中,在故障诊断中应用较少,并且大部分是在信号处理中的应用。
参考资料
  • 1.  龚明.基于局部线性嵌入算法的柴油机故障诊断研究[硕士学位论文].太原:中北大学,2013.
  • 2.  Mclvor R T ,Humphreys P K.A case-based reasoning approach to the make or buy decision[J].Integrated Manufacturing Systems,2000(5):295~310.
  • 3.  王珏,周志华,周傲英.机器学习及其应用:清华大学出版社,2000:135-139
  • 4.  罗方琼.LLE 流行学习的若干问题分析[J].研究与开发,2012(3):13~16

八.等距特征映射

实际上,6和7属于流形算法了。这里简单解释下等距特征映射,也叫ISOMAP。

例如,一张人脸100*100的二维图像,在这种情况下,每张人脸是10000维空间里的一个点。假设随着一个人由右向左慢慢转动头部,我们取一系列照片,我们得到这个面部图像序列沿着10000维空间中的一条轨迹,而这条曲线不是线性的。现在考虑许多人的人脸,随着他们转动头部,所有人的人脸轨迹定义了10000维空间中的一个流形,并且这就是我们想要建模的。实际上,两张人脸的相似性,不能简单用像素差的和表示,而是应该计算沿着流形的距离,这称作测地距离。等距特征映射就是估计这种距离,使用多维定标进行纬度归约。

九.小波变换

离散小波变换(DWT)是一种线性信号处理技术,当用于数据向量X时,将它变换成数值上不同的小波系数向量X’。两个向量具有相同的长度。当这种技术用于数据归约时,每个元组看作一个n维数据向量,即X = (x1, x2, ., xn),描述n个数据库属性在元组上的n个测量值。

“如果小波变换的数据与原数据的长度相等,这种技术如何能够用于数据归约?”关键在在我们的记号中,代表向量的变量用粗斜体,描述向量的度量用斜体。

于小波变换后的数据可以截短。仅存放一小部分最强的小波系数,就能保留近似的压缩数据。

例如,保留大于用户设定的某个阈值的所有小波系数,其他系数置为0。这样,结果数据表示非常稀疏,使得如果在小波空间进行计算,利用数据稀疏特点的操作计算得非常快。该技术也能用于消除噪声,而不会光滑掉数据的主要特征,使得它们也能有效地用于数据清理。给定一组系数,使用所用的DWT的逆,可以构造原数据的近似。

DWT与离散傅里叶变换(DFT)有密切关系,DFT是一种涉及正弦和余弦的信号处理技术。然而一般地说,DWT是一种更好的有损压缩。也就是说,对于给定的数据向量,如果DWT和DFT保留相同数目的系数,DWT将提供原数据的更准确的近似。因此,对于等价的近似,DWT比DFT需要的空间小。不像DFT,小波空间局部性相当好,有助于保留局部细节。

只有一种DFT,但有若干族DWT。图2-16显示了一些小波族。流行的小波变换包括Haar-2, Daubechies-4和Daubechies-6变换。应用离散小波变换的一般过程使用一种分层金字塔算法(pyramid algorithm),它在每次迭代将数据减半,导致很快的计算速度。该方法如下:

图2-16 小波族的例子。小波名后的数是小波的消失瞬间。这是系数必须满足的数学联系集,并且与小波系数的个数有关

(1)输入数据向量的长度L必须是2的整数幂。必要时(L≥n),通过在数据向量后添加0, 这一条件可以满足。
(2)每个变换涉及应用两个函数。第一个使用某种数据光滑,如求和或加权平均。第二个进行加权差分,产生数据的细节特征。
(3)两个函数作用于X中的数据点对,即用于所有的测量对(x2i, x2i+1)。这导致两个长度为L/2的数据集。一般,它们分别代表输入数据的光滑后的版本或低频版本和它的高频内容。
(4)两个函数递归地作用于前面循环得到的数据集,直到得到的数据集长度为2。
(5)由以上迭代得到的数据集中选择值,指定其为数据变换的小波系数。等价地,可以将矩阵乘法用于输入数据,以得到小波系数。所用的矩阵依赖于给定的DWT。矩阵必须是标准正交的,即列是单位向量并相互正交,使得矩阵的逆是它的转置。尽管受篇幅限制,这里我们不再讨论,但这种性质允许由光滑和光滑-差数据集重构数据。通过将矩阵因子分解成几个稀疏矩阵,对于长度为n的输入向量,“快速DWT”算法的复杂度为O (n)。

小波变换可以用于多维数据,如数据立方体。可以按以下方法做:首先将变换用于第一个维,然后第二个,如此下去。计算复杂性关于立方体中单元的个数是线性的。对于稀疏或倾斜数据和具有有序属性的数据,小波变换给出很好的结果。据报道,小波变换的有损压缩比当前的商业标准JPEG压缩好。小波变换有许多实际应用,包括指纹图像压缩、计算机视觉、时间序列数据分析和数据清理。

阅读更多
想对作者说点什么?

博主推荐

换一批

没有更多推荐了,返回首页