01:PCA最大方差理论
在机器学习领域中 ,我们对原始数据进行特征提取,有时会得到 比较高维的特征向量 。 在这些向量所处的高维空间中 , 包含很多的冗余和噪声 。 我们希望通过降维的方式来寻找数据内部的特性 , 从而提升特征表达能力 , 降低训练复杂度 。 主成分分析( Principal Components Analys is, P CA )作为降维中最经典的方法,至今已有 100 多年的历史,属于一种线性、非监督、全局的降维算法,是面试中经常被问到的问题 。
知识点:PCA ,线性代数
问题:如何定义主成分?从这种定义出发,如何设计目标函数使得降维达到提取主成份的目的?针对这个目标函数,如何对 PCA 问题进行求解?
分析与解答:
PCA 旨在找到数据中的主成分,并利用这些主成分表征原始数据,从而达到降维的目的。举一个简单的例子,在三维空间中有一系列数据点,这些点分布在一个过原点的平面上。如果我们用自然坐标系
x
,
y
,
z
x, y, z
x,y,z三个轴来表示数据,就需要使用三个维度。而实际上,这些点只出现在一个二维平面上, 如果我们通过坐标系旋转变换使得数据所在平面与
x
,
y
x, y
x,y平面重合,那么我们就可以通过
x
′
,
y
′
x^{\prime}, y^{\prime}
x′,y′ 两个维度表达原始数据,并且没有任何损失,这样就完成了数据的降维。而
x
′
,
y
′
x^{\prime}, y^{\prime}
x′,y′ 两个轴所包含的信息就是我们要找到的主成分。
但在高维空间中,我们往往不能像刚才这样直观地想象出数据的分布形式 ,也就更难精确地找到主成分对应的轴是哪些。不妨,我们先从最简单的二维数据来看看 PCA 究竟是如何工作的,如国 4.1 所示。
图4.1( a ) 是二维空间中经过中心化的一组数据,我们很容易看出主成分所在的轴(以下称为主轴 ) ) ) 的大致方向,即图 4.1 ( b ) 中黄线所处的轴。因为在黄线所处的轴上,数据分布得更为分散,这也意味着数据在这个方向上方差更大。在信号处理领域,我们认为信号具有较大方差,噪声具有较小方差,信号与噪声之比称为信噪比。信噪比越大意味着数据的质量越好,反之,信噪比越小意味着数据的质量越差。由此我们不难引出 PCA 的目标,即最大化投影方差,也就是让数据在主轴上投影的方差最大。
对于给定的一组数据点
{
v
1
,
v
2
,
…
,
v
n
}
,
\left\{v_{1}, v_{2}, \ldots, v_{n}\right\},
{v1,v2,…,vn}, 其中所有向量均为列向量, 中心化后的表示为
{
x
1
,
x
2
,
…
,
x
n
}
=
{
v
1
−
μ
,
v
2
−
μ
,
…
,
v
n
−
μ
}
,
\left\{x_{1}, x_{2}, \ldots, x_{n}\right\}=\left\{v_{1}-\mu, v_{2}-\mu, \ldots, v_{n}-\mu\right\},
{x1,x2,…,xn}={v1−μ,v2−μ,…,vn−μ}, 其中
μ
=
1
n
∑
i
=
1
n
v
i
∘
\mu=\frac{1}{n} \sum_{i=1}^{n} v_{i \circ}
μ=n1∑i=1nvi∘ 我们知道,向量内积在几何上表示为第一个向量投影到第二个向量上的长度,因此向量
x
i
x_{i}
xi 在
ω
(
\omega(
ω( 单位方向向量 ) 上的投影坐标可以表示为
(
x
i
,
ω
)
=
x
i
T
ω
∘
\left(x_{i},\omega\right)=x_{i}^{\mathrm{T}} \omega_{\circ}
(xi,ω)=xiTω∘ 所以目标是找到一个投影方向
ω
,
\omega,
ω, 使得
x
1
,
x
2
,
…
,
x
n
x_{1}, x_{2}, \ldots, x_{n}
x1,x2,…,xn 在
ω
\omega
ω 上的投影方差尽可能大。易知,投影之后均值为
0
\mathbf{0}
0 ( 因为
μ
′
=
1
n
∑
i
=
1
n
x
i
T
ω
=
(
1
n
∑
i
=
1
n
x
i
T
)
ω
=
0
,
\boldsymbol{\mu}^{\prime}=\frac{1}{n} \sum_{i=1}^{n} \boldsymbol{x}_{i}^{\mathrm{T}} \omega=\left(\frac{1}{n} \sum_{i=1}^{n} \boldsymbol{x}_{i}^{\mathrm{T}}\right) \omega=0,
μ′=n1∑i=1nxiTω=(n1∑i=1nxiT)ω=0, 这也是我们进行中心化的意义
)
)
) ,
因此投影后的方差可以表示为
D
(
x
)
=
1
n
∑
i
=
1
n
(
x
i
T
ω
)
2
=
1
n
∑
i
=
1
n
(
x
i
T
ω
)
T
(
x
i
T
ω
)
=
1
n
∑
i
=
1
n
ω
T
x
i
x
i
T
ω
=
ω
T
(
1
n
∑
i
=
1
n
x
i
x
i
T
)
ω
\begin{aligned} D(x)=\frac{1}{n} \sum_{i=1}^{n}\left(x_{i}^{\mathrm{T}} \omega\right)^{2} &=\frac{1}{n} \sum_{i=1}^{n}\left(x_{i}^{\mathrm{T}} \omega\right)^{\mathrm{T}}\left(x_{i}^{\mathrm{T}} \omega\right) \\ &=\frac{1}{n} \sum_{i=1}^{n} \omega^{\mathrm{T}} x_{i} x_{i}^{\mathrm{T}} \omega \\ &=\omega^{\mathrm{T}}\left(\frac{1}{n} \sum_{i=1}^{n} x_{i} x_{i}^{\mathrm{T}}\right) \omega \end{aligned}
D(x)=n1i=1∑n(xiTω)2=n1i=1∑n(xiTω)T(xiTω)=n1i=1∑nωTxixiTω=ωT(n1i=1∑nxixiT)ω
仔细一看,
(
1
n
∑
i
=
1
n
ω
T
x
i
x
i
T
ω
)
\left(\frac{1}{n} \sum_{i=1}^{n} \omega^{\mathrm{T}} x_{i} x_{i}^{\mathrm{T}} \omega\right)
(n1∑i=1nωTxixiTω) 其实就是样本协方差矩阵,我们将其写作
Σ
\Sigma
Σ 。另外,由于
ω
\omega
ω 是单位方向向量,即有
ω
T
ω
=
1
\omega^{\mathrm{T}} \omega=\mathbf{1}
ωTω=1 。因此我们要求解一个最大化问题,可表示为
{
max
{
ω
T
Σ
ω
}
s.t.
ω
T
ω
=
1
\left\{\begin{array}{l} \max \left\{\omega^{\mathrm{T}} \Sigma \omega\right\} \\ \text {s.t. } \quad \omega^{\mathrm{T}} \omega=1 \end{array}\right.
{max{ωTΣω}s.t. ωTω=1
引入拉格朗日乘子, 并对
ω
\omega
ω 求导令其等于
0
,
0,
0, 便可以推出
Σ
ω
=
λ
ω
,
\Sigma \omega=\lambda \omega,
Σω=λω, 此时
D
(
x
)
=
ω
T
Σ
ω
=
λ
ω
T
ω
=
λ
D(x)=\omega^{\mathrm{T}} \Sigma \omega=\lambda \omega^{\mathrm{T}} \omega=\lambda
D(x)=ωTΣω=λωTω=λ
熟悉线性代数可以马上就会发现,原来,
x
x
x 投影后的方差就是协方差矩阵的特征值。我们要找到最大的方差也就是协方差矩阵最大的特征值,最佳投影方向就是最大特征值所对应的特征向量。次佳投影方向位于最佳投影方向的正交空间中,是第二大特征值对应的特征向量,以此类推。至此,我们得到以下几种 PCA 的求解方法。
(1 )对样本数据进行中心化处理。
(2 ) 求样本协方差矩阵。
(3) 对协方差矩阵进行特征值分解,将特征值从大到小排列。
(4 ) 取特征值前
d
d
d 大对应的特征向量
ω
1
,
ω
2
,
…
,
ω
d
,
\omega_{1}, \omega_{2}, \ldots, \omega_{d},
ω1,ω2,…,ωd, 通过以下映射将
n
n
n 维样本映射到
d
d
d 维
x
i
′
=
[
ω
1
T
x
i
ω
2
T
x
i
⋮
ω
d
T
x
i
]
\boldsymbol{x}_{i}^{\prime}=\left[\begin{array}{c} \omega_{1}^{\mathrm{T}} \boldsymbol{x}_{i} \\ \omega_{2}^{\mathrm{T}} \boldsymbol{x}_{i} \\ \vdots \\ \omega_{d}^{\mathrm{T}} \boldsymbol{x}_{i} \end{array}\right]
xi′=⎣⎢⎢⎢⎡ω1Txiω2Txi⋮ωdTxi⎦⎥⎥⎥⎤
新的
x
i
′
x_{i}^{\prime}
xi′ 的第
d
d
d 维就是
x
i
x_{i}
xi 在第
d
d
d 个主成分
ω
d
\omega_{d}
ωd 方向上的投影,通过选取最
大的
d
d
d 个特征值对应的特征向量,我们将方差较小的特征(噪声 ) 抛弃,使
得每个
n
n
n 维列向量
x
i
x_{i}
xi 被映射为
d
d
d 维列向量
x
i
′
,
x_{i}^{\prime},
xi′, 定义降维后的信息占比为
η
=
∑
i
=
1
d
λ
i
2
∑
i
=
1
n
λ
i
2
\eta=\sqrt{\frac{\sum_{i=1}^{d} \lambda_{i}^{2}}{\sum_{i=1}^{n} \lambda_{i}^{2}}}
η=∑i=1nλi2∑i=1dλi2
总结与拓展
至此,我们从最大化投影方差的角度解释了 PCA 的原理、目标函数和求解方法。其实,PCA还可以用其他思路进行分析,比如从最小回归误差的角度得到新的目标函数。但最终我们会发现其对应的原理和求解方法与本文中的是等价的。另外,由于 PCA 是一种线性降维方法,虽然经典,但具有一定的局限性。我们可以通过核映射对 PCA 进行扩展得到核主成分分析 (KPCA), 也可以通过流形映射的降维方法, 比如等距映射、局部线性嵌入、拉普拉斯特征映射等,对一些 PCA 效果不好的复杂数据集进行非线性降维操作。
PCA最小平方误差理论
上一节介绍了从最大方差的角度解释 PCA 的原理 、 目标函数和求解方法。本节将通过最小平方误差的思路对 PCA 进行推导 。
知识点:线性代数,最小平方误差
问题:PCA 求解的真实是最佳投影方向,即一条直线,这与数学中线性回归问题的目标不谋而合,能否从回归的角度定义 PCA 的目标并相应地求解问题呢?
分析与解答:我们还是考虑二维空间中的样本点,如图4.2 所示 。 上一节求解得到一条直线使得样本点投影到该直线上的方差最大。从求解直线的思路出发,很容易联想到数学中的线性回归问题,其目标也是求解一个线性函数使得对应直线能够更好地拟合样本点集合 。 如果我们从这个角度定义PCA 的目标,那么问题就会转化为一个回归问题 。
顺着这个思路,在高维空间中,我们实际上是要找到 一个 d 维超平面,使得数据点到这个超平面的距离平方和最小 。 以 d= 1 为例,超平面退化为直线 , 即把样本点投影到最佳直线,最小化的就是所有点到直线的距离平方之和, 如图 4 .3 所示。
数据集中每个点
x
k
x_{k}
xk 到
d
d
d 维超平面
D
D
D 的距离为
distance
(
x
k
,
D
)
=
∥
x
k
−
x
~
k
∥
2
\operatorname{distance}\left(x_{k}, D\right)=\left\|x_{k}-\widetilde{x}_{k}\right\|_{2}
distance(xk,D)=∥xk−x
k∥2
其中
x
k
~
\widetilde{x_{k}}
xk
表示
x
k
x_{k}
xk 在超平面
D
D
D 上的投影向量。如果该超平面由
d
d
d 个标准正交基
W
=
{
ω
1
,
ω
2
,
…
,
ω
d
}
W=\left\{\omega_{1}, \omega_{2}, \ldots, \omega_{d}\right\}
W={ω1,ω2,…,ωd} 构成,根据线性代数理论
x
k
x_{k}
xk 可以由这组基
线性表示
x
k
~
=
∑
i
=
1
d
(
ω
i
T
x
k
)
ω
i
\widetilde{\boldsymbol{x}_{k}}=\sum_{i=1}^{d}\left(\boldsymbol{\omega}_{i}^{\mathrm{T}} \boldsymbol{x}_{k}\right) \boldsymbol{\omega}_{i}
xk
=i=1∑d(ωiTxk)ωi
其中
ω
i
T
x
k
\omega_{i}^{\mathrm{T}} x_{k}
ωiTxk 表示
x
k
x_{k}
xk 在
ω
i
\omega_{i}
ωi 方向上投影的长度。因此,
x
~
k
\widetilde{x}_{k}
x
k 实际上就是
x
k
x_{k}
xk 在
W
W
W 这组标准正交基下的坐标。而 PCA 要优化的目标为
{
arg
min
ω
1
,
…
,
ω
d
∑
k
=
1
n
∥
x
k
−
x
k
~
∥
2
2
s.t.
ω
i
T
ω
j
=
δ
i
j
=
{
1
,
i
=
j
∀
i
,
j
\left\{\begin{array}{l} \underset{\omega_{1}, \ldots, \omega_{d}}{\arg \min } \sum_{k=1}^{n}\left\|\boldsymbol{x}_{k}-\widetilde{\boldsymbol{x}_{k}}\right\|_{2}^{2} \\ \text {s.t. } \quad \omega_{i}^{\mathrm{T}} \boldsymbol{\omega}_{j}=\delta_{i j}=\left\{\begin{array}{l} 1, i=j \\ \forall i, j \end{array}\right. \end{array}\right.
⎩⎪⎨⎪⎧ω1,…,ωdargmin∑k=1n∥xk−xk
∥22s.t. ωiTωj=δij={1,i=j∀i,j
由向量内积的性质,我们知道
x
k
T
x
k
~
=
x
k
~
T
x
k
,
x_{k}^{\mathrm{T}} \widetilde{x_{k}}=\widetilde{x_{k}}^{\mathrm{T}} x_{k},
xkTxk
=xk
Txk, 于是将式 (4.8) 中的每一个距离展开
∥
x
k
−
x
k
~
∥
2
2
=
(
x
k
−
x
k
~
)
T
(
x
k
−
x
~
k
)
=
x
k
T
x
k
−
x
k
T
x
k
~
−
x
k
~
T
x
k
+
x
k
~
T
x
k
~
=
x
k
T
x
k
−
2
x
k
T
x
k
~
+
x
k
~
T
x
k
~
\begin{aligned}\left\|x_{k}-\widetilde{x_{k}}\right\|_{2}^{2} &=\left(x_{k}-\widetilde{x_{k}}\right)^{\mathrm{T}}\left(x_{k}-\tilde{x}_{k}\right) \\ &=x_{k}^{\mathrm{T}} x_{k}-x_{k}^{\mathrm{T}} \widetilde{x_{k}}-\widetilde{x_{k}}^{\mathrm{T}} x_{k}+\widetilde{x_{k}}^{\mathrm{T}} \widetilde{x_{k}} \\ &=x_{k}^{\mathrm{T}} x_{k}-2 x_{k}^{\mathrm{T}} \widetilde{x_{k}}+\widetilde{x_{k}}^{\mathrm{T}} \widetilde{x_{k}} \end{aligned}
∥xk−xk
∥22=(xk−xk
)T(xk−x~k)=xkTxk−xkTxk
−xk
Txk+xk
Txk
=xkTxk−2xkTxk
+xk
Txk
其中第一项 x k T x k x_{k}^{\mathrm{T}} x_{k} xkTxk 与选取的 W W W 无关, 是个常数。将式 (4.7) 代入式 (4.9)的第二项和第三项可得到:
x
k
T
x
k
~
=
x
k
T
∑
i
=
1
d
(
ω
i
T
x
k
)
ω
i
x_{k}^{\mathrm{T}} \widetilde{x_{k}}=x_{k}^{\mathrm{T}} \sum_{i=1}^{d}\left(\omega_{i}^{\mathrm{T}} x_{k}\right) \omega_{i}
xkTxk
=xkT∑i=1d(ωiTxk)ωi
=
∑
i
=
1
d
(
ω
i
T
x
k
)
x
k
T
ω
i
=\sum_{i=1}^{d}\left(\omega_{i}^{\mathrm{T}} x_{k}\right) x_{k}^{\mathrm{T}} \omega_{i}
=∑i=1d(ωiTxk)xkTωi
=
∑
i
=
1
d
ω
i
T
x
k
x
k
T
ω
i
=\sum_{i=1}^{d} \omega_{i}^{\mathrm{T}} x_{k} x_{k}^{\mathrm{T}} \omega_{i}
=∑i=1dωiTxkxkTωi
x
k
~
T
x
k
~
=
(
∑
i
=
1
d
(
ω
i
T
x
k
)
ω
i
)
T
(
∑
j
=
1
d
(
ω
j
T
x
k
)
ω
j
)
\tilde{x_{k}}^{\mathrm{T}} \tilde{x_{k}}=\left(\sum_{i=1}^{d}\left(\omega_{i}^{\mathrm{T}} x_{k}\right) \omega_{i}\right)^{\mathrm{T}}\left(\sum_{j=1}^{d}\left(\omega_{j}^{\mathrm{T}} x_{k}\right) \omega_{j}\right)
xk~Txk~=(∑i=1d(ωiTxk)ωi)T(∑j=1d(ωjTxk)ωj)
=
∑
i
=
1
d
∑
j
=
1
d
(
(
ω
i
T
x
k
)
ω
i
)
T
(
(
ω
j
T
x
k
)
ω
j
)
=\sum_{i=1}^{d} \sum_{j=1}^{d}\left(\left(\omega_{i}^{\mathrm{T}} x_{k}\right) \omega_{i}\right)^{\mathrm{T}}\left(\left(\omega_{j}^{\mathrm{T}} x_{k}\right) \omega_{j}\right)
=∑i=1d∑j=1d((ωiTxk)ωi)T((ωjTxk)ωj)
注意到,其中
ω
i
T
x
k
\omega_{i}^{\mathrm{T}} x_{k}
ωiTxk 和
ω
j
T
x
k
\omega_{j}^{\mathrm{T}} x_{k}
ωjTxk 表示投影长度, 都是数字。且当
i
≠
j
i \neq j
i=j 时,
ω
i
T
ω
j
=
0
,
\omega_{i}^{\mathrm{T}} \omega_{j}=0,
ωiTωj=0, 因此式 (4.11) 的交叉项中只剩下
d
d
d 项
x
~
k
T
x
k
~
=
∑
i
=
1
d
(
(
ω
i
T
x
k
)
ω
i
)
T
(
(
ω
i
T
x
k
)
ω
i
)
=
∑
i
=
1
d
(
ω
i
T
x
k
)
(
ω
i
T
x
k
)
=
∑
i
=
1
d
(
ω
i
T
x
k
)
(
x
k
T
ω
i
)
=
∑
i
=
1
d
ω
i
T
x
k
x
k
T
ω
i
\begin{aligned} \widetilde{x}_{k}^{\mathrm{T}} \tilde{x_{k}} &=\sum_{i=1}^{d}\left(\left(\omega_{i}^{\mathrm{T}} \boldsymbol{x}_{k}\right) \omega_{i}\right)^{\mathrm{T}}\left(\left(\omega_{i}^{\mathrm{T}} \boldsymbol{x}_{k}\right) \omega_{i}\right)=\sum_{i=1}^{d}\left(\omega_{i}^{\mathrm{T}} \boldsymbol{x}_{k}\right)\left(\omega_{i}^{\mathrm{T}} \boldsymbol{x}_{k}\right) \\ &=\sum_{i=1}^{d}\left(\omega_{i}^{\mathrm{T}} \boldsymbol{x}_{k}\right)\left(\boldsymbol{x}_{k}^{\mathrm{T}} \boldsymbol{\omega}_{i}\right)=\sum_{i=1}^{d} \boldsymbol{\omega}_{i}^{\mathrm{T}} \boldsymbol{x}_{k} \boldsymbol{x}_{k}^{\mathrm{T}} \boldsymbol{\omega}_{i} \end{aligned}
x
kTxk~=i=1∑d((ωiTxk)ωi)T((ωiTxk)ωi)=i=1∑d(ωiTxk)(ωiTxk)=i=1∑d(ωiTxk)(xkTωi)=i=1∑dωiTxkxkTωi
注意到
,
∑
i
=
1
d
ω
i
T
x
k
x
k
T
ω
i
, \sum_{i=1}^{d} \omega_{i}^{\mathrm{T}} \boldsymbol{x}_{k} \boldsymbol{x}_{k}^{\mathrm{T}} \omega_{i}
,∑i=1dωiTxkxkTωi 实际上就是矩阵
W
T
x
k
x
k
T
W
\boldsymbol{W}^{\mathrm{T}} \boldsymbol{x}_{k} \boldsymbol{x}_{k}^{\mathrm{T}} \boldsymbol{W}
WTxkxkTW 的迹(对角线
元素之和 ),于是可以将式(4.9)继续化简
∥
x
k
−
x
k
~
∥
2
2
=
−
∑
i
=
1
d
ω
i
T
x
k
x
k
T
ω
i
+
x
k
T
x
k
=
−
t
r
(
W
T
x
k
x
k
T
W
)
+
x
k
T
x
k
\begin{aligned} \| x_{k}-\widetilde{x_{k}} & \|_{2}^{2}=-\sum_{i=1}^{d} \omega_{i}^{\mathrm{T}} x_{k} x_{k}^{\mathrm{T}} \omega_{i}+x_{k}^{\mathrm{T}} x_{k} \\ &=-t r\left(W^{\mathrm{T}} x_{k} x_{k}^{\mathrm{T}} W\right)+x_{k}^{\mathrm{T}} x_{k} \end{aligned}
∥xk−xk
∥22=−i=1∑dωiTxkxkTωi+xkTxk=−tr(WTxkxkTW)+xkTxk
因此式 (4.8) 可以写成
arg
min
W
∑
k
=
1
n
∥
x
k
−
x
k
~
∥
2
2
=
∑
k
=
1
n
(
−
tr
(
W
T
x
k
x
k
T
W
)
+
x
k
T
x
k
)
=
−
∑
k
=
1
n
t
r
(
W
T
x
k
x
k
T
W
)
+
C
\begin{aligned} \arg \min _{W} \sum_{k=1}^{n}\left\|\boldsymbol{x}_{k}-\widetilde{\boldsymbol{x}_{k}}\right\|_{2}^{2} &=\sum_{k=1}^{n}\left(-\operatorname{tr}\left(\boldsymbol{W}^{\mathrm{T}} \boldsymbol{x}_{k} \boldsymbol{x}_{k}^{\mathrm{T}} \boldsymbol{W}\right)+\boldsymbol{x}_{k}^{\mathrm{T}} \boldsymbol{x}_{k}\right) \\ &=-\sum_{k=1}^{n} t r\left(\boldsymbol{W}^{\mathrm{T}} \boldsymbol{x}_{k} \boldsymbol{x}_{k}^{\mathrm{T}} \boldsymbol{W}\right)+C \end{aligned}
argWmink=1∑n∥xk−xk
∥22=k=1∑n(−tr(WTxkxkTW)+xkTxk)=−k=1∑ntr(WTxkxkTW)+C
根据矩阵乘法的性质
∑
k
x
k
x
k
T
=
X
X
T
,
\sum_{k} x_{k} x_{k}^{\mathrm{T}}=X X^{\mathrm{T}},
∑kxkxkT=XXT, 因此优化问题可以转化为
arg
max
W
∑
k
=
1
n
tr
(
W
T
x
k
x
k
T
W
)
,
\arg \max _{W} \sum_{k=1}^{n} \operatorname{tr}\left(W^{\mathrm{T}} x_{k} x_{k}^{\mathrm{T}} W\right),
argmaxW∑k=1ntr(WTxkxkTW), 这等价于求解带约束的优化问题
{
arg
max
W
tr
(
W
T
X
X
T
W
)
s.t.
W
T
W
=
I
\left\{\begin{array}{l} \arg \max _{W} \operatorname{tr}\left(W^{\mathrm{T}} X X^{\mathrm{T}} W\right) \\ \text { s.t. } \quad W^{\mathrm{T}} \boldsymbol{W}=I \end{array}\right.
{argmaxWtr(WTXXTW) s.t. WTW=I
如果我们对
W
W
W 中的
d
d
d 个基
ω
1
,
ω
2
,
…
,
ω
d
\omega_{1}, \omega_{2}, \ldots, \omega_{d}
ω1,ω2,…,ωd 依次求解,就会发现和最大方差理论的方法完全等价。比如当
d
=
1
d=1
d=1 时,我们实际求解的问题是
{
arg
max
U
ω
T
X
X
T
ω
s.t.
ω
T
ω
=
1
\left\{\begin{array}{l} \arg \max _{\mathscr{U}} \omega^{\mathrm{T}} X X^{\mathrm{T}} \omega \\ \text {s.t.} \quad \omega^{\mathrm{T}} \omega=1 \end{array}\right.
{argmaxUωTXXTωs.t.ωTω=1
最佳直线
ω
\omega
ω 与最大方差法求解的最佳投影方向一致,即协方差矩阵的最大特征值所对应的特征向量,差别仅是协方差矩阵
Σ
\Sigma
Σ 的一个倍数,以及常数
∑
k
=
1
n
x
k
T
x
k
\sum_{k=1}^{n} x_{k}^{\mathrm{T}} x_{k}
∑k=1nxkTxk 偏差,但这并不影响我们对最大值的优化。
总结与拓展
至此, 我们从最小平万误差的角度解释了 PCA 的原理、目标函数和求解方法。不难发这与最大方差角度殊途同归,从不同的目标函数出发,得到了相同的求解方法。
03:线性判别分析
线性判别分析( Linear Discriminant Analysis, LDA )是一种有监督学习算法,同时经常被用来对数据进行降维 。 它是Ronald Fisher 在 1936 年发明的,有些资料上也称之为Fisher LDA (Fisher ’ s Linear Discriminant Ana l ysis ) 。 LDA 是目前机器学习、数据拮据领域中经典旦热门的一种算法 。
相比于 PCA , LDA 可以作为一种有监督的降维算法 。 在 PCA 中, 算法没有考虑数据的标签(类别),只是把原数据映射到一些方差比较大的方向上而已。
假设用不同的颜色标注
C
1
,
C
2
C_{1}, C_{2}
C1,C2 两个不同类别的数据,如图 4.4 所示。根据
P
C
A
\mathrm{PCA}
PCA 算法,数据应该映射到方差最大的那个方向,亦即
y
y
y 轴方向。但是,
C
1
,
C
2
C_{1}, C_{2}
C1,C2 两个不同类别的数据就会完全混合在一起,很难区分开。所以,使用 PCA 算法进行降维后再进行分类的效果会非常差。但是,如果使用 LDA 算法,数据会映射到
x
x
x 轴方向。那么,LDA 算法究竟是如何做到这一点的呢?
知识点:线性代数, LDA
问题:对于具有类别标签的数据,应当如何设计目标函数使得降维的过程中不损失类别信息?在这种目标下,应当如何避行求解?
LDA 首先是为了分类服务的,因此只要找到一个投影方向
ω
\omega
ω,使得投影后的样本尽可能按照原始类别分开。我们不妨从一个简单的二分类问题出发,有
C
1
、
C
2
C_{1} 、 C_{2}
C1、C2 两个类别的样本,两类的均值分别为
μ
1
=
1
N
1
∑
x
∈
C
1
x
,
μ
2
=
1
N
2
∑
x
∈
C
2
x
∘
\mu_{1}=\frac{1}{N_{1}} \sum_{x \in C_{1}} x, \mu_{2}=\frac{1}{N_{2}} \sum_{x \in C_{2}} x_{\circ}
μ1=N11∑x∈C1x,μ2=N21∑x∈C2x∘ 我们希望投影之后两类之间的距离尽可能大,距离表示为
D
(
C
1
,
C
2
)
=
∥
μ
~
1
−
μ
~
2
∥
2
2
D\left(C_{1}, C_{2}\right)=\left\|\widetilde{\mu}_{1}-\widetilde{\mu}_{2}\right\|_{2}^{2}
D(C1,C2)=∥μ
1−μ
2∥22
其中
μ
~
1
,
μ
~
2
\tilde{\mu}_{1}, \widetilde{\mu}_{2}
μ~1,μ
2 表示两类的中心在
ω
\omega
ω 方向上的投影向量,
μ
~
1
=
ω
T
μ
1
\tilde{\mu}_{1}=\omega^{\mathrm{T}} \mu_{1}
μ~1=ωTμ1,
μ
2
~
=
ω
T
μ
2
,
\widetilde{\mu_{2}}=\omega^{\mathrm{T}} \mu_{2},
μ2
=ωTμ2, 因此需要优化的问题为
{
max
ω
∥
ω
T
(
μ
1
−
μ
2
)
∥
2
2
s.t.
ω
T
ω
=
1
\left\{\begin{array}{l} \max _{\omega}\left\|\omega^{\mathrm{T}}\left(\mu_{1}-\mu_{2}\right)\right\|_{2}^{2} \\ \text { s.t. } \quad \omega^{\mathrm{T}} \omega=1 \end{array}\right.
{maxω∥∥ωT(μ1−μ2)∥∥22 s.t. ωTω=1
容易发现,当 ω \omega ω 方向与 ( μ 1 − μ 2 ) \left(\mu_{1}-\mu_{2}\right) (μ1−μ2) 一致的时候,该距离达到最大值,例如对图 4.5 ( a ) 的黄棕两种类别的样本点进行降维时,若按照最大化两类投影中心距离的准则,会将样本点投影到下方的黑线上。但是原本可以被线性划分的两类样本,经过投影后有了一定程度的重叠,这显然不能使我们满意。
我们希望得到的投影结果如图 4.5 ( b ) 4.5(\mathrm{b}) 4.5(b) 所示,虽然两类的中心在投影之后的距离有所减小,但确使投影之后样本的可区分性提高了。仔细观察两种投影方式的区别,可以发现,在图 4.5 ( b ) 中,投影后的样本点似乎在每一类中分布得更为集中了,用数学化的语言描述就是每类内部的方差比左图中更小。这就引出了 LDA 的中心思想最大化类间距离和最小化类内距离。
在前文中我们已经找到了使得类间距离尽可能大的投影方式,现在只需要同时优化类内方差,使其尽可能小。我们将整个数据集的类内方差定义为各个类分别的方差之和,将目标函数定义为类间距离和类内距离的比值,于是引出我们需要最大化的目标
max
ω
J
(
ω
)
=
∥
ω
T
(
μ
1
−
μ
2
)
∥
2
2
D
1
+
D
2
\max _{\omega} J(\omega)=\frac{\left\|\omega^{\mathrm{T}}\left(\mu_{1}-\mu_{2}\right)\right\|_{2}^{2}}{D_{1}+D_{2}}
ωmaxJ(ω)=D1+D2∥∥ωT(μ1−μ2)∥∥22
其中
ω
\omega
ω 为单位向量,
D
1
,
D
2
D_{1}, D_{2}
D1,D2 分别表示两类投影后的方差
D
1
=
∑
x
∈
C
1
(
ω
T
x
−
ω
T
μ
1
)
2
=
∑
x
∈
C
1
ω
T
(
x
−
μ
1
)
(
x
−
μ
1
)
T
ω
D
2
=
∑
x
∈
C
2
ω
T
(
x
−
μ
2
)
(
x
−
μ
2
)
T
ω
\begin{array}{c} D_{1}=\sum_{x \in C_{1}}\left(\omega^{\mathrm{T}} x-\omega^{\mathrm{T}} \mu_{1}\right)^{2}=\sum_{x \in C_{1}} \omega^{\mathrm{T}}\left(x-\mu_{1}\right)\left(x-\mu_{1}\right)^{\mathrm{T}} \omega \\ D_{2}=\sum_{x \in C_{2}} \omega^{\mathrm{T}}\left(x-\mu_{2}\right)\left(x-\mu_{2}\right)^{\mathrm{T}} \omega \end{array}
D1=∑x∈C1(ωTx−ωTμ1)2=∑x∈C1ωT(x−μ1)(x−μ1)TωD2=∑x∈C2ωT(x−μ2)(x−μ2)Tω
因此
J
(
ω
)
J(\omega)
J(ω) 可以写成
J
(
ω
)
=
ω
T
(
μ
1
−
μ
2
)
(
μ
1
−
μ
2
)
T
ω
∑
x
∈
C
i
ω
T
(
x
−
μ
i
)
(
x
−
μ
i
)
T
ω
J(\omega)=\frac{\omega^{\mathrm{T}}\left(\mu_{1}-\mu_{2}\right)\left(\mu_{1}-\mu_{2}\right)^{\mathrm{T}} \omega}{\sum_{x \in C_{i}} \omega^{\mathrm{T}}\left(x-\mu_{i}\right)\left(x-\mu_{i}\right)^{\mathrm{T}} \omega}
J(ω)=∑x∈CiωT(x−μi)(x−μi)TωωT(μ1−μ2)(μ1−μ2)Tω
定义类 间 散 度 矩 阵
S
B
=
(
μ
1
−
μ
2
)
(
μ
1
−
μ
2
)
T
,
S_{B}=\left(\mu_{1}-\mu_{2}\right)\left(\mu_{1}-\mu_{2}\right)^{\mathrm{T}},
SB=(μ1−μ2)(μ1−μ2)T, 类 内 散 度 矩 阵
S
w
=
∑
x
∈
C
i
(
x
−
μ
i
)
(
x
−
μ
i
)
T
∘
S_{w}=\sum_{x \in C_{i}}\left(x-\mu_{i}\right)\left(x-\mu_{i}\right)^{\mathrm{T}} \circ
Sw=∑x∈Ci(x−μi)(x−μi)T∘ 则式 (4.22) 可以写为
J
(
ω
)
=
ω
T
S
B
ω
ω
T
S
w
ω
J(\omega)=\frac{\omega^{\mathrm{T}} \boldsymbol{S}_{B} \omega}{\omega^{\mathrm{T}} \boldsymbol{S}_{w} \omega}
J(ω)=ωTSwωωTSBω
我们要最大化
J
(
ω
)
J(\omega)
J(ω) ,只需对
ω
\omega
ω 求偏导,并令导数等于零
∂
J
(
ω
)
∂
ω
=
(
∂
ω
T
S
B
ω
∂
ω
ω
T
S
w
ω
−
∂
ω
T
S
w
ω
∂
ω
ω
T
S
B
ω
)
(
ω
T
S
w
ω
)
2
⋅
=
0
(
4.24
)
\frac{\partial J(\omega)}{\partial \omega}=\frac{\left(\frac{\partial \omega^{\mathrm{T}} S_{B} \omega}{\partial \omega} \omega^{\mathrm{T}} S_{w} \omega-\frac{\partial \omega^{\mathrm{T}} S_{w} \omega}{\partial \omega} \omega^{\mathrm{T}} S_{B} \omega\right)}{\left(\omega^{\mathrm{T}} S_{w} \omega\right)^{2} \cdot}=0 \quad(4.24)
∂ω∂J(ω)=(ωTSwω)2⋅(∂ω∂ωTSBωωTSwω−∂ω∂ωTSwωωTSBω)=0(4.24)
于是得出,
(
ω
T
S
w
ω
)
S
B
ω
=
(
ω
T
S
B
ω
)
S
w
ω
\left(\omega^{\mathrm{T}} S_{w} \omega\right) S_{B} \omega=\left(\omega^{\mathrm{T}} S_{B} \omega\right) S_{w} \omega
(ωTSwω)SBω=(ωTSBω)Swω
由于在简化的二分类问题中
ω
T
S
w
ω
\omega^{\mathrm{T}} S_{w} \omega
ωTSwω 和
ω
T
S
B
ω
\omega^{\mathrm{T}} S_{B} \omega
ωTSBω 是两个数,我们令
λ
=
J
(
ω
)
=
ω
T
S
B
ω
ω
T
S
w
ω
,
\lambda=J(\omega)=\frac{\omega^{\mathrm{T}} S_{B} \omega}{\omega^{\mathrm{T}} S_{w} \omega},
λ=J(ω)=ωTSwωωTSBω, 于是可以把式 (4.25) 写成如下形式
S
B
ω
=
λ
S
w
ω
S_{B} \omega=\lambda S_{w} \omega
SBω=λSwω
整理得,
S
w
−
1
S
B
ω
=
λ
ω
S_{w}^{-1} S_{B} \omega=\lambda \omega
Sw−1SBω=λω
从这里我们可以看出, 我们最大化的目标对应了一个矩阵的特征值,于是 LDA 降维变成了一个求矩阵特征向量的问题。
J
(
ω
)
J(\omega)
J(ω) 就对应了矩阵
S
w
−
1
S
B
S_{w}^{-1} S_{B}
Sw−1SB 最大的特征值,而投影方向就是这个特征值对应的特征向量。对于二分类这一问题,由于
S
B
=
(
μ
1
−
μ
2
)
(
μ
1
−
μ
2
)
T
,
S_{B}=\left(\mu_{1}-\mu_{2}\right)\left(\mu_{1}-\mu_{2}\right)^{\mathrm{T}},
SB=(μ1−μ2)(μ1−μ2)T, 因此
S
B
ω
S_{B} \omega
SBω 的方向始终与
(
μ
1
−
μ
2
)
\left(\mu_{1}-\mu_{2}\right)
(μ1−μ2) 一致,如果只考虑
ω
\omega
ω 的方向,不考虑其长度,可以得到
ω
=
S
w
−
1
(
μ
1
−
μ
2
)
∘
\omega=S_{w}^{-1}\left(\mu_{1}-\mu_{2}\right)_{\circ}
ω=Sw−1(μ1−μ2)∘ 换句话说,我们只需要求样本的均值和类内方差,就可以马上得出最佳的投影方向
ω
∘
\omega_{\circ}
ω∘ 这便是 Fisher 在 1936 年提出的线性判别分析。
总结与拓展
至此,我们从最大化类间距离、最小化类内距离的思想出发,推导出了 LDA 的优化目标以及求解方法。Fisher LDA 相比 PCA 更善于对有类别信息的数据进行降维处理,但它对数据的分布做了一些很强的假设,例如,每个类数据都是高斯分布、各个类的协方差相等。尽管这些假设在实际中并不一定完全满足,但 LDA 已被证明是非常有效的一种降维方法。主要是因为线性模型对于噪声的鲁棒性比较好,但由于模型简单,表达能力有一定局限性,我们可以通过引入核函数扩展 LDA 方法以处理分布较为复杂的数据。
04:线性判别分析与主成分分析
场景描述:同样作为线性降维方法, PCA 是有监督的降维维算法,而 LDA 是无监督的降维算法 。虽然在原理或应用方面二者有一定的区别,但是从这两种方法的数学本质出发,我们不难发现二者有很多共通的特性 。
知识点:线性代数、PCA、LDA
问题:LDA 和 PCA 作为经典的降维算法,如何从应用的角度分析其原理的异同?从数学推导的角度,两种降维算法在目标函数上有何区别与联系?
首先将 LDA 扩展到多类高维的情况,以和问题 1 中
P
C
A
\mathrm{PCA}
PCA 的求解对应。假设有
N
N
N 个类别,并需要最终将特征降维至
d
d
d 维。因此,我们要找到一个
d
d
d 维投影超平面
W
=
{
ω
1
,
ω
2
,
…
,
ω
d
}
,
W=\left\{\omega_{1}, \omega_{2}, \ldots, \omega_{d}\right\},
W={ω1,ω2,…,ωd}, 使得投影后的样本点满足 LDA 的目标一一最大化类间距离和最小化类内距离。
回顾两个散度矩阵,类内散 度 矩 阵
S
w
=
∑
x
∈
C
i
(
x
−
μ
i
)
(
x
−
μ
i
)
T
S_{w}=\sum_{x \in C_{i}}\left(x-\mu_{i}\right)\left(x-\mu_{i}\right)^{\mathrm{T}}
Sw=∑x∈Ci(x−μi)(x−μi)T
在类别增加至
N
N
N 时仍满足定义,而之前两类问题的类间散度矩阵
S
b
=
(
μ
1
−
μ
2
)
(
μ
1
−
μ
2
)
T
S_{b}=\left(\mu_{1}-\mu_{2}\right)\left(\mu_{1}-\mu_{2}\right)^{\mathrm{T}}
Sb=(μ1−μ2)(μ1−μ2)T 在类别增加后就无法按照原始定义。图 4.6 是三类样本的分布情况, 其中
μ
1
,
μ
2
,
μ
3
\mu_{1}, \mu_{2}, \mu_{3}
μ1,μ2,μ3 分别表示棕绿黄三类样本的中心,
μ
\mu
μ 表示这三个中心的均值 ( 也即全部样本的中心
)
,
S
w
i
), S_{w i}
),Swi 表示第
i
i
i 类的类内散度。我们可以定义一个新的矩阵
S
t
,
S_{t},
St, 来表示全局整体的散度,称为全局散度矩阵
S
t
=
∑
i
=
1
n
(
x
i
−
μ
)
(
x
i
−
μ
)
T
S_{t}=\sum_{i=1}^{n}\left(x_{i}-\mu\right)\left(x_{i}-\mu\right)^{\mathrm{T}}
St=i=1∑n(xi−μ)(xi−μ)T
如果把全局散度定义为类内散度与类间散度之和,即
S
i
=
S
b
+
S
w
,
S_{i}=S_{b}+S_{w},
Si=Sb+Sw, 那
么类间散度矩阵可表示为
S
b
=
S
t
−
S
w
=
∑
i
=
1
n
(
x
i
−
μ
)
(
x
i
−
μ
)
T
−
∑
x
∈
C
i
(
x
−
μ
i
)
(
x
−
μ
i
)
T
=
∑
j
=
1
N
(
∑
x
∈
C
j
(
x
−
μ
)
(
x
−
μ
)
T
−
∑
x
∈
C
j
(
x
−
μ
j
)
(
x
−
μ
j
)
T
)
=
∑
j
=
1
N
m
j
(
μ
j
−
μ
)
(
μ
j
−
μ
)
T
\begin{aligned} S_{b} &=S_{t}-S_{w} \\ &=\sum_{i=1}^{n}\left(x_{i}-\mu\right)\left(x_{i}-\mu\right)^{\mathrm{T}}-\sum_{x \in C_{i}}\left(x-\mu_{i}\right)\left(x-\mu_{i}\right)^{\mathrm{T}} \\ &=\sum_{j=1}^{N}\left(\sum_{x \in C_{j}}(x-\mu)(x-\mu)^{\mathrm{T}}-\sum_{x \in C_{j}}\left(x-\mu_{j}\right)\left(x-\mu_{j}\right)^{\mathrm{T}}\right) \\ &=\sum_{j=1}^{N} m_{j}\left(\mu_{j}-\mu\right)\left(\mu_{j}-\mu\right)^{\mathrm{T}} \end{aligned}
Sb=St−Sw=i=1∑n(xi−μ)(xi−μ)T−x∈Ci∑(x−μi)(x−μi)T=j=1∑N⎝⎛x∈Cj∑(x−μ)(x−μ)T−x∈Cj∑(x−μj)(x−μj)T⎠⎞=j=1∑Nmj(μj−μ)(μj−μ)T
其中
m
j
m_{j}
mj 是第
j
j
j 个类别中的样本个数,
N
N
N 是总的类别个数。从式 (4.29)可以看出,类间散度表示的就是每个类别中心到全局中心的一种加权距离。我们最大化类间散度实际上优化的是每个类别的中心经过投影后离全局中心的投影足够远。
根据 LDA 的原理,可以将最大化的目标定义为
J
(
W
)
=
t
r
(
W
T
S
b
W
)
t
r
(
W
T
S
w
W
)
J(W)=\frac{t r\left(W^{\mathrm{T}} S_{b} W\right)}{t r\left(W^{\mathrm{T}} S_{w} W\right)}
J(W)=tr(WTSwW)tr(WTSbW)
其中
W
W
W 是需要求解的投影超平面
,
W
T
W
=
I
,
, W^{\mathrm{T}} W=I,
,WTW=I, 根据问题 2 和问题 3 中的部分结论,我们可以推导出最大化
J
(
W
)
J(W)
J(W) 对应了以下广义特征值求解的问题
S
b
ω
=
λ
S
w
ω
S_{b} \omega=\lambda S_{w} \omega
Sbω=λSwω
求解最佳投影平面
W
=
{
ω
1
,
ω
2
,
…
,
ω
d
}
W=\left\{\omega_{1}, \omega_{2}, \ldots, \omega_{d}\right\}
W={ω1,ω2,…,ωd} 即求解
S
w
−
1
S
b
S_{w}^{-1} S_{b}
Sw−1Sb 矩阵特征值前
d
d
d 大对应的特征向量组成的矩阵,这就将原始的特征空间投影到了新的
d
d
d 维空间中。至此我们得到了与 PCA 步骤类似,但具有多个类别标签
高维数据的 LDA 求解方法。
( 1 )计算数据集中每个类别样本的均值向量
μ
j
,
\mu_{j},
μj, 及总体均值向量
μ
\mu
μ 。
(2 )计算类内散度矩阵
S
w
,
S_{w},
Sw, 全局散度矩阵
S
t
,
S_{t},
St, 并得到类间散度矩
阵
S
b
=
S
t
−
S
w
∘
S_{b}=S_{t}-S_{w \circ}
Sb=St−Sw∘
(3) 对矩阵
S
w
−
1
S
b
S_{w}^{-1} S_{b}
Sw−1Sb 进行特征值分解,将特征值从大到小排列。
(4 ) 取特征值前
d
d
d 大的对应的特征向量
ω
1
,
ω
2
,
…
,
ω
d
,
\omega_{1}, \omega_{2}, \ldots, \omega_{d},
ω1,ω2,…,ωd, 通过以下映
射将
n
n
n 维样本映射到
d
d
d 维
x
i
′
=
[
ω
1
T
x
i
ω
2
T
x
i
⋮
ω
d
T
x
i
]
\boldsymbol{x}_{i}^{\prime}=\left[\begin{array}{c} \omega_{1}^{\mathrm{T}} \boldsymbol{x}_{i} \\ \omega_{2}^{\mathrm{T}} \boldsymbol{x}_{i} \\ \vdots \\ \omega_{d}^{\mathrm{T}} \boldsymbol{x}_{i} \end{array}\right]
xi′=⎣⎢⎢⎢⎡ω1Txiω2Txi⋮ωdTxi⎦⎥⎥⎥⎤
从
P
C
A
\mathrm{PCA}
PCA 和
L
D
A
\mathrm{LDA}
LDA 两种降维方法的求解过程来看,它们确实有着很
大的相似性,但对应的原理却有所区别。首先从目标出发,PCA 选择的是投影后数据方差最大的方向。由于它是无监督的,因此 PCA 假设方差越大,信息量越多,用主成分来表示原始数据可以去除咒余的维度, 达到降维。而 LDA 选择的是投影后类内方差小、类间方差大的方向。其用到了类别标签信息,为了找到数据中具有判别性的维度,使得原始数据在这些方向上投影后,不同类别尽可能区分开。举一个简单的例子,在语音识别中,我们想从一段音频中提取出人的语音信号,这时可以使用 PCA 先进行降维,过滤掉一些固定频率(方差较小
)
)
) 的背景噪声。但如果我们的需求是从这段音频中区分出声音属于哪个人,那么我们应该使用 LDA 对数据进行降维,使每个人的语音信号具有区分性。另外,在人脸识别领域中,PCA 和 LDA 都会被频繁使用。基于PCA 的人脸识别方法也称为特征脸(Eigenface ) 方法,该方法将人脸图像按行展开形成一个高维向量,对多个人脸特征的协方差矩阵做特征值分解,其中较大特征值对应的特征向量具有与人脸相似的形状,故称为特征脸。Eigenface for Recognition 一文中将人脸用 7 个特征脸表示(见图 4.7 ),于是可以把原始 65536 维的图像特征瞬间降到 7 维,人脸识别在降维后的空间上进行。然而由于其利用 PCA 进行降维,一般情况下保留的是最佳描述特征(主成分 ),而非分类特征。如果我们想要达到更好的人脸识别效果,应该用 LDA 方法对数据集进行降维,使得不同人脸在投影后的特征具有一定区分性。
从应用的角度,我们可以掌握一个基本的原则——对无监督的任务使用 PCA 进行降维 ,对有监督的则应用 LDA。
总结与展望
至此,我们从数学原理、优化目标以及应用场景的角度对比了 PCA 和 LDA 这两种经典的结性降维方法,对于非线性数据,可以通过核映射等方法对二者分别进行扩展以得到更好的降维效果 。 关于特征脸这一降维应用,有兴趣的读者可以拜读最经典的 Eigenface论文[4],更好地理解降维算法的实际应用。