主成分分析（PCA）是一种基于变量协方差矩阵对数据进行压缩降维、去噪的有效方法，PCA的思想是将n维特征映射到k维上（k&n），这k维特征称为主元，是旧特征的线性组合，这些线性组合最大化样本方差，尽量使新的k个特征互不相关。
介绍一个PCA的教程：
1.协方差&Covariance
　　变量X和变量Y的协方差公式如下，协方差是描述不同变量之间的相关关系，协方差&0时说明 X和 Y是正相关关系，协方差&0时 X和Y是负相关关系，协方差为0时 X和Y相互独立。
　　协方差的计算是针对两维的，对于n维的数据集，可以计算C(n,2)种协方差。&n维数据的协方差矩阵的定义如下：　　& & & Dim(x)表示第x维。
& & & 对于三维(x,y,z)，其协方差矩阵如下，可看出协方差矩阵是一个对称矩阵（symmetrical），其对角线元素为每一维的方差：& 　&
2.特征向量和特征值　
　　若，则称是A的特征值，X是对应的特征向量。可以这样理解：矩阵A作用在它的特征向量X上，仅仅使得X的长度发生了变化，缩放比例就是相应的特征值。特征向量只能在方阵中找到，而且并不是所有的方阵都有特征向量，并且如果一个n*n的方阵有特征向量，那么就有n个特征向量。一个矩阵的所有特征向量是正交的，即特征向量之间的点积为0，一般情况下，会将特征向量归一化，即向量长度为1。
　　第一步，获取数据，下图中的Data为原始数据，一共有两个维度，可看出二维平面上的点。
　　下图是Data在二维坐标平面上的散点图：
　　第二步，减去平均值，对于Data中的每一维数据分别求平均值，并减去平均值，得到DataAdjust数据。
　　第三步，计算DataAdjust的协方差矩阵
　　第四步，计算协方差矩阵的特征向量和特征值，选取特征向量
　　特征值0.对应的特征向量是（-0., 0.），这里的特征向量是正交的、归一化的，即长度为1。
　　下图展示DataAdjust数据和特征向量的关系：
　　正号表示预处理后的样本点，斜着的两条线就分别是正交的特征向量（由于协方差矩阵是对称的，因此其特征向量正交），特征值较大的那个特征向量是这个数据集的主要成分（principle component）。通常来说，当从协方差矩阵计算出特征向量之后，下一步就是通过特征值，对特征向量进行从大到小的排序，这将给出成分意义的顺序。成分的特征值越小，其包含的信息量也就越少，因此可以适当选择。　&
　　如果数据中有n维，计算出n个特征向量和特征值，选择前k个特征向量，然后最终的数据集合只有k维，取的特征向量命名为FeatureVector。&
　　这里特征值只有两个，我们选择其中最大的那个，1.，对应的特征向量是。
　　第五步，将样本点投影到选取的特征向量上，得到新的数据集
　　假设样例数为m，特征数为n，减去均值后的样本矩阵为DataAdjust(m*n)，协方差矩阵是n*n，选取的k个特征向量组成的矩阵为EigenVectors(n*k)。那么投影后的数据FinalData为
　　这里是FinalData(10*1) = DataAdjust(10*2矩阵)&特征向量
　　得到结果为
　　下图是FinalData根据最大特征值对应的特征向量转化回去后的数据集形式，可看出是将DataAdjust样本点分别往特征向量对应的轴上做投影：
　　如果取的k=2，那么结果是
　　可见，若使用了所有特征向量得到的新的数据集，转化回去之后，与原来的数据集完全一样（只是坐标轴旋转）。
Python实现PCA
　　将数据转化成前K个主成分的伪码大致如下：　
减去平均数
计算协方差矩阵
计算协方差矩阵的特征值和特征向量
将特征值从大到小排序
保留最大的K个特征向量
将数据转换到上述K各特征向量构建的新空间中
　　代码实现如下：
from numpy import *
def loadDataSet(fileName, delim='\t'):
fr = open(fileName)
stringArr = [line.strip().split(delim) for line in fr.readlines()]
datArr = [map(float,line) for line in stringArr]
return mat(datArr)
def pca(dataMat, topNfeat=999999):
meanVals = mean(dataMat, axis=0)
DataAdjust = dataMat - meanVals
#减去平均值
covMat = cov(DataAdjust, rowvar=0)
eigVals,eigVects = linalg.eig(mat(covMat)) #计算特征值和特征向量
#print eigVals
eigValInd = argsort(eigVals)
eigValInd = eigValInd[:-(topNfeat+1):-1]
#保留最大的前K个特征值
redEigVects = eigVects[:,eigValInd]
#对应的特征向量
lowDDataMat = DataAdjust * redEigVects
#将数据转换到低维新空间
reconMat = (lowDDataMat * redEigVects.T) + meanVals
#重构数据，用于调试
return lowDDataMat, reconMat
　　测试数据由1000个数据点组成。下面对数据进行降维，并用matplotlib模块将降维后的数据和原始数据一起绘制出来。
import matplotlib
import matplotlib.pyplot as plt
dataMat = loadDataSet('testSet.txt')
lowDMat, reconMat = pca(dataMat,1)
print "shape(lowDMat): ",shape(lowDMat)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(dataMat[:,0].flatten().A[0],dataMat[:,1].flatten().A[0],marker='^',s=90)
ax.scatter(reconMat[:,0].flatten().A[0],reconMat[:,1].flatten().A[0],marker='o',s=50,c='red')
plt.show()
　　结果如下图：
Python环境
1.编译环境：win8.1 + 32位 + python2.7
2.相关模块安装：
　　（1）Numpy和Scipy：numpy用于存储和处理大型矩阵，进行数值计算。scipy是基于numpy的科学和工程计算工具，用于处理多维数组向量、矩阵、图形（图形图像是像素的二维数组）、表格。下载地址：
　　（2）Matplotlib：python的一个图形框架，用于数据绘图。win32的安装文件下载地址：
　　（3）Dateutil 和 Pyparsing模块：安装配置matplotlib包时需要。win32的安装文件下载地址：
3.编译遇到问题：
　　（1）提示&No module name six&，将\Python27\Lib\site-packages\scipy\lib中的six.py six.pyc six.pyo三个文件拷贝到\Python27\Lib\site-packages目录下即可。
　　（2）提示 &ImportError:&six&1.3&or&later&is&&you&have&1.2.0&，说明six.py版本过低，下载最新的版本，将其中的six.py替换\Python27\Lib\site-packages中的six.py即可，下载地址：
　　注：six模块是Python 2和3的兼容工具
PCA可以从数据中识别其主要特征，它是通过沿着数据最大方差方向旋转坐标轴来实现的。其中的矩阵运算，求特征值、特征向量等不是很熟悉，仍需进一步学习。原文链接:主成分分析PCA（Principal Component Analysis）介绍 - Ckary的博客
- 博客频道 - CSDN.NET
3706人阅读
很久之前，有一次做人脸识别的时候用过PCA，大概记得是降维用的，然后前段时间用到LDA的时候顺带看到PCA才发现忘的差不多了，干脆把一些资料整理一下吧。
一.K-L变换
说PCA的话，必须先介绍一下K-L变换了。
K-L变换是Karhunen-Loeve变换的简称，是一种特殊的正交变换。它是建立在统计特性基础上的一种变换，有的文献也称其为霍特林（Hotelling）变换，因为他在1933年最先给出将离散信号变换成一串不相关系数的方法。
K-L变换的突出优点是它能去相关性，而且是均方误差（Mean Square Error，MSE）意义下的最佳变换。
下面就简单的介绍一下K-L变换了。
设，随机向量X ∈Rn(n阶列向量)，它的均值向量为mX，则其协方差矩阵可以表示为
Cx= E{(X-mx)*(X-mx)T}&&&&&&&&&&&&&&&&&&&&&&&&&& （2.1）
Cx是一个n*n阶的实对称阵。
K-L变换定义了一正交变换A ∈Rn*n，将X ∈Rn的向量映射到用Y ∈Rn代表的向量，并且使Y向量中各分量间不相关：
Y = A*(X-mx)&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& （2.2）
因为Y的各分量间不相关，则其协方差矩阵Cy为对角阵，即
Cy = diag(λ1,λ2,...,λn)
（CSDN的编辑器不知道怎么打矩阵，diag这个符号忘了的同学自己百度吧）
而矩阵A总是可以找到的，因为对于实对称阵，总能找到一个正交阵A，使得ACxAT的运算结果为对称阵。K-L变换中，将A的每一行取为Cx的特征向量，并且将这些特征向量按对应的特征值大小进行降序排序，使最大特征值对应的特征向量在A的第一行，而最小特征值对应的特征向量在A的最后一行。而Cy是Cx对角化后的结果，所以两个矩阵的特征值是一致的（λ1,λ2,...,λn）。
这样就可以通过矩阵A实现由随机向量X到随机向量Y的K-L变换了，而由
X = ATY +mx&&&&&&&&&&&& & & & & & & & & & & & & & && （2.3）
就可以实现Y反变换到X。
若选择的最大k个特征值对应的k个特征向量，组成k×n的转换矩阵A，则变换后Y降为k维的，则由Y对X的恢复公式如下：
X‘ = AKY +mx&&&&&&&&&&&&& & & & & & & & & & && & &
&& （2.4）
这时候Cy = diag(λ1,λ2,...,λk)，X与X’之间的均方误差可以由下式表达：
λk+1+.λk+2...+λn&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& （2.5）&&&&&&&&&&&&&&&&&&&&&&&&&&& （没有公式编辑器啊）
上面我们提到了对于特征值λ是从大到小排序的，那么这时候通过式子2.5可以表明通过选择k个具有最大特征值的特征向量来降低误差。因此，从可以将向量X和它的近似X‘之间的均方误差降至最小这方面来说，K-L变换是最佳变换。
二.PCA，主成分分析
在二十世纪九十年代初，Kirby和Sirovich开始讨论利用PCA技术进行人脸图像的最优表示问题。并且由M.Turk和A.Pentland将此技术用于人脸识别中，并称为特征脸方法。M.Turk和A.Pentland将m×n的人脸图像，重新排列为m *n维的列向量。则所有的训练图像经此变换后得到一组列向量：{ xi }，xi∈Rm*n，其中N代表训练样本集中图像的个数。将图像看成一随机列向量，并通过训练样本对其均值向量和协方差矩阵进行估计。
均值向量μ通过下式估计：
μ = (1/N)*((x1+x2+...+xN)&&&&&&&&&&&&&&&&&&&&&&& （3.1）
协方差矩阵
ST& = E{(xi-u)*(xi-u)T}& = X'X'T &&&&&&&&&&&&& （3.2）
其中X’ = [x1-μ, x2-μ,...., xN-μ]
则将投影变换矩阵A取为ST的前k个最大特征值对应的特征向量。利用K-L变换式对原图像进行去相关并降维：
Y = AK*(X-mx)&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& （3.3）
因为ST =X'X'T，而X‘为（m*n）*N矩阵，但是因为X’为N阶矩阵，所以ST的秩最大为N-1，这样只要计算出ST的特征向量就可以计算出K-L变换矩阵了。
但是因为ST是（m*n）*(m*n)阶的矩阵，所以计算它的特征向量比较复杂，这里使用了一个技巧：
XTXvi=δivi&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& （3.4）
(XXT)(Xvi)=δi(Xvi) &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& （3.5）
根据式子3.4与3.5可以看出，只要计算出XTX的特征值和特征向量δi与vi，然后就可以计算出XXT的&特征值和特征向量δi与Xvi，而XTX为N*N阶的矩阵，计算起来比较容易，除此以外，也可以使用SVD，这里就不提了。
三.PCA流程整理
PCA的整个变换过程整理了一下，如下：
1.将mxn的训练图像重新排列为m *n维的列向量。计算均值向量，并利用均值向量将所有样本中心化。
2.利用中心化后的样本向量，根据式（3.2）计算其协方差矩阵；对其特征值分解，并将特征向量按其对应的特征值大小进行降序排列。
3.选取第2步所得的k ≤N-1个最大特征值对应的特征向量组成投影矩阵A，将每幅已中心化的训练图像(x1-μ, x2-μ,...., xN-μ)，向矩阵A投影，得到每幅训练图像的降维表示为(y1-μ, y2-μ,...., yN)
4.对测试图像中心化，并投影到矩阵A，得到测试图像的降维表示。
5.选择合适的分类器，对测试图像进行分类。
* 以上用户言论只代表其个人观点，不代表CSDN网站的观点或立场
访问：115799次
积分：1946
积分：1946
排名：第9824名
原创：78篇
评论：41条
(1)(1)(2)(1)(2)(1)(2)(3)(1)(5)(4)(1)(2)(3)(2)(4)(2)(3)(4)(2)(3)(13)(19)(1)再谈协方差矩阵之主成分分析
再谈协方差矩阵之主成分分析& 1 / 4
再谈协方差矩阵之主成分分析
自从上次谈了协方差矩阵之后，感觉写这种科普性文章还不错，那我就再谈一把协方差
矩阵吧。上次那篇文章在理论层次介绍了下协方差矩阵，没准很多人觉得这东西用处不大，
其实协方差矩阵在好多学科里都有很重要的作用，比如多维的正态分布，再比如今天我们今
天的主角——主成分分析(Principal Component Analysis，简称PCA)。结合PCA相信能对协
方差矩阵有个更深入的认识~
PCA 大概是198x年提出来的吧，简单的说，它是一种通用的降维工具。在我们处理高
维数据的时候，为了能降低后续计算的复杂度，在“预处理”阶段通常要先对原始数据进行降
维，而 PCA 就是干这个事的。本质上讲，PCA 就是将高维的数据通过线性变换投影到低维
空间上去，但这个投影可不是随便投投，要遵循一个指导思想，那就是：找出最能够代表原
始数据的投影方法。这里怎么理解这个思想呢？“最能代表原始数据”希望降维后的数据不能
失真，也就是说，被 PCA 降掉的那些维度只能是那些噪声或是冗余的数据。这里的噪声和
冗余我认为可以这样认识：
噪声：我们常说“噪音污染”，意思就是“噪声”干扰我们想听到的真正声音。同样，
假设样本中某个主要的维度 A，它能代表原始数据，是“我们真正想听到的东西”，
它本身含有的“能量”(即该维度的方差，为啥？别急，后文该解释的时候就有啦
~)本来应该是很大的，但由于它与其他维度有那么一些千丝万缕的相关性，受到这
些个相关维度的干扰，它的能量被削弱了，我们就希望通过 PCA处理后，使维度 A
与其他维度的相关性尽可能减弱，进而恢复维度 A应有的能量，让我们“听的更清
冗余：冗余也就是多余的意思，就是有它没它都一样，放着就是占地方。同样，假如样
本中有些个维度，在所有的样本上变化不明显(极端情况：在所有的样本中该维度
都等于同一个数)，也就是说该维度上的方差接近于零，那么显然它对区分不同的
样本丝毫起不到任何作用，这个维度即是冗余的，有它没它一个样，所以 PCA 应该
去掉这些维度。
这么一分析，那么 PCA 的最终目的就是“降噪”和消灭这些“冗余”的维度，以使降低维
度的同时保存数据原有的特征不失真。后面我们将结合例子继续讨论。
协方差矩阵——PCA实现的关键
前面我们说了，PCA的目的就是“降噪”和“去冗余”。“降噪”的目的就是使保留下来的维
度间的相关性尽可能小，而“去冗余”的目的就是使保留下来的维度含有的“能量”即方差尽可
能大。那首先的首先，我们得需要知道各维度间的相关性以及各个维度上的方差啊！那有什
么数据结构能同时表现不同维度间的相关性以及各个维度上的方差呢？自然是非协方差矩
阵莫属。回忆下《浅谈协方差矩阵》的内容，协方差矩阵度量的是维度与维度之间的关系，
而非样本与样本之间。协方差矩阵的主对角线上的元素是各个维度上的方差(即能量)，其他
元素是两两维度间的协方差(即相关性)。我们要的东西协方差矩阵都有了，先来看“降噪”，再谈协方差矩阵之主成分分析&
让保留下的不同维度间的相关性尽可能小，也就是说让协方差矩阵中非对角线元素都基本为
零。达到这个目的的方式自然不用说，线代中奖的很明确——矩阵对角化。而对角化后得到
的矩阵，其对角线上是协方差矩阵的特征值，它还有两个身份：首先，它还是各个维度上的
新方差；其次，它是各个维度本身应该拥有的能量(能量的概念伴随特征值而来)。这也就是
我们为何在前面称“方差”为“能量”的原因。也许第二点可能存在疑问，但我们应该注意到这
个事实，通过对角化后，剩余维度间的相关性已经减到最弱，已经不会再受“噪声”的影响了，
故此时拥有的能量应该比先前大了。看完了“降噪”，我们的“去冗余”还没完呢。对角化后的
协方差矩阵，对角线上较小的新方差对应的就是那些该去掉的维度。所以我们只取那些含有
较大能量(特征值)的维度，其余的就舍掉即可。PCA的本质其实就是对角化协方差矩阵。下
面就让我们跟着上面的感觉来推推公式吧。假设我们有一个样本集 X，里面有N个样本，每
个样本的维度为 d。即：
X={X1,…,XN}, Xi=(xi1,…,xid)∈Rd
, i=1,…,N
将这些样本组织成样本矩阵的形式，即每行为一个样本，每一列为一个维度，得到样本
矩阵 S：S∈RN&d
。我们先将样本进行中心化，即保证每个维度的均值为零，只需让矩阵的
每一列除以对应的均值即可。很多算法都会先将样本中心化，以保证所有维度上的偏移都是
以零为基点的。然后，对样本矩阵计算其协方差矩阵，按照《浅谈协方差矩阵》里末尾的
update，我们知道，协方差矩阵可以简单的按下式计算得到：
&&& C∈Rd&d
下面，根据我们上文的推理，将协方差矩阵 C 对角化。注意到，这里的矩阵 C 是是对
称矩阵，对称矩阵对角化就是找到一个正交矩阵 P，满足：PT
CP=Λ。具体操作是：先对 C
进行特征值分解，得到特征值矩阵(对角阵)即为 Λ，得到特征向量矩阵并正交化即为 P。显
然，P,Λ∈Rd&d
。假如我们取最大的前 p(p
成了新的对角阵 Λ1∈Rp&p
，对应的p个特征向量组成了新的特征向量矩阵P1∈Rd&p
实际上，这个新的特征向量矩阵 P1就是投影矩阵，为什么这么说呢？假设PCA 降维后
的样本矩阵为$S_1$，显然，根据PCA的目的，S1中的各个维度间的协方差基本为零，也就
是说，S1的协方差矩阵应该为Λ1。即满足：
而我们又有公式：
CP=Λ →& P1
代入可得：
→S1=SP1&& S1∈RN&p
由于样本矩阵 SN&d的每一行是一个样本，特征向量矩阵P1(d&p)的每一列是一个特征向量。
右乘 P1 相当于每个样本以 P1 的特征向量为基进行线性变换，得到的新样本矩阵 S1∈RN&p
中每个样本的维数变为了p，完成了降维操作。实际上，P1中的特征向量就是低维空间新的
坐标系，称之为“主成分”。这就是“主成分分析”的名称由来。同时，S1的协方差矩阵 Λ1为
近对角阵，说明不同维度间已经基本独立，噪声和冗余的数据已经不见了。至此，整个PCA
的过程已经结束。小小总结一下：
1.形成样本矩阵，样本中心化
2.计算样本矩阵的协方差矩阵
3.对协方差矩阵进行特征值分解，选取最大的 p个特征值对应的特征向量组成投影矩阵
再谈协方差矩阵之主成分分析& 3 / 4
4.对原始样本矩阵进行投影，得到降维后的新样本矩阵
Matlab 中 PCA 实战
首先，随机产生一个 10*3 维的整数矩阵作为样本集，10 为样本的个数，3 为样本的维
&& S = fix(rand(10,3)*50)
计算协方差矩阵：
&& S=S-repmat(mean(S),10,1);
&& C=(S'*S)./(size(S,1)-1)
&& C=cov(S)
& 257.7333&
-92.1778&& 57.8889
& -92.1778&
145.7333&& 23.2222
23.2222& 295.1111
对协方差矩阵进行特征值分解：
&& [P,Lambda]=eig(C)
-0.5061&& -0.6947
-0.8282&&&
0.6985&& -0.6778
82.4013&&&&&&&&
269.9907&&&&&&&&
0& 346.1857
这里由于三个方差没有明显特别小的，所以我们都保留下来，虽然维度没有降，但可以
将Lambda(即PCA后的样本协方差矩阵)和C(即原始的协方差矩阵)，可以发现，3 个维度上
的方差都有增大，也就是能量都比原来增大了，这就是“降噪”的功劳。最后我们得到降维后
的样本矩阵：
-0.8991&&&
7.7506 再谈协方差矩阵之主成分分析& 4 / 4
-8.0308&&&
2.7188&& 28.1140
-16.6805&&&
-6.2561&&&
6.5610& -27.8154
&& -4.1605&
-33.1376& -11.9527
21.9852&& -1.9115
-8.2271& -11.1225
19.9819& -18.6612
-3.0967&&&
0.5027&& 28.9765
-9.3995&&&
为了验证，我们调用 matlab 自带的主成分分析函数princomp：
&& [COEFF,SCORE]=princomp(S)%COEFF 表示投影矩阵，SCORE 表示投影后新样
-0.5061&&&
-0.2409&&&
0.6985&& -0.2295
-0.8991& -21.2650
-28.1140&&&
&& -0.9157&
-16.6805&& -1.1844
27.8154&&&
&& 11.9527&
-33.1376&&&
21.9852&& -5.8717
-8.2271&& -5.2616
19.9819&&&
-28.9765&&&
-5.7065&&&
对比，可以发现，SCORE和S1在不考虑维度顺序和正负的情况下是完全吻合的，之所
以我们计算的S1的维度顺序不同，是因为通常都是将投影矩阵P 按能量(特征值)的降序排列
的，而刚才我们用 eig 函数得到的结果是升序。另外，在通常的应用中，我们一般是不使
用matlab的princomp函数的，因为它不能真正的降维(不提供相关参数，还是我没发现？)。
一般情况下，我们都是按照协方差矩阵分解后特征值所包含的能量来算的，比如取 90%的
能量，那就从最大的特征值开始加，一直到部分和占特征值总和的 90%为止，此时部分和
含有的特征值个数即为 p。
经过了一番推公式加敲代码的过程，相信大家对主成分分析应该不陌生了吧，同时对协
方差矩阵也有了更深层次的认识了吧，它可不只是花花枪啊。我个人觉得 PCA 在数学上的
理论还是很完备的，想必这也是它能在多种应用中博得鳌头的原因吧。
已投稿到：
以上网友发言只代表其个人观点，不代表新浪网的观点或立场。基于DWT和PCA及LDA的人耳识别
&&&&2014, Vol. 36 Issue (5): 389-393
吕秀丽, 郑江红, 段吉国, 赵丽华, 于波. 基于DWT和PCA及LDA的人耳识别[J]. 光学仪器 ,)389-393.&&
LÜ Xiuli, ZHENG Jianghong, DUAN Jiguo, ZHAO Lihua, YU Bo.
Ear recognition based on DWT,PCA and LDA[J].Optical Instruments ,)389-393.&&
基于DWT和PCA及LDA的人耳识别
吕秀丽1, 郑江红1, 段吉国2, 赵丽华1, 于波1&&&&
1. 东北石油大学 电子科学学院, 黑龙江 大庆 163318;2. 东北石油大学 化学化工学院, 黑龙江 大庆 163318
基金项目：黑龙江省教育厅科学技术项目()
作者简介：吕秀丽(1971-), 女, 教授, 主要从事数字图像处理、模式识别、生物特征识别技术方面的研究. E-mail:lxl_
摘要：根据主成分分析(PCA)和线性判别分析(LDA)在人耳识别过程中存在识别率不高的问题,提出一种基于离散小波变换(DWT)和PCA及LDA的人耳识别算法.将人耳图像进行二维DWT,选择包含图像大部分信息的低频子带,先利用PCA再利用LDA提取最优样本映射空间,最后利用最近邻法则进行人耳图像的分类.实验结果表明,该方法识别效果优于基于PCA及LDA的方法.
人耳识别&&&&
离散小波变换&&&&
主成分分析&&&&
线性判别分析&&&&
Ear recognition based on DWT,PCA and LDA
LÜ Xiuli1, ZHENG Jianghong1, DUAN Jiguo2, ZHAO Lihua1, YU Bo1&&&&
1. College of Electronic Science, Northeast Petroleum University, Daqing 163318, C2. School of Chemistry and Chemical Engineering, Northeast Petroleum University, Daqing 163318, China
Abstract: According to the principal component analysis (PCA) and linear discriminant analysis (LDA) identification accuracy is not high in the ear recognition process. Ear recognition algorithm based on discrete wavelet transform (DWT), PCA and LDA is proposed. The algorithm is decompose the human ear image with two-dimensional DWT, select low-frequency sub-band that contains the most image information, use PCA and LDA in succession to extract the optimal sample mapping space, and use the nearest neighbor rule to classify the human ear image. The experimental results show that this method is superior to identification method of PCA and LDA.
Key words:
ear recognition&&&&
discrete wavelet transform&&&&
principal component analysis&&&&
linear discriminant analysis&&&&
随着科学技术的日新月异,人类对智能化的要求越来越高,人耳识别作为生物识别技术的一个新方向可以单独应用,也可以将人耳和人脸结合进行组合识别,应用前景广泛[, ]。典型人耳识别算法有基于主成分分析(principal component analysis,PCA)、线性判别分析(linear discriminant analysis,LDA)、独立成分分析(independent component analysis,ICA),二维主成分分析(2DPCA)、二维线性判别分析(2DLDA)等,单独使用这些方法会存在识别率低或适应条件差的问题[, , , , ]。本文提出将DWT和PCA及LDA相结合的人耳识别算法。该方法先对人耳图像进行二维小波分解,选择其中的低频子带,并在低频子带利用PCA及LDA降维提取人耳特征,最后用最近邻法则进行分类。实验结果表明,该方法可提高识别率。
1 主成分分析(PCA)
主成分分析[, ] 是一种基于统计分析的特征提取方法,可以将图像的多数特征映射为少数几个特征,本质是将图像的高维数据投影到低维数据空间得到低维数据,使得到的低维数据仍然能够代表原图像的大部分能量,能够保留所需要的识别信息。主成分又称主分量,主成分分析的数学原理如下。
假设N个大小为m×n人耳图像x组成训练样本集{x1,x2,…,xN},将其每列串接构成大小为D=m×n维的列向量,即全部训练样本组成矩阵XD×N,所有训练样本的均值μ为
则样本的协方差矩阵可按照下式计算
主成分分析的最优投影向量Z1,Z2,…,Zd由协方差矩阵S的前d个特征值λ1≥λ2≥…≥λd所对应的特征向量组成。d的选择使前d个主分量的累计贡献率达到某个阈值Te,即满足
对于给定的任意人耳图像向量x,则其主成分为
实际计算中,由于人耳图像训练样本的维数D远大于训练样本总数N,协方差矩阵S较大(大小为D×D),直接由S得到其特征值和特征向量的计算量大,甚至无法计算出结果。所以,根据奇异值分解(SVD)定理,可以通过求解一个大小为N×N的矩阵间接求特征向量[]。
假设N个人耳差值图像=x-μ每列串接构成矩阵D×N,构造矩阵T(大小为N×N),求其特征值λ1≥λ2≥…≥λN及其对应的标准特征向量v1,v2,…,vN,然后根据下式得到前d个最大特征值对应的标准正交投影向量u1,u2,…,ud
给定任意人耳图像向量,则其主成分为
2 线性判别分析(LDA)
线性判别分析[, ]的基本思想是把高维的模式样本投影到最佳鉴别特征空间,从而达到抽取分类信息和压缩特征空间维数的目的。投影后保证模式样本在新的子空间类间距离最大且类内距离最小,即模式在该空间中的可分离性最佳。线性判别分析的数学原理如下。
假设N个d维数据x组成的训练样本集{x1,x2,…,xN}分别属于C类ω1,ω2,…,ωc,类ωi的样本数目为Ni,则类ωi的样本均值μi为
假设Sb表示类间离散度,Sw表示类内离散度,则
LDA作为一个分类算法,希望它所分的类间耦合度低,类内的聚合度高,这样具有较好的分类效果。因此,需要寻找一个最佳投影方向Wopt,使如下的Fisher准则函数的取值最大
在Sw非奇异时,线性判别分析的最佳投影方向Wopt是由矩阵S－1wSb的最大特征值对应的特征向量组成。通过变换
将样本x投影到最佳投影方向上得到投影向量y。实际应用时通常取矩阵S－1wSb的前K个最大特征值对应的特征向量组成最佳投影向量集。
3 离散小波变换(DWT)
小波变换[, ]是现代谱分析工具,与传统的傅里叶变换相比,小波变换是时-频局部分化分析,它通过伸缩平移变换对信号逐步进行多尺寸分析,从而能聚焦到信号的任意细节,成为继傅里叶之后在科学方法上的最大突破。连续小波的定义为:设ψ(t)是基本小波,ψa,b(t)是连续小波函数,对于f(t)∈L2(R),其连续小波为
式中,a(≠0)、b、t均为连续变量,ψ*(t)表示ψ(t)的复共轭。不难发现,连续小波变换具有线性和平移不变性的性质。连续小波的逆小波变换为:如果ψ是允许小波,则对于任何f(t)∈L2(R)和f连续的点x∈R,其中,a,b∈R,a≠0
离散小波变换是对尺度因子a和位移因子b的离散化,一般令参数a=2－j,b=k2－j,其中j,k∈Z,则离散小波为
对应于离散小波ψj,k(t),f(t)∈L2(R)的离散小波变换为
图像经过二维离散小波分解的示意图如图 1所示,其中图 1(b)为对图像I进行一层小波分解后获得的4个不同频带的子带图。LL1子带图表示原始图像的低频信息,包含原图像的主要特征;HL1表示原始图像垂直方向的高频信息;LH1表示原始图像水平方向的高频信息;HH1表示原始图像对角方向的高频信息。图 1(c)是二层小波分解的示意图,其中LL2、HL2、LH2、HH2是将一层小波低频子图进一步分解得到的。图像经过小波分解后,高频分量取值较小,大部分数值趋近于零,低频分量取值较大,集中了图像的大部分能量。本文采用Haar小波作为基函数对待识别图像进行预处理。
图 1Fig. 1
图 1 二维离散小波分解示意图
Fig. 1 Schematic of two-dimensional discrete wavelet decomposition
4 基于DWT和PCA及LDA人耳识别算法
运用PCA、LDA进行人耳特征提取时,二维人耳图像矩阵需要降维成一维向量,通常一维向量的维数远高于样本个数,直接求解投影特征向量困难。可采用DWT对原始图像进行预处理,再将PCA和LDA结合进行人耳识别。本文算法原理框图如图 2所示,算法步骤如下。
图 2Fig. 2
图 2 基于DWT和PCA及LDA人耳识别算法原理框图
Fig. 2 The schematic based on DWT,PCA and LDA algorithm
(1)先将所有人耳图像进行二维DWT,提取低频分量;选择低频分量图像的差值图像做为训练样本集和测试样本集。
(2)由训练样本集计算T矩阵的特征值和特征向量,根据能量阈值Te确定标准特征向量u1,u2,…,ud,将所有样本根据公式P=(u1,u2,…,ud)T降维处理得到其主成分P。
(3)根据训练样本集的主成分求取特征向量类间离散度Sb和类内离散度Sw,计算S－1wSb矩阵的特征值和特征向量,取前K个最大特征值对应的特征向量组成最优投影向量Wopt,将所有样本根据公式=WToptP降维处理得到最佳人耳特征向量。
(4)最后,计算测试样本和训练样本的最佳人耳特征向量之间的欧式距离,并利用最近邻法则分类人耳图像并统计识别率。
5 实验结果及分析
实验采用MATLAB软件编写,在Intel(R)Core(TM)i3 CPU 2.53 GHz,内存2 GB的计算机上运行。图像库采用北京科技大学人耳图像库,该图像库有60人的180张人耳图像,每一张图像为80×150像素。每个人3张人耳图像,3张图像分别为正面人耳图像及角度微变化人耳图像及光照微变化人耳图像。本实验用每个人的正面人耳图像和角度微变化人耳图像作为训练样本,每个人的光照微变化人耳图像作为测试样本。实验对原始人耳图像进行一层小波分解和二层小波分解,分别提取低频分量,再将其进行PCA降维和LDA映射求取最优样本空间,最后利用最近邻法则分类。表 1给出了采用一层小波和二层小波分解再进行PCA及LDA和直接进行PCA及LDA方法的正确识别率实验数据,能量阈值Te取0.94,选取的最大特征值的个数K取5,6,…,17,提出的方法与PCA及LDA方法的识别率比较如图 3所示。
表 1(<span
class="TabText" id="Table 1">Table 1)
表 1 实验测试数据
Tab. 1 The experimental test data
特征值K一层DWT+PCA+LDA二层DWT+PCA+LDAPCA+LDA
50.800 00.766 70.816 7
60.866 70.850 00.833 3
70.866 70.866 70.800 0
80.866 70.850 00.850 0
90.883 30.900 00.900 0
100.900 00.950 00.883 3
110.983 30.933 30.950 0
120.983 30.950 00.950 0
130.950 00.933 30.966 7
140.950 00.933 30.950 0
150.950 00.950 00.916 7
160.950 00.950 00.916 7
170.966 70.950 00.916 7
表 1 实验测试数据
Tab. 1 The experimental test data
图 3Fig. 3
图 3 人耳图像识别率对照图
Fig. 3 The contrast chart of ear recognition rate
由表 1和图 3可知,在相同人耳图像库测试条件下,PCA及LDA加入DWT的方法总体识别率高于基于PCA及LDA的方法;且加入一层DWT前提下,当K取11、12时达到最高识别率98.33%,高于PCA及LDA方法的最高识别率96.67%。实验结果表明,提出的方法识别效果优于PCA及LDA方法。
在人耳图像识别中将离散小波变换应用于图像的预处理,通过提取包含图像主要信息的低频分量,减少了包含人耳图像细节信息的高频噪声分量,降低了原始图像的冗余信息,也大大降低了人耳图像的维数。实验结果表明,采用DWT和PCA及LDA相结合的方法,识别率最高达到了98.33%,显示了本文算法的优越性。
苑玮琦,柯 丽,白 云.生物特征识别技术[M]..
袁 立,穆志纯,徐正光,等.基于人耳生物特征的身份识别[J]..
VICTOR B,BOWYER K,SARKAR S..
徐正光,武 楠,穆志纯.基于独立分量分析的人耳识别方法[J]..
王 杨,程 科.基于PCA与Fisherface互补双特征提取的人耳图像识别[J]..
赵海龙,穆志纯.基于小波变换和LDA/FKT及SVM的人耳识别[J]..
吕秀丽,申屠红峰,等.基于2DLDA与FSVM的人耳识别[J]..
THEODORIDIS S,KOUTROUMBAS K.模式识别[M].4版..
华顺刚,周 羽,刘 婷.基于PCA+LDA的热红外成像人脸识别[J]..
孙延奎.小波变换与图像、图形处理技术[M]..
邸继征.小波分析原理[M]..