PCA(Principal Component Analysis)是一种常用的数据分析方法。PCA通过线性变换将原始数据变换为一组各维度线性无关的表示,可用于提取数据的主要特征分量,常用于高维数据的降维。
数据的向量表示及降维问题
一般情况下,在数据挖掘和机器学习中,数据被表示为向量。例如某个淘宝店2012年全年的流量及交易情况可以看成一组记录的集合,其中每一天的数据是一条记录,格式如下:
(日期, 浏览量, 访客数, 下单数, 成交数, 成交金额)
其中“日期”是一个记录标志而非度量值,而数据挖掘关心的大多是度量值,因此如果我们忽略日期这个字段后,我们得到一组记录,每条记录可以被表示为一个五维向量,其中一条看起来大约是这个样子:
( 500 , 240 , 25 , 13 , 2312.15 ) T (500,240,25,13,2312.15)^T (500,240,25,13,2312.15)T
我们当然可以对这一组五维向量进行分析和挖掘,不过我们知道,很多机器学习算法的复杂度和数据的维数有着密切关系,甚至与维数呈指数级关联。当然,这里区区五维的数据,也许还无所谓,但是实际机器学习中处理成千上万甚至几十万维的情况也并不罕见,在这种情况下,机器学习的资源消耗是不可接受的,因此我们必须对数据进行降维。
降维当然意味着信息的丢失,不过鉴于实际数据本身常常存在的相关性,我们可以想办法在降维的同时将信息的损失尽量降低。
举个例子,假如某学籍数据有两列M和F,其中M列的取值是如何此学生为男性取值1,为女性取值0;而F列是学生为女性取值1,男性取值0。此时如果我们统计全部学籍数据,会发现对于任何一条记录来说,当M为1时F必定为0,反之当M为0时F必定为1。在这种情况下,我们将M或F去掉实际上没有任何信息的损失,因为只要保留一列就可以完全还原另一列。
当然上面是一个极端的情况,在现实中也许不会出现,不过类似的情况还是很常见的。例如上面淘宝店铺的数据,从经验我们可以知道,“浏览量”和“访客数”往往具有较强的相关关系,而“下单数”和“成交数”也具有较强的相关关系。这里我们非正式的使用“相关关系”这个词,可以直观理解为“当某一天这个店铺的浏览量较高(或较低)时,我们应该很大程度上认为这天的访客数也较高(或较低)”。后面的章节中我们会给出相关性的严格数学定义。
这种情况表明,如果我们删除浏览量或访客数其中一个指标,我们应该期待并不会丢失太多信息。因此我们可以删除一个,以降低机器学习算法的复杂度。
上面给出的是降维的朴素思想描述,可以有助于直观理解降维的动机和可行性,但并不具有操作指导意义。例如,我们到底删除哪一列损失的信息才最小?亦或根本不是单纯删除几列,而是通过某些变换将原始数据变为更少的列但又使得丢失的信息最小?到底如何度量丢失信息的多少?如何根据原始数据决定具体的降维操作步骤?
向量的表述
要准确描述向量,首先要确定一组基,然后给出在基所在的各个直线上的投影值,就可以了。只不过我们经常省略第一步,而默认以(1,0)和(0,1)为基,也就是自然基。
一个向量在一组单位正交基上的投影,也就是找到向量在对应单位正交基上的坐标。
上述(4)式表明 p r o j w y proj_wy projwy是 U U U的列的线性组合,且对应权值分别为 y ⋅ u 1 ⋯ y ⋅ u p y\cdot u_1 \cdots y\cdot u_p y⋅u1⋯y⋅up,是个标量,也就是在对应正交基下的坐标值。
举个例子,自然基下的向量 ( 3 , 2 ) T (3,2)^T (3,2)T,同时规定了一组新的单位正交基,
( 1 2 1 2 ) , ( − 1 2 1 2 ) \begin{pmatrix} \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}} \end{pmatrix}, \begin{pmatrix} \frac{-1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}} \end{pmatrix} (2121),(2−121)
那么利用上面的投影公式,可以求出在新基下的坐标,用矩阵表示如下,
( 1 2 1 2 − 1 2 1 2 ) ( 3 2 ) = ( 5 2 − 1 2 ) \begin{pmatrix} \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}}\\ \frac{-1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{pmatrix} \begin{pmatrix} 3\\ 2 \end{pmatrix}= \begin{pmatrix} \frac{5}{\sqrt{2}}\\ \frac{-1}{\sqrt{2}} \end{pmatrix} (212−12121)(32)=(252−1)
协方差矩阵及优化目标
如果我们有一组N维向量,现在要将其降到K维(K小于N),那么我们应该如何选择K个基才能最大程度保留原有的信息?
假设我们的数据由五条记录组成,将它们表示成矩阵形式:
( 1 1 2 4 2 1 3 3 4 4 ) \begin{pmatrix} 1&1&2&4&2\\ 1&3&3&4&4 \end{pmatrix} (1113234424)
其中每一列为一条数据记录,而一行为一个字段。为了后续处理方便,我们首先将每个字段内所有值都减去字段均值,其结果是将每个字段都变为均值为0(这样做的道理和好处后面会看到)。
我们看上面的数据,第一个字段均值为2,第二个字段均值为3,所以变换后:
( − 1 − 1 0 2 0 − 2 0 0 1 1 ) \begin{pmatrix} -1&-1&0&2&0\\ -2&0&0&1&1 \end{pmatrix} (−1−2−10002101)
我们可以看下五条数据在平面直角坐标系内的样子:
现在问题来了:如果我们必须使用一维来表示这些数据,又希望尽量保留原始的信息,你要如何选择?
通过上一节对基变换的讨论我们知道,这个问题实际上是要在二维平面中选择一个方向,将所有数据都投影到这个方向所在直线上,用投影值表示原始记录。这是一个实际的二维降到一维的问题。
那么如何选择这个方向(或者说基)才能尽量保留最多的原始信息呢?一种直观的看法是:希望投影后的投影值尽可能分散。
以上图为例,可以看出如果向x轴投影,那么最左边的两个点会重叠在一起,中间的两个点也会重叠在一起,于是本身四个各不相同的二维点投影后只剩下两个不同的值了,这是一种严重的信息丢失,同理,如果向y轴投影最上面的两个点和分布在x轴上的两个点也会重叠。所以看来x和y轴都不是最好的投影选择。我们直观目测,如果向通过第一象限和第三象限的斜线投影,则五个点在投影后还是可以区分的。
下面,我们用数学方法表述这个问题。
方差
上文说到,我们希望投影后投影值尽可能分散,而这种分散程度,可以用数学上的方差来表述。此处,一个字段的方差可以看做是每个元素与字段均值的差的平方和的均值,即:
V a r ( a ) = 1 m ∑ i = 1 m ( a i − μ ) 2 Var(a)=\frac{1}{m}\sum_{i=1}^m(a_i-\mu)^2 Var(a)=m1i=1∑m(ai−μ)2
于上面我们已经将每个字段的均值都化为0了,因此方差可以直接用每个元素的平方和除以元素个数表示:
V a r ( a ) = 1 m ∑ i = 1 m ( a i ) 2 Var(a)=\frac{1}{m}\sum_{i=1}^m(a_i)^2 Var(a)=m1i=1∑m(ai)2
于是上面的问题被形式化表述为:寻找一个一维基,使得所有数据变换为这个基上的坐标表示后,方差值最大。
协方差
对于上面二维降成一维的问题来说,找到那个使得方差最大的方向就可以了。不过对于更高维,还有一个问题需要解决。考虑三维降到二维问题。与之前相同,首先我们希望找到一个方向使得投影后方差最大,这样就完成了第一个方向的选择,继而我们选择第二个投影方向。
如果我们还是单纯只选择方差最大的方向,很明显,这个方向与第一个方向应该是“几乎重合在一起”,显然这样的维度是没有用的,因此,应该有其他约束条件。从直观上说,让两个字段尽可能表示更多的原始信息,我们是不希望它们之间存在(线性)相关性的,因为相关性意味着两个字段不是完全独立,必然存在重复表示的信息。
数学上可以用两个字段的协方差表示其相关性,由于已经让每个字段均值为0,则:
C o v ( a , b ) = 1 m ∑ i = 1 m a i b i Cov(a,b)=\frac{1}{m}\sum_{i=1}^m a_ib_i Cov(a,b)=m1i=1∑maibi
可以看到,在字段均值为0的情况下,两个字段的协方差简洁的表示为其内积除以元素数m。
当协方差为0时,表示两个字段完全独立。为了让协方差为0,我们选择第二个基时只能在与第一个基正交的方向上选择。因此最终选择的两个方向一定是正交的。
至此,我们得到了降维问题的优化目标:将一组N维向量降为K维(K大于0,小于N),其目标是选择K个单位(模为1)正交基,使得原始数据变换到这组基上后,各字段两两间协方差为0,而字段的方差则尽可能大(在正交的约束下,取最大的K个方差)。
协方差矩阵
上面我们导出了优化目标,但是这个目标似乎不能直接作为操作指南(或者说算法),因为它只说要什么,但根本没有说怎么做。所以我们要继续在数学上研究计算方案。
我们看到,最终要达到的目的与字段内方差及字段间协方差有密切关系。因此我们希望能将两者统一表示,仔细观察发现,两者均可以表示为内积的形式,而内积又与矩阵相乘密切相关。于是我们来了灵感:
假设我们只有a和b两个字段,那么我们将它们按行组成矩阵 X X X:
X = ( a 1 a 2 … a m b 1 b 2 … b m ) X= \begin{pmatrix} a_1 & a_2 & \dots & a_m \\ b_1 & b_2 & \dots & b_m \end{pmatrix} X=(a1b1a2b2……ambm)
然后我们用 X X X乘以 X X X的转置,并乘上系数 1 m \frac{1}{m} m1:
1 m X X T = ( 1 m ∑ i = 1 m a i 2 1 m ∑ i = 1 m a i b i 1 m ∑ i = 1 m a i b i 1 m ∑ i = 1 m b i 2 ) \frac{1}{m}XX^T= \begin{pmatrix} \frac{1}{m}\sum_{i=1}^m a_i^2 & \frac{1}{m}\sum_{i=1}^m a_ib_i \\ \frac{1}{m}\sum_{i=1}^m a_ib_i & \frac{1}{m}\sum_{i=1}^m b_i^2 \end{pmatrix} m1XXT=(m1∑i=1mai2m1∑i=1maibim1∑i=1maibim1∑i=1mbi2)
奇迹出现了!这个矩阵对角线上的两个元素分别是两个字段的方差,而其它元素是a和b的协方差。两者被统一到了一个矩阵的。
根据矩阵相乘的运算法则,这个结论很容易被推广到一般情况:
设我们有m个n维数据记录,将其按列排成n乘m的矩阵X,设 C = 1 m X X T C=\frac{1}{m}XX^T C=m1XXT,则C是一个对称矩阵,其对角线分别个各个字段的方差,而第i行j列和j行i列元素相同,表示i和j两个字段的协方差。
协方差矩阵对角化
根据上述推导,我们发现要达到优化目的,等价于将协方差矩阵对角化:即除对角线外的其它元素化为0,并且在对角线上将元素按大小从上到下排列,这样我们就达到了优化目的。这样说可能还不是很明晰,我们进一步看下原矩阵与基变换后矩阵协方差矩阵的关系:
设原始数据矩阵X对应的协方差矩阵为C,而P是一组基按行组成的矩阵, Y = P X Y=PX Y=PX,则 Y Y Y为 X X X对 P P P做基变换后的数据。设 Y Y Y的协方差矩阵为 D D D,我们推导一下 D D D与 C C C的关系:
D = 1 m Y Y T = 1 m ( P X ) ( P X ) T = P ( 1 m X X T ) P T = P C P T \begin{aligned} D &= \frac{1}{m}YY^T \\ &= \frac{1}{m}(PX)(PX)^T\\ &= P(\frac{1}{m}XX^T)P^T\\ &= PCP^T \end{aligned} D=m1YYT=m1(PX)(PX)T=P(m1XXT)PT=PCPT
现在事情很明白了!我们要找的 P P P不是别的,而是能让原始协方差矩阵对角化的 P P P。换句话说,优化目标变成了寻找一个矩阵 P P P,满足 P C P T PCP^T PCPT是一个对角矩阵,并且对角元素按从大到小依次排列,那么 P P P的前 K K K行就是要寻找的基,用 P P P的前 K K K行组成的矩阵乘以 X X X就使得 X X X从 N N N维降到了 K K K维并满足上述优化条件。
算法及实例
总结一下PCA的算法步骤:
设有m条n维数据。
- 将原始数据按列组成n行m列矩阵X
- 将X的每一行(代表一个属性字段)进行零均值化,即减去这一行的均值
- 求出协方差矩阵 C = 1 m X X T C=\frac{1}{m}XX^T C=m1XXT
- 求出协方差矩阵的特征值及对应的特征向量
- 将特征向量按对应特征值大小从上到下按行排列成矩阵,取前k行组成矩阵P
- Y = P X Y=PX Y=PX即为降维到k维后的数据
实例
这里以上文提到的
( − 1 − 1 0 2 0 − 2 0 0 1 1 ) \begin{pmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 &1 &1 \end{pmatrix} (−1−2−10002101)
为例,我们用PCA方法将这组二维数据其降到一维。
因为这个矩阵的每行已经是零均值,这里我们直接求协方差矩阵
C = 1 5 ( − 1 − 1 0 2 0 − 2 0 0 1 1 ) ( − 1 − 2 − 1 0 0 0 2 1 0 1 ) = ( 6 5 4 5 4 5 6 5 ) C=\frac{1}{5} \begin{pmatrix} -1 & -1 & 0 & 2 & 0\\ -2 & 0 & 0 & 1 & 1\\ \end{pmatrix} \begin{pmatrix} -1 & -2\\ -1 & 0 \\ 0 & 0\\ 2 & 1\\ 0 & 1 \end{pmatrix} =\begin{pmatrix} \frac{6}{5} & \frac{4}{5} \\ \frac{4}{5} & \frac{6}{5} \end{pmatrix} C=51(−1−2−10002101)⎝⎜⎜⎜⎜⎛−1−1020−20011⎠⎟⎟⎟⎟⎞=(56545456)
然后求其特征值和特征向量,求解后特征值为:
λ
1
=
2
,
λ
2
=
2
/
5
\lambda_1=2 , \lambda_2=2/5
λ1=2,λ2=2/5
其对应的特征向量分别是:
c 1 ( 1 1 ) , c 2 ( − 1 1 ) c1 \begin{pmatrix} 1\\ 1 \end{pmatrix}, c2 \begin{pmatrix} -1\\ 1 \end{pmatrix} c1(11),c2(−11)
其中对应的特征向量分别是一个通解, c 1 c1 c1和 c 2 c2 c2可取任意实数。那么标准化后的特征向量为:
( 1 2 1 2 ) , ( − 1 2 1 2 ) \begin{pmatrix} \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\\ \end{pmatrix}, \begin{pmatrix} \frac{-1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}} \end{pmatrix} (2121),(2−121)
因此我们的矩阵P是:
P = ( 1 2 1 2 − 1 2 1 2 ) P= \begin{pmatrix} \frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\\ \frac{-1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{pmatrix} P=(212−12121)
可以验证协方差矩阵 C C C的对角化:
P C P T = ( 1 2 1 2 − 1 2 1 2 ) ( 6 5 4 5 4 5 6 5 ) ( 1 2 − 1 2 1 2 1 2 ) = ( 2 0 0 2 5 ) PCP^T= \begin{pmatrix} \frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\\ \frac{-1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{pmatrix} \begin{pmatrix} \frac{6}{5} & \frac{4}{5} \\ \frac{4}{5} & \frac{6}{5} \end{pmatrix} \begin{pmatrix} \frac{1}{\sqrt{2}}& \frac{-1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{pmatrix} =\begin{pmatrix} 2 & 0\\ 0 & \frac{2}{5} \end{pmatrix} PCPT=(212−12121)(56545456)(21212−121)=(20052)
最后我们用P的第一行乘以数据矩阵,就得到了降维后的表示:
Y = ( 1 2 1 2 ) ( − 1 − 1 0 2 0 − 2 0 0 1 1 ) = ( − 3 2 − 1 2 0 3 2 − 1 2 ) Y= \begin{pmatrix} \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{pmatrix} \begin{pmatrix} -1 & -1 & 0 & 2 & 0\\ -2 & 0 & 0 & 1 & 1\\ \end{pmatrix} =\begin{pmatrix} \frac{-3}{\sqrt{2}} & \frac{-1}{\sqrt{2}} & 0 & \frac{3}{\sqrt{2}} & \frac{-1}{\sqrt{2}} \end{pmatrix} Y=(2121)(−1−2−10002101)=(2−32−10232−1)
降维投影结果如下图: