【学习笔记】一些数据降维方法

810 阅读9分钟

数据降维即将多维数据降为低维,代表性的方法有主成因分析(PCA,principal component analysis)、局部线性嵌入(LLE,locally linear embedding)和拉普拉斯特征映射(Laplacian eigen-map)。

主成因分析(PCA)

在介绍PCA之前,首先需要了解一些线性代数的基本概念:

特征值与特征向量

对于一个n阶矩阵A和实数λ\lambda,如果能找到一个非零向量x满足:

Ax=λxAx = \lambda x

则将λ\lambda称为A的特征值,x称为A的特征向量。

从几何意义上来说,对于矩阵A,存在一些向量x,使得Ax和x的方向没有发生变化,只有长度发生了变化,长度发生的变化可以用系数λ\lambda来度量。

协方差矩阵

协方差

协方差用于衡量两个变量间的误差,方差是协方差的一种特殊情况(两个变量一致时),方差用于衡量单个随机变量的离散程度,其公式如下所示:

σx2=1n1i=1n(xixˉ)2\sigma^2_x=\frac{1}{n-1}\sum^n_{i=1}{(x_i-\bar{x})^2}

协方差则度量两个随机变量的相似程度(变化趋势),其公式如下:

σ(x,y)=1n1i=1n(xixˉ)(yiyˉ)\sigma(x,y)=\frac{1}{n-1}\sum^n_{i=1}{(x_i-\bar{x})(y_i-\bar{y})}
Cov(X,Y)=E[(Xμx)(Yμy)]Cov(X,Y)=E[(X-\mu_x)(Y-\mu_y)]

协方差的计算中,两个变量和自身的期望做差再相乘,当两个变量的取值都大于/小于自身期望时,计算得到的协方差为正值,即两个变量的变化趋势相同,反之若一大一小,则协方差为负值。因此协方差可以反映变量间的相关程度。

在现实中描述一个物体的特征通常是多维的,如描述一个人会包括他的身高、性别、体重等等。在对多维数据进行分析时,就需要协方差矩阵来描述,例如对于2维数据,任意两个维度间求其协方差,则可以得到一个包括4个协方差的协方差矩阵:

Σ=[σ(x,x)σ(x,y)σ(y,x)σ(y,y)](2)\Sigma = \left[ \begin{matrix} \sigma(x,x) & \sigma(x,y)\\ \sigma(y,x) & \sigma(y,y)\\ \end{matrix} \right] \tag{2}

可以看到,主对角线上的元素实际为两个维度的方差(特殊的协方差)。下图为4个不同协方差矩阵的数据分布情况:

image.png

主对角线上的值表示数据在x和y方向的离散程度,值越小越集中。而σ(x,y)\sigma(x,y)σ(y,x)\sigma(y,x)的正负以及值的大小决定了x与y的相关性,协方差的绝对值越大,相关关系越明显。

PCA基本原理

PCA的主要目的是将n维特征映射到比n小的k维上,并最大程度降低其信息损失。其方法是从原始的空间顺序中找到一组相互正交的坐标轴,然后第一个新坐标轴的选择时原始数据中方差最大的方向,第二个是与第一个坐标轴正交的平面中方差最大的方向,以此类推。研究发现,大部分方差都包含在前k个坐标轴中,后面包含的很少,因此可以只报了前k个含有绝大部分方差的坐标轴,从而实现数据特征的降维。

通过计算数据矩阵的协方差矩阵,再计算协方差矩阵的特征值和特征向量,选取特征值最大(即方差最大)的k个特征所对应的特征向量组成新矩阵,就能实现上面的降维方法。

计算协方差矩阵特征值和特征向量有两种方法,因此PCA算法的实现也有两种:

基于特征值分解协方差矩阵

对于n维数据集X={x1,x2,...,xn}X=\{x_1,x_2,...,x_n\}进行降维主要步骤如下所示:

  • 零均值化(中心化):每一个特征减去各自的均值。(关于这一步的解释见链接
  • 计算协方差矩阵1nXXT\frac{1}{n}XX^T。(推导过程可见链接
  • 用特征值分解法求协方差矩阵的特征和特征向量

特征值分解矩阵 若矩阵A是一个m×mm\times m的实对称矩阵(A=ATA=A^T),矩阵A通过特征值分解可以化为下式:

A=QΣQ1=QΣQTA=Q\Sigma Q^{-1}=Q\Sigma Q^{T}

其中Q是矩阵A特征向量组成的矩阵,Σ\Sigma是一个对角阵,对角线上的元素为特征值,QQT=IQQ^T=I。(对于正交矩阵 QT=Q1Q^T=Q^{-1}

上述式子的推导过程见链接

  • 将求得的特征值从大到小排列并选取其中最大的k个,将其对应的k个特征向量作为行向量组成特征向量矩阵P。

  • 将数据转换到k个特征向量组成的新空间中:Y=PXY=PX

基于SVD(奇异值)分解矩阵 基于SVD的方法则是通过奇异值分解来求协方差矩阵的特征值和特征向量。特征值分解要求矩阵A必须为实对称矩阵,但遇到一般性矩阵(如m×nm\times n的矩阵A),则需要用SVD的方法进行分解。

对于矩阵A可以有以下分解形式:

A=UΣVTA=U\Sigma V^T

其中,U和V均为单位正交矩阵,即UUT=I,VVT=IUU^T=I,VV^T=IΣ\Sigma仅在主对角线上有值,称为奇异值。且URm×m,ΣRm×n,VRn×nU\in R^{m\times m}, \Sigma \in R^{m\times n}, V \in R^{n\times n}

推导过程见链接

局部线性嵌入(LLE)

PCA为线性降维方法,但现实中存在很多非线性高维数据,为了有效处理非线性数据,提出了LLE方法。

LLE基本原理

LLE方法基于假设:局部原始数据近似位于一张超平面上,局部内某一数据可以由其领域数据线性表示。数学表达为:

xi=wi1xi1+wi2xi2+...+wikxikx_i = w_{i1} * x_{i1} + w_{i2} * x_{i2} + ...+ w_{ik} * x_{ik}

其中的xix_i表示原始数据集中第i个元素,xikx_{ik}表示x领域内的第k个数据,wikw_{ik}表示对应的权重。在xix_i确定的情况下,可以使用K-NN等算法求其领域数据。

降维后,样本在低维上的投影也应该满足其原有的线性关系,y是x在低维空间上的投影,则y应该满足:

yi=wi1yi1+wi2yi2+...+wikyiky_i = w_{i1} * y_{i1} + w_{i2} * y_{i2} + ...+ w_{ik} * y_{ik}

LLE算法的基本步骤为:

  • 在高维空间中根据样本的近邻样本求权重系数
  • 在权重系数和关系不变的情况下计算样本低维投影

权重计算

对m个n维样本{x1,x2,...,xm}\{x_1,x_2,...,x_m\}xjx_jxix_i的近邻样本,通过最小化损失函数(均方差)来求权重系数:

l(wi)=xij=1kwijxij2l(w_i)=\lVert x_i - \sum^k_{j=1}w_{ij}x_{ij}\rVert^2

对权重添加约束条件:j=1kwij=1\sum^k_{j=1}w_{ij}=1

代入上式变为:

l(wi)=j=1kwijxij=1kwijxij2=j=1kwij(xixij)2(3)=[xixi1,...,xixik]Wi2(4)=WiTQiTQiWil(w_i)=\lVert \sum^k_{j=1}w_{ij}x_i - \sum^k_{j=1}w_{ij}x_{ij}\rVert^2 \\ = \lVert \sum^k_{j=1}w_{ij}(x_i-x_{ij}) \rVert^2 \\ (3)= \lVert [xi-x_{i1},...,x_i-x_{ik}]W_i\rVert ^2 \\ (4)= W_i^TQ_i^TQ_iW_i

其中,Wi=[wi1,...,wik]W_i=[w_{i1},...,w_{ik}]Qi=[xixi1,...,xixik]Q_i =[xi-x_{i1},...,x_i-x_{ik}]

关于公式中(3)→(4)的个人思考,损失函数是在计算均方差,之前说到方差是协方差的特殊情况而协方差又可以表示为1nXXT\frac{1}{n}XX^T,因此(3)可以转化为(4)式。

Ai=QiTQiA_i=Q^T_iQ_i,则可以得到:l(Wi)=WiTAiWil(W_i)=W^T_iA_iW_i

将权重的约束矩阵化为:j=1kwij=WiT1k=1\sum^k_{j=1}w_{ij}=W^T_i1_k=1,其中1k1_k为k维的全1向量。

对于带约束的优化问题,通常用拉格朗日乘子法(Lagrange Multiplier),将约束式用一个系数与原式合并为一个式子(拉格朗日函数),从而消除约束,再对拉格朗日函数求导求极值极值验证得最优解。(具体见链接

使用拉格朗日乘子法得到:

l(Wi)=WiTAiWi+λ(WiT1k1)l(W_i)=W^T_iA_iW_i+\lambda(W^T_i1_k-1)

对上式的W求偏导得:

lWi=2AiWi+λ1kT\frac{\partial l}{\partial W_i} = 2A_iW_i+\lambda 1_k^T\\

为求极值点,令偏导数为0,可得:Wi=12λAi11kW_i=-\frac{1}{2}\lambda A^{-1}_i1_k

因为 WiT1k=1W^T_i1_k=1,对WiW_i归一化后上式可以化为:

Wi=Ai11k1kTZi11kW_i=\frac{A^{-1}_i1_k}{1^T_kZ^{-1}_i1_k}

WiW_i为一个向量,若原始数据集中有m个数据,则可求得m个WiW_i,组成一个k×mk\times m的矩阵WW

W=[Wi1...Wm1.........Wik....Wmk]W= \left[ \begin{matrix} W_{i1} & ... & W_{m1}\\ ... & ... & ... \\ W_{ik} & .... & W_{mk} \end{matrix} \right]

求解低维数据

在求得WiW_i后,根据WiW_i可求低维数据,低维数据yiy_i也应满足以下的关系:

yi=wi1yi1+wi2yi2+...+wikyiky_i = w_{i1} * y_{i1} + w_{i2} * y_{i2} + ...+ w_{ik} * y_{ik}

同样,通过最小化损失函数(均方差)来求解,高维的式子中高维数据已知因此求WW,低维则反过来在已知WW的情况下求低维数据。

l(y)=i=1myij=1kwijyij2=i=1m(yIiymi)2l(y)=\sum^m_{i=1}\lVert y_i - \sum^k_{j=1}w_{ij}y_{ij}\rVert^2\\ = \sum^m_{i=1}\lVert (yI_i-ym_i) \rVert^2

其中,IiI_i表示从低维数据集y中取的第i个数据,mim_i则是由wikw_ik组成的稀疏矩阵。进一步化简得:

l(y)=i=1my(Iimi)2=yRRTyTl(y)=\sum^m_{i=1}\lVert y(I_i-m_i)\rVert^2\\ =yRR^Ty^T

其中,R=(Iimi)R=(I_i-m_i)

为了得到标准化的低维数据,添加两个约束:

1mi=1myiyiT=Idi=1myi=0\frac{1}{m}\sum^m_{i=1}y_iy_i^T=I_d\\ \sum^m_{i=1}y_i=0

其中,IdI_d为d维的单位矩阵。第一个约束表示低维数据的协方差矩阵之和为单位矩阵的m倍,第二个约束表示所有低维数据的和为0。

第一个约束条件可化为:

yyT=nIdyy^T=nI_d

其中,y=[y1,y2,...,ym]y=[y_1,y_2,...,y_m]

同样使用拉格朗日乘子法得:

l(y)=yRRTyT+λ(yyTnId)l(y)=yRR^Ty^T+\lambda(yy^T-nI_d)

对y求偏导,令偏导数为0得:

ly=2RRTyT+λyTRRTyT=12λyT\frac{\partial l}{\partial y} = 2RR^Ty^T+\lambda y^T\\ RR^Ty^T=-\frac{1}{2}\lambda y^T

σ=12λ,M=RRT\sigma=-\frac{1}{2}\lambda,M=RR^T,上式化为:

MyT=σyTMy^T=\sigma y^T

由上文的特征向量和特征值公式可得,yTy^T为矩阵MM的特征向量构成的矩阵。要得到最小d维数据集,只需求出矩阵MM最小的d个特征向量组成的矩阵Y=(y1,y2,...,yd)TY=(y_1,y_2,...,y_d)^T

拉普拉斯特征映射(Laplacian Eigen-map, LE)

LE方法的基本思想和LLE类似,但区别在于其使用了拉普拉斯矩阵。

拉普拉斯矩阵

拉普拉斯矩阵L的定义为:L=DWL=D-W。其中D为度矩阵,W为邻接矩阵。

度矩阵和邻接矩阵

度矩阵和邻接矩阵的概念与图有关,对于图 G(V,E)G(V,E),,V=(v1,v2,...,vn)V=(v_1,v_2,...,v_n)是图中所有点的集合,E是边的集合,wijw_{ij}是点viv_ivjv_j间的权重,没有边连接的两个点间权重为0。

对于图中任意一点viv_i,它的度did_i定义为和该点相连的所有边的权重和:

di=j=1nwijd_i=\sum^n_{j=1}w_{ij}

根据这点,可以得到一个n×nn\times n的度矩阵DD,为一个对角矩阵,如下所示:

[d1d2dn]\left[ \begin{matrix} d_1 & \cdots & \cdots& \cdots\\ \cdots & d_2 & \cdots& \cdots\\ \vdots & \vdots &\ddots & \vdots\\ \cdots & \cdots & \cdots & d_n \end{matrix} \right]

此外,还可以得到图的n×nn\times n邻接矩阵WW,第i行的第j个值表示对应wijw_{ij},可以看出DDWW都是对称矩阵。

拉普拉斯矩阵的性质

  • 拉普拉斯矩阵是对称矩阵(因为DDWW都是对称矩阵)
  • 所有特征值都是实数。(对称矩阵的性质)
  • 通过其定义易证明,对任意向量ff有:
fTLf=12i,j=1nwij(fifj)2f^TLf=\frac{1}{2}\sum^n_{i,j=1}w_{ij}(f_i-f_j)^2
  • 为半正定矩阵,且有0=λ1λ2...λn0=\lambda_1\leq\lambda_2\leq...\lambda_n,最小的特征值为0。

LE的思想

与LLE类似,LE首先在高维空间中描述样本结构。假设原始数据集为X=[x1,x2,...,xn],XRmX=[x_1,x_2,...,x_n],X\in R^m,为一个n×mn\times m的集合。

对于样本点xix_i,若存在点xjx_j为其近邻则有:

Wij=expxixj2/tW_{ij}=exp^{\lVert x_i-x_j\rVert ^2 /t}

t为指定的某个常数,若xjx_j不是近邻则Wij=0W_{ij}=0。(这里是使用K-NN算法构建邻接矩阵)

通过高维计算出WW后,在低维空间中需要保持样本的结构,低维数据Y=[y1,y2,...,yn],YRkY=[y_1,y_2,...,y_n],Y \in R^k,为一个n×kn\times k的数据集。在低维空间中要让xix_i尽可能接近xjx_j,则需要最小化目标函数:

mini,jyiyj2Wij=i=1nj=1n(yiTyi2yiTyj+yjTyj)Wij=2i=1nDiiyiTyi2i=1nj=1nyiTyiWijmin \sum_{i,j}\lVert y_i-y_j\rVert^2 W_{ij}\\ = \sum^n_{i=1}\sum^n_{j=1}(y_i^Ty_i-2y_i^Ty_j+y_j^Ty_j)W_{ij}\\ =2\sum^n_{i=1}D_{ii}y_i^Ty_i-2\sum^n_{i=1}\sum^n_{j=1}y^T_iy_iW_{ij}

此处的Dii=j=1nWijD_{ii}=\sum^n_{j=1}W_{ij},即图的度矩阵。令Y=i=1nyiY=\sum^n_{i=1}y_i,上式可继续化为:

2trace(YTDY)2i=1nyiT(j=1njiWij)=2trace(YTDY)2i=1nyiT(YW)i=2trace(YTDY)2trace(TTWY)=2trace(YT(DW)Y)2trace(Y^TDY)-2\sum^n_{i=1}y^T_i(\sum^n_{j=1}j_iW_{ij})\\ = 2trace(Y^TDY) - 2\sum^n_{i=1}y^T_i(YW)_i\\ = 2trace(Y^TDY) - 2trace(T^TWY)\\ = 2trace(Y^T(D-W)Y)

由拉普拉斯矩阵的定义L=DWL=D-W,可得最终转化的目标函数:

mintrace(YTLY)s.t.YTDY=Imin trace(Y^TLY) s.t. Y^TDY=I

YTDY=IY^TDY=I的限制条件是为了消除低维空间中的缩放因子,并保证优化问题有解。

同样使用拉格朗日乘子法求解,其中Λ\Lambda为一个对角矩阵:

f(Y)=trace(YTLY)+trace(Λ(YTDYI))f(Y)=trace(Y^TLY)+trace(\Lambda(Y^TDY-I))

对Y求偏导,用到了矩阵的迹求导

f(Y)Y=LY+LTY+DTYΛT+DYΛ=2LY+2DYΛ=0\frac{\partial f(Y)}{\partial Y} = LY+L^TY+D^TY\Lambda^T+DY\Lambda\\ = 2LY+2DY\Lambda = 0

即:LY=DYΛLY=-DY\Lambda,对于单个样本值yy,式子可化为D1Ly=λyD^{-1}Ly=\lambda y,可以看出这变成了广义特征值分解问题,通过对D1LD^{-1}L进行特征值分解取其第二小到第m小(最小值为0,故舍弃)的特征值所对应的特征向量即为降维后的YY