主成分分析(Principal components analysis,以下简称PCA)是最重要的降维方法之一。在数据压缩消除冗余和数据噪音消除等领域都有广泛的应用。
在理解特征提取与处理时,涉及高维特征向量的问题往往容易陷入维度灾难。随着数据集维度的增加,算法学习需要的样本数量呈指数级增加。有些应用中,遇到这样的大数据是非常不利的,而且从大数据集中学习需要更多的内存和处理能力。另外,随着维度的增加,数据的稀疏性会越来越高。在高维向量空间中探索同样的数据集比在同样稀疏的数据集中探索更加困难。
主成分分析也称为卡尔胡宁-勒夫变换(Karhunen-Loeve Transform),是一种用于探索高维数据结构的技术。PCA通常用于高维数据集的探索与可视化。还可以用于数据压缩,数据预处理等。PCA可以把可能具有相关性的高维变量合成线性无关的低维变量,称为主成分( principal components)。新的低维数据集会尽可能的保留原始数据的变量。
假如你是一本养花工具宣传册的摄影师,你正在拍摄一个水壶。水壶是三维的,但是照片是二维的,为了更全面的把水壶展示给客户,你需要从不同角度拍几张图片。下图是你从四个方向拍的照片
第一张图里水壶的背面可以看到,但是看不到前面。第二张图是拍前面,可以看到壶嘴,这张图可以提供了第一张图缺失的信息,但是壶把看不到了。从第三张俯视图里无法看出壶的高度。第四张图是你真正想要的,水壶的高度,顶部,壶嘴和壶把都清晰可见。
首先考虑一个问题:对于正交属性空间中的样本点,如何用一个超平面(直线的高维推广)对所有样本进行恰当的表达?
可以想到,若存在这样的超平面,那么它大概具有这样的性质:
- 最近重构性:样本点到这个超平面的距离足够近
- 最大可分性:样本点在这个超平面上的投影能尽可能的分开
基于最近重构性和最大可分性能分别得到主成分分析的两种等价推到,我们这里主要考虑最大可分性,并且一步一步推到出最终PCA。
1.PCA最大可分性的思想
PCA顾名思义,就是找出数据里最主要的方面,用数据里最主要的方面来代替原始数据。具体的,假如我们的数据集是 n n n维的,共有 m m m个数据 ( x 1 , x 2 , … , x m ) \left(x_{1}, x_{2}, \ldots, x_{m}\right) (x1,x2,…,xm) 。我们希望将这 m m m个数据的维度从 n n n维降到 n ′ n^{\prime} n′ 维,希望这 m m m 个 n ′ n^{\prime} n′ 维的数据集尽可能的代表原始数据集。我们知道数据从 n n n维降到 n ′ n^{\prime} n′维肯定会有损失,但是我们希望损失尽可能的小。那么如何让这 n ′ n^{\prime} n′维的数据尽可能表示原来的数据呢?
我们先看看最简单的情况,也就是
n
=
2
n = 2
n=2,
n
′
=
1
n^{\prime} = 1
n′=1,也就是将数据从二维降维到一维。数据如下图。我们希望找到某一个维度方向,它可以代表这两个维度的数据。图中列了两个向量方向,
u
1
u_{1}
u1 和
u
2
u_{2}
u2,那么哪个向量可以更好的代表原始数据集呢?
从直观上也可以看出,
u
1
u_{1}
u1 比
u
2
u_{2}
u2 好,这就是我们所说的最大可分性。
2.基变换
一般来说,欲获得原始数据新的表示空间,最简单的是对原始数据进行线性变换(基变换):
Y = P X Y = PX Y=PX
其中
Y
Y
Y 是样本在新空间的表达,
P
P
P是基向量,
X
X
X是原始样本。我们可知选择不同的基可以对一组数据给出不同的表示,同时当基的数量少于原始样本本身的维数则可达到降维的效果,矩阵表示如下:
(
p
1
p
2
⋮
p
R
)
(
a
1
a
2
⋯
a
M
)
=
(
p
1
a
1
p
1
a
2
⋯
p
1
a
M
p
2
a
1
p
2
a
2
⋯
p
2
a
M
⋮
⋮
⋱
⋮
p
R
a
1
p
R
a
2
⋯
p
R
a
M
)
\left(\begin{array}{c}p_{1} \\ p_{2} \\ \vdots \\ p_{R}\end{array}\right)\left(\begin{array}{cccc}a_{1} & a_{2} & \cdots & a_{M}\end{array}\right)=\left(\begin{array}{cccc}p_{1} a_{1} & p_{1} a_{2} & \cdots & p_{1} a_{M} \\ p_{2} a_{1} & p_{2} a_{2} & \cdots & p_{2} a_{M} \\ \vdots & \vdots & \ddots & \vdots \\ p_{R} a_{1} & p_{R} a_{2} & \cdots & p_{R} a_{M}\end{array}\right)
⎝⎜⎜⎜⎛p1p2⋮pR⎠⎟⎟⎟⎞(a1a2⋯aM)=⎝⎜⎜⎜⎛p1a1p2a1⋮pRa1p1a2p2a2⋮pRa2⋯⋯⋱⋯p1aMp2aM⋮pRaM⎠⎟⎟⎟⎞
其中 p i ∈ { p 1 , p 2 , … , p R } , p i ∈ R 1 ∗ N p_{i} \in\left\{p_{1}, p_{2}, \ldots, p_{R}\right\}, \quad p_{i} \in \mathbb{R}^{1 * N} pi∈{p1,p2,…,pR},pi∈R1∗N 是一个行向量,表示第 i i i个基 a j ∈ { a 1 , a 2 , … a M } , a i ∈ R N ∗ 1 a_{j} \in\left\{a_{1}, a_{2}, \ldots a_{M}\right\}, \quad a_{i} \in \mathbb{R}^{N * 1} aj∈{a1,a2,…aM},ai∈RN∗1是一个列向量,表示第 j j j个原始数据记录。特别要注意的是,这里 R R R可以小于 N N N,而 R R R 决定了变换后数据的维数。也就是说,我们可以将一个 N N N维数据变换到更低维度的空间中去,变换后的维度取决于基的数量。从原本 X ∈ R N ∗ M X \in \mathbb{R}^{N * M} X∈RN∗M降维到 Y ∈ R R ∗ M Y \in \mathbb{R}^{R * M} Y∈RR∗M 。因此这种矩阵相乘的表示也可以表示降维变换。
最后,上述分析同时给矩阵相乘找到了一种物理解释:两个矩阵相乘的意义是将右边矩阵中的每一列列向量变换到左边矩阵中每一行行向量为基所表示的空间中去。更抽象的说,一个矩阵可以表示一种线性变换。很多同学在学线性代数时对矩阵相乘的方法感到奇怪,但是如果明白了矩阵相乘的物理意义,其合理性就一目了然了。
3.方差
那么考虑,如何选择一个方向或 者基才是最优的?观察下图
我们将所有的点分别向两条直线做投影,基于前面PCA最大可分思想,我们要找的方向是降维后损失最小,可以理解为投影后的数据尽可能的分开,那么这种分散程度可以用数学上的方差来表示,方差越大数据越分散。方差公式如下:
Var ( a ) = 1 m ∑ i = 1 m ( a i − μ ) 2 \operatorname{Var}(a)=\frac{1}{m} \sum_{i=1}^{m}\left(a_{i}-\mu\right)^{2} Var(a)=m1∑i=1m(ai−μ)2
对数据进行了中心化后(可以方便后面的操作):
Var ( a ) = 1 m ∑ i = 1 m a i 2 \operatorname{Var}(a)=\frac{1}{m} \sum_{i=1}^{m} a_{i}^{2} Var(a)=m1∑i=1mai2
现在我们已经知道了以下几点:
- 对原始样本进行(线性变换)基变换可以对原始样本给出不同的表示
- 基的维度小于数据的维度可以起到降维的效果
- 对基变换后新的样本求其方差,选取使其方差最大的基
- 那么在下面我们来考虑一个新的问题
上面我们导出了优化目标,但是这个目标似乎不能直接作为操作指南(或者说算法),因为它只说要什么,但根本没有说怎么做。所以我们要继续在数学上研究计算方案。
4.协方差
从二维降到一维可以使用方差最大来选出能使基变换后数据分散最大的方向(基),但如果遇到高维的变换,当完成第一个方向(基)选择后,第二个投影方向应该与第一个“几乎重合在一起”,这显然是没有用的,因此要有其它的约束条件。我们希望两个字段尽可能表示更多的信息,使其不存在相关性。
数学上用协方差表示其相关性:
Cov ( a , b ) = 1 m ∑ i = 1 m a i b i \operatorname{Cov}(a, b)=\frac{1}{m} \sum_{i=1}^{m} a_{i} b_{i} Cov(a,b)=m1∑i=1maibi
当 Cov ( a , b ) = 0 \operatorname{Cov}(a, b)=0 Cov(a,b)=0时,表示两个字段完全独立,这也是我们的优化目标。
5.协方差矩阵
我们想达到的目标与字段内方差及字段间协方差有密切关系,假如只有 a a a 、 b b b 两个字段,那么我们将它们按行组成矩阵 X X X ,表示如下:
X = ( a 1 a 2 … a m b 1 b 2 … b m ) X=\left(\begin{array}{llll}a_{1} & a_{2} & \dots & a_{m} \\ b_{1} & b_{2} & \dots & b_{m}\end{array}\right) X=(a1b1a2b2……ambm)
然后我们用 X X X 乘以 X X X 的转置,并乘上系数 1 m \frac{1}{m} m1 :
1 m X X ⊤ = ( 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} X X^{\top}=\left(\begin{array}{cc}\frac{1}{m} \sum_{i=1}^{m} a_{i}^{2} & \frac{1}{m} \sum_{i=1}^{m} a_{i} b_{i} \\ \frac{1}{m} \sum_{i=1}^{m} a_{i} b_{i} & \frac{1}{m} \sum_{i=1}^{m} b_{i}^{2}\end{array}\right) m1XX⊤=(m1∑i=1mai2m1∑i=1maibim1∑i=1maibim1∑i=1mbi2)
可见,协方差矩阵是一个对称的矩阵,而且对角线是各个维度的方差,而其它元素是 [公式] 和 [公式] 的协方差,然后会发现两者被统一到了一个矩阵的。
6.协方差矩阵对角化
我们的目标是使 1 m ∑ i = 1 m a i b i = 0 \frac{1}{m} \sum_{i=1}^{m} a_{i} b_{i}=0 m1∑i=1maibi=0,根据上述推倒,可以看出我们的优化目标 C = 1 m X X ⊤ C=\frac{1}{m} X X^{\top} C=m1XX⊤等价于协方差矩阵对角化。即除对角线外的其它元素( 如 1 m ∑ i = 1 m a i b i \frac{1}{m} \sum_{i=1}^{m} a_{i} b_{i} m1∑i=1maibi )化为0,并且在对角线上将元素按大小从上到下排列,这样我们就达到了优化目的。这样说可能还不是很明晰,我们进一步看下原矩阵与基变换后矩阵协方差矩阵的关系:
设原始数据矩阵 X X X 对应的协方差矩阵为 C C C ,而 P P 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
⊤
D=\frac{1}{m} Y Y^{\top}
D=m1YY⊤
=
1
m
(
P
X
)
(
P
X
)
⊤
=\frac{1}{m}(P X)(P X)^{\top}
=m1(PX)(PX)⊤
=
1
m
P
X
X
⊤
P
⊤
=\frac{1}{m} P X X^{\top} P^{\top}
=m1PXX⊤P⊤
=
P
(
1
m
X
X
⊤
)
P
⊤
=P\left(\frac{1}{m} X X^{\top}\right) P^{\top}
=P(m1XX⊤)P⊤
=
P
C
P
⊤
=P C P^{\top}
=PCP⊤
=
P
(
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
)
P
⊤
=P\left(\begin{array}{cc}\frac{1}{m} \sum_{i=1}^{m} a_{i}^{2} & \frac{1}{m} \sum_{i=1}^{m} a_{i} b_{i} \\ \frac{1}{m} \sum_{i=1}^{m} a_{i} b_{i} & \frac{1}{m} \sum_{i=1}^{m} b_{i}^{2}\end{array}\right) P^{\top}
=P(m1∑i=1mai2m1∑i=1maibim1∑i=1maibim1∑i=1mbi2)P⊤
成了寻找一个矩阵 P , P, P, 满足 P C P ⊤ P C P^{\top} PCP⊤ 是一个对角矩阵, 并且对角元素按从大到小依次排列, 那么 P P P 的前 K K K 行就是要寻找的基, 用 P P P 的前 K K K 行组成的矩阵乘以 X X X 就使得 X X X 从 N N N 维降 到了 K 维并满足上述优化条件。
我们希望的是投影后的方差最大化,于是我们的优化目标可以写为:
max
P
tr
(
P
C
P
⊤
)
\max _{P} \operatorname{tr}\left(P C P^{\top}\right)
maxPtr(PCP⊤)
s.t.
P
P
⊤
=
I
P P^{\top}=I
PP⊤=I
利用拉格朗日函数可以得到:
J ( P ) = tr ( P C P ⊤ ) + λ ( P P ⊤ − I ) J(P)=\operatorname{tr}\left(P C P^{\top}\right)+\lambda\left(P P^{\top}-I\right) J(P)=tr(PCP⊤)+λ(PP⊤−I)
对 P P P 求导有 C P ⊤ + λ P ⊤ = 0 C P^{\top}+\lambda P^{\top}=0 CP⊤+λP⊤=0, 整理下即为:
C P ⊤ = ( − λ ) P ⊤ C P^{\top}=(-\lambda) P^{\top} CP⊤=(−λ)P⊤
于是, 只需对协方差矩阵 C 进行特征分解,对求得的特征值进行排序, 再对
P
⊤
=
(
P
1
,
P
2
,
…
,
P
R
)
P^{\top}=\left(P_{1}, P_{2}, \ldots, P_{R}\right) \quad
P⊤=(P1,P2,…,PR) 取前
K
K
K 列组成的矩阵乘以原始数据矩阵X, 就得到了我们需要
的降维后的数据矩阵Y。
7.PCA算法流程
从上面两节我们可以看出, 求样本
x
i
x_{i}
xi 的
n
′
n^{\prime}
n′ 维的主成分其实就是求样本集的协方差矩阵
1
m
X
X
⊤
\frac{1}{m} X X^{\top}
m1XX⊤ 的前
n
′
n^{\prime}
n′ 个特征值对应特征向量矩阵
P
,
P,
P, 然后对于每个样本
x
i
x_{i}
xi, 做如下变换
y
i
=
P
x
i
,
y_{i}=P x_{i},
yi=Pxi, 即达到降维的PCA目的。
下面我们看看具体的算法流程:
输入:
n
n
n 维样本集
X
=
(
x
1
,
x
2
,
…
,
x
m
)
,
X=\left(x_{1}, x_{2}, \ldots, x_{m}\right),
X=(x1,x2,…,xm), 要降维到的维数
n
′
n^{\prime}
n′
输出:降维后的样本集
Y
Y
Y
- 1.对所有的样本进行中心化 x i = x i − 1 m ∑ j = 1 m x j x_{i}=x_{i}-\frac{1}{m} \sum_{j=1}^{m} x_{j} xi=xi−m1∑j=1mxj
- 2.计算样本的协方差矩阵 C = 1 m X X ⊤ C=\frac{1}{m} X X^{\top} C=m1XX⊤
- 3.求出协方差矩阵的特征值及对应的特征向量
- 4.将特征向量按对应特征值大小从上到下按行排列成矩阵,取前k行组成矩阵P
- 5.Y=PX即为降维到k维后的数据
注意:
有时候, 我们不指定降维后的
n
′
n^{\prime}
n′ 的值,而是换种方式, 指定一个降维到的主成分比重阅值
t
t
t 这个闰值t在 (0,1] 之间。假如我们的
n
n
n 个特征值为
λ
1
≥
λ
2
≥
…
≥
λ
n
\lambda_{1} \geq \lambda_{2} \geq \ldots \geq \lambda_{n}
λ1≥λ2≥…≥λn,则n’可以通过下 式得到:
∑
i
=
1
n
′
λ
i
∑
i
=
1
n
λ
i
≥
t
\frac{\sum_{i=1}^{n^{\prime}} \lambda_{i}}{\sum_{i=1}^{n} \lambda_{i}} \geq t
∑i=1nλi∑i=1n′λi≥t
8.PCA算法总结
这里对PCA算法做一个总结。作为一个非监督学习的降维方法,它只需要特征值分解,就可以对数据进行压缩,去噪。因此在实际场景应用很广泛。为了克服PCA的一些缺点,出现了很多PCA的变种,比如为解决非线性降维的KPCA,还有解决内存限制的增量PCA方法Incremental PCA,以及解决稀疏数据降维的PCA方法Sparse PCA等。
PCA算法的主要优点有:
- 仅仅需要以方差衡量信息量,不受数据集以外的因素影响。
- 各主成分之间正交,可消除原始数据成分间的相互影响的因素。
- 计算方法简单,主要运算是特征值分解,易于实现。
PCA算法的主要缺点有:
- 主成分各个特征维度的含义具有一定的模糊性,不如原始样本特征的解释性强。
- 方差小的非主成分也可能含有对样本差异的重要信息,因降维丢弃可能对后续数据处理有影响。