Principal Components Analysis

Principal Components Analysis(主成分分析)

2017年12月09日

PCA出现的原因

我们经常提取到许多特征，希望通过提高维度来增强分类器的分类性能，但是这种“提高”是有一定界限的。超过了这个界限后，分类器的分类性能不升反降，出现了过拟合，且使得计算变得非常复杂。基于此，我们需要一种可以将数据的维度降低(dimensionality reduction)，提高数据的密度，但是其“关键”的维度仍然得以保留的算法——降维算法。PCA就是一种典型的降维算法。理论上说，这样一种降维算法有两种可能：特征选择和特征抽取。特征选择是指选择全部特征的一个子集作为特征向量；特征抽取是通过已有的特征，重新建立一个新的特征子集——PCA技术就是一种特征抽取方法。

如何实现PCA？——PCA的数学推导

如何将这个问题转化为数学语言？

将这个问题转化为数学语言，那就是：如何在低维样本空间中较好地代表高维样本空间中的样本。那么问题来了：如何定量地描述：“较好地代表” 呢？一种思路是考虑线性代数中的特征基变换，使用前 n 个较大的 pairs of eigenvalues and eigenvectors，作为降维基变换矩阵进行 change of basis ，这样将原来高维样本空间中的样本映射到低维的以特征基张成的样本空间(principal subspace)中，达到了降维的效果。但是这样的做法是“最好地代表”吗？考虑这样一个超平面(实际上是一个principal subspace)，使得样本点到这个超平面的距离都足够近，这样在用超平面上的点”代表“样本点时，误差就会足够小。这个性质叫做 Minimum-error formulation，是”较好地代表“的一个要求。这个条件如何用数学语言来描述呢？

我们先来表示出样本点在低维空间中的投影

假设所有的样本数据已经进行了中心化：

∑i=1nxi=0 (样本中心化)

且通过投影变换后，得到的principal subspace 的 basis 是一组标准正交基：

{w1,w2,⋯,wD},符合wiwTj=0,(i,j=1,2,⋯,D),∥wi∥=1将这些基向量写作基变换矩阵W,由于正交矩阵的性质,我们可以等价写为WTW=I.

我们的最终目标就是要找到这个基变换矩阵W:

Wn×D,使得Wn×DWTD×nxn×1=Tn×nxn×1=[Tx]n×1:=zn×1总落在Span{Col(W)}中,其中dim(Span{Col(W)})=D<n=dim(Col(X)),这样我们就可以用Span{Col(W)}的D个基向量即w1,w2,⋯,wD来表示[Tx]n×1，这样我们就只要储存线性组合的D个系数，将原来的n维降到了D维。

即（这一步转化可以参见我的线性代数笔记Chapter6 正交性和最小二乘法 中的 正交投影 一节）我们将数据投影到 D-dimensional subspace 中得到的新样本点为：

$$

\mathbf{z}{n}=\sum{i=1}^{D}\left ( \mathbf{x}{n}^{T}\mathbf{w}{i} \right ) \mathbf{w}{i}=WW^{T}\mathbf{x}{n},\text{注意这里的}\mathbf{z}_{n}\text{仍然是一个}n\text{维向量}
$$

表示 Minimum-error formulation 这个要求——化为最优化问题

考虑到我们需要 ”误差足够小“，即 原来的样本点 与 投影后得到的 新样本点 之间的距离最小。这转化为数学语言就是：

min{∑i=1n∥xn−zn∥2}由WTW=I，XTX是常量得LHS=min{(WWTX−X)T(WWTX−X)}由矩阵乘积的转置得=min{(XWWT−XT)(WWTX−X)}展开得=min{XTWWTWWTX−XTWWTX−XTWWTX+XTX}=min{−XTWWTX}

\其中矩阵X是由高维样本点作为列向量的矩阵

由此，我们就将这个Minimum-error formulation，转化成了一个最优化问题：

min{−tr(WTXXTW)}条件:WTW=I.

解决这个最优化问题

这个最优化问题解决起来是比较容易的：只要对上式使用 Lagrange multiplier method（可以在我高等数学的笔记 Chapter 6 拉格朗日乘数法与条件极值问题 中找到），就可以得到：

我们需要找的基矩阵W应该满足其每个列向量wi符合这样的关系:XXTwi=λiwi

这正好是矩阵XX^{T}的特征值与特征向量的定义式。对于矩阵XX^{T}作 eigenvalues decomposition。由矩阵分析的知识（可以在我的线性代数笔记 Chapter 5 特征值与特征向量 一章中的 charactoristc polynomial 一节中找到）我们知道XX^{T}具有n个复特征值，我们将它们降序排列，取出前D个 pais of eigenvalues and eigenvectors，将 eigenvectors 组成矩阵 W，这就是我们要求的矩阵W。至此，我们已经解决了这个最优化问题，得到的矩阵W——投影基矩阵，将它左乘在高维样本向量x上，就可以得到我们通过 PCA 降维得到的低维样本向量 z。这个向量z满足其在高维空间中与自己的原象x是最近的——满足了 Minimum-error formulation 。

推导中得到的惊喜——等价于 Maximum variance formulation

我们发现矩阵XX^{T}正比于样本的协方差矩阵，我们在取 XX^{T}的前 D 大 eigenvalues 时，实际上是取出了协方差矩阵的D个方向，这些方向上数据的方差要大于其余方向上数据的方差。因此，PCA技术除了满足 Minimum-error formulation 外，自然而然也满足了 Maximum variance formulation，即 最小误差等价于最大方差。

反思推导中的问题

PCA与线性最小二乘法的区别

PCA的优化目标是 压缩后的像与原像 （垂直）距离最近，而线性最小二乘法的优化目标是 预测出来的像与实际值 （铅直）距离最近。如下图所示：

PCA的哲学观点

PCA实际上是将数据投影到了D个互相正交的特征方向上，这种 ”正交“ 导致了得到的各个方向之间是互相无关的。但是，数据真正的特征之间是否是无关的呢？并不尽然，我们这种 ”正交“ 投影，实际上使得数据的各个特征之间互相正交——无关，所以如果数据的特征是相关的，那么我们使用 PCA 就会产生较大的误差。另外，PCA除了参数 D 可以由我们来选择，其余的地方并没有参数，无法进行人为调整，方法比较”僵化“。为此，我提出了一种猜想：通过改变n维样本空间上内积的定义，来给PCA加上权值。（类似于加权 Linear Least Square）这样我们就可以人为指定哪些维度上的信息比较”重要“——在进行Minimum-error formulation时，比较重要的维度上的数据，应该与原来的数据保留地更加紧密。事实上，这正是Weight-PCA(加权PCA)算法。（当我得知真的有加权PCA的时候我觉得自己简直出生晚了）

待补充

度量 PCA 算法所丢失的数据

我们采用:

tr(W)tr(X)=t,其中t是PCA所保留的特征维度占原式数据维度的比例，叫做特征贡献率。

来度量PCA算法所保留的数据，损失的数据用1减去即可。

总结PCA算法的步骤

1. 对所有样本进行中心化
2. 计算样本的协方差矩阵XX^{T}
3. 对协方差矩阵XX^{T}进行 eigenvalues decomposition。
4. 取前D大特征值所对应的特征向量组成投影基矩阵W。

如何做到PCA？——PCA的计算机算法

PCA的 Python 实现

这是一个相对来讲比较复杂的PCA实现：

首先加载相关的库

import numpy as npimport pandas as pdimport matplotlib.pyplot as plt

定义一个均值函数，以便进行中心化：

#计算均值,要求输入数据为numpy的矩阵格式，行表示样本数，列表示特征    def meanX(dataX):    return np.mean(dataX,axis=0)#axis=0表示按照列来求均值，如果输入list,则axis=1

PCA的主要部分：

"""参数：    - XMat：传入的是一个numpy的矩阵格式，行表示样本数，列表示特征        - k：PCA唯一的参数——目标维度返回值：    - finalData：进行中心化后(减去均值)后的低维数据矩阵    - reconData：进行逆中心化(加上均值)后的低维数据矩阵"""def pca(XMat, k):    average = meanX(XMat)               #求样本均值    m, n = np.shape(XMat)               #使用shape方法获得原始数据矩阵的行数m和列数n    data_adjust = []    avgs = np.tile(average, (m, 1))     #将中心化后的矩阵average按照m行1列复制(以便m个样本向量可以进行矩阵减法)    data_adjust = XMat - avgs           #样本减去均值，以达到中心化    covX = np.cov(data_adjust.T)   #计算协方差矩阵    featValue, featVec=  np.linalg.eig(covX)  #求解协方差矩阵的特征值和特征向量    index = np.argsort(-featValue) #按照featValue进行从大到小排序    finalData = []    if k > n:        print "k must lower than feature number"        return    else:        #注意特征向量为列向量，而numpy的二维矩阵(数组)a[m][n]中，a[1]表示第1行值        selectVec = np.matrix(featVec.T[index[:k]]) #所以这里需要进行转置——取前k大eigenvectors按行拼成矩阵然后转置为按列拼成矩阵。        finalData = data_adjust * selectVec.T       #进行投影        reconData = (finalData * selectVec) + average      return finalData, reconData

反思

这个算法的时间复杂度非常高，主要在于协方差矩阵的特征对角化。除此之外，如果只是要取出前k大的eigenvalues对应的eigenvectors，大可不必使用o(M^3)的特征对角化，如果使用 The Power Method，只要使用 o(MD^2)就可以获得前k大的eigenvectors.

计算特征值和特征向量——The Power Method.

这个算法可以给以O(M*D^2)的时间复杂度给出方阵的前 k 大的 sequence of vectors（且以标准正交形式给出）.其算法伪代码如下（见 《Matrix Computation》8.2.1 The Power Method）：for k=1,2,...    z^{k}=Aq^{k-1}    q^{k}=z^{k}/||z^{k}||_{2}    lambda^{k}=[q^{k}]^{T}Aq^{k}end

使用Python可以这样来实现 The Power Method:

import numpy as npdef power_iteration(A, num_simulations):    # Ideally choose a random vector    # To decrease the chance that our vector    # Is orthogonal to the eigenvector    b_k = np.random.rand(A.shape[0])    for _ in range(num_simulations):        # calculate the matrix-by-vector product Ab        b_k1 = np.dot(A, b_k)        # calculate the norm        b_k1_norm = np.linalg.norm(b_k1)        # re normalize the vector        b_k = b_k1 / b_k1_norm    return b_k

使用The Power Method可以加快一些PCA的运行速度。

PCA的拓展与应用

Weight-PCA(加权PCA)

哈尔滨工程大学自动化学院的王科俊等曾将WPCA用于步态识别方面，取得了很好的识别效果。用于图像时，与传统的PCA技术不同。传统的PCA技术要先将2维图像矩阵转化为1维图像向量(PCA中高维样本空间中的样本向量x)，然后再进行PCA分析，但是将2维图像矩阵转化为1维图像向量，这个图像向量的维度非常高，使得计算量过于巨大。因此王俊科等采用了2DPCA算法，即直接利用图像矩阵计算图像协方差矩阵的特征向量。

考虑图像矩阵Am×n,将其投影到Xn×d,即得到图像A的投影特征向量YYm×d=Am×nXn×d

我们的最终目标就是找到这个X。考虑PCA的最大方差准则，即我们的优化目标是：矩阵投影后得到的矩阵Y的协方差矩阵S的trace最大。我们可以用数学语言表达我们的优化目标为：

tr{S}=tr{E[(Y−E(Y))(Y−E(Y))T]}=tr{E[(AX−E(AX))(AX−E(AX)T)]}=tr{XTE[(A−EA)T(A−EA)]X}定义E[(A−EA)T(A−EA)]=Gn×n为图像的协方差矩阵，其中E为取均值。

通常，我们可以根据训练样本直接计算矩阵G：

假设有M幅m×n大小的测试样本图像矩阵Ai,其平均图像矩阵为:A¯=1M∑k=1MAk那么Gn×n可以被估计为:G=1M∑k=1M(Ak−A¯)T(Ak−A¯)

这样我们就可以将我们的优化目标写为:

max{J(X)}=max{tr(XTGX)},其中X为归一正交化列向量(只有一个列向量时的投影矩阵).此式即为广义总体散布准则

X是使J(X)达到最优时的最佳投影轴——图像总体散布矩阵最大eigenvalue对应的单位eigenvector

选取协方差矩阵G的前d大eigenvectors x1,x2,⋯,xd，按列组成矩阵X.

这样，我们就获得了我们要找的矩阵X，从而获得了Y，称为图像A的特征矩阵或特征图。但是，2DPCA的数据仍然较高（可见传统PCA处理的数据量有多么巨大），所以王科俊等想到了通过行列相结合的方式进行(2D)^{2}PCA分析。先来看行方向上的2DPCA

记Ak=⎡⎣⎢⎢⎢⎢⎢a(1)ka(2)k⋮a(m)k⎤⎦⎥⎥⎥⎥⎥,A¯=⎡⎣⎢⎢⎢⎢⎢a¯(1)ka¯(2)k⋮a¯(m)k⎤⎦⎥⎥⎥⎥⎥我们就可以将协方差矩阵G表示为:Gt=1M∑k=1M∑t=1m(a(t)k−a¯(t))T(a(t)k−a¯(t))接下来与前面的2DPCA一样，我们可以通过G的前d大特征向量得到最优投影矩阵X

列方向上的2DPCA与行方向上的同理。现在我们将行和列结合起来——假设行方向上的2DPCA产生的投影矩阵为X，列方向产生的为Z。那么我们将矩阵A同时投影到X和Z上，即：

Cq×d=ZTq×mAm×nXn×d，其中C称为特征矩阵

疑问：为什么可以做到这样的同时投影？这是什么数学原理在支持？再来考虑对于一个方向（行方向）上的特征进行加权——W(2D)^{2}PCA算法：

Ỹ =Yλ=[λω1y1,λω2y2,⋯,λωdyd]=[λω1Ax1,λω2Ax2,⋯,λωdAxd]=[A(λω1x1),A(λω2x2),⋯,A(λωdxd)]=AXλω其中，λ1⩾λ2⩾⋯⩾λd为G的前d大eigenvalues.x1,x2,⋯,xd为其对应的eigenvectors.ω∈(0,1)是权值

可以看到，权值的表现形式就是在AX后面乘以omega作为权的对角矩阵lambda.通过引入权值，PCA变成了一个可以根据经验人工调整的方法了。但是，这种幂指数加权的方法存在着以下的问题：

1. 经加权后，大的特征值所对应特征向量 会 淫灭 小特征值所对应的特征向量 所起的作用。
2. 对于(2D)^{2}PCA这种双向投影的算法，仅仅对一个方向进行加权可能会破坏图像的二维结构信息。
   所以，我们可以对于行方向和列方向两个方向进行双向（归一化）加权，这里归一化指的是：给较大的特征值给出较小的权，给较小的特征值较大的权；双向指的是在行方向和列方向同时进行投影，尽可能不破坏图像的二维结构信息。事实上，之所以提出2DPCA，一个重要原因就是传统PCA技术将二维图像变为一维，破坏了图像的二维结构信息。
   首先来解决归一化：
   
   λ′∗i={λ∗iλ∗ωiλ∗i<1λ∗i⩾1,其中λ∗i=λiλd′,(i=1,2,⋯,d)其中λd′可以通过累计特征贡献率t获取。将这些λ′∗i组成对角矩阵取代原来的λω就可以进行归一化加权了。
   
   再来解决双向：只要在列方向上进行类似的操作即可。
   WPCA与2DPCA的反思
   采用(2D)^{2}PCA，其时间复杂度为o(max{m,n})，远低于普通PCA的o(M^3)。这启示我在考虑将一个抽象的方法（如PCA方法）运用到工程领域的实际问题中时，应该考虑的是算法的思想，这样才能意识到普通PCA的样本点是一维向量，而在CV中的样本点是二维图像矩阵，如果压缩为一维就会带来大量的时间复杂度。这就需要在二维矩阵作为样本点的情形下，重新推导PCA方法。由此可见，掌握算法的本质和推导过程是十分重要的。

Quaternion PCA(四元数主成分分析法)

浙江大学电气工程学院的黎云汉等用四元数表示红、黄、蓝，将一幅m×n图像转化为一幅纯四元数图像。在四元数域中，直接对彩色图像进行分析，而不进行灰度化。疑问：目前关于彩色图像处理，都大致有哪些方向呢？目前用于表示彩色图像的彩色模型主要有：RGB模型，CMY（K）模型，HSI模型，HSV模型，YUV模型，YIQ模型和Lab模型等。对于它们的处理方式大致可以分为两类：对于每个分量单独处理和使用四元数进行处理。

参考资料

专著类：

• Christopher Bishop《Pattern Recognition and Machine Learning》
• Springer 《The Elements of Statistical Learning》
• Ian Goodfellow & Yoshua Bengio & Aaron Courveille 《深度学习》
• Golub & Van Loan 《Matrix Computations》
• 《数字图像处理与机器视觉》

课程类：

• 林轩田 《机器学习基石》
• Andrew Ng 《Machine Learning》
• 周志华 《机器学习》

网站类：

• 机器学习笔记033|主成分分析法(PCA)
• PRML读书会第十二章
• Principal Components Analysis
• PCA主成分分析Python实现
• Power Intration

论文类：

• 王科俊等《基于步态能量图像和2维主成分分析的步态识别方法》——中国图象图形学报2009年第14卷第12期
• 黎云汉等《四元数主成分分析的人脸识别算法》——信号处理2007年第23卷第2期
• 管凤旭等《归一双向加权2D2PCA的手指静脉识别方法》