模式识别作业三——kl变换

实验二：KL 变换实验学时：4学时实验目的：1. 掌握特征提取的基本方法。

2. 掌握基于KL 变换的特征提取的方法。

3. 培养学生灵活使用KL 变换进行模式识别的能力。

实验内容：给出ORL 人脸数据库，共有400幅人脸图像（40人，每人10幅，大小为92*112象素）。

其中第一个人的图像如下图：选取数据库中的部分样本（每个人的前5张图片）作为训练样本，其余作为未知的测试样本。

从训练样本中得到KL 变换矩阵，然后对训练样本和测试样本都进行变换，用变换后的数据作最近邻识别，距离可以为对应灰度值之差的平方和，统计识别率。

KL 变换方法简介：设图像数据库图像大小为Width ⨯Height ，令d = Width ⨯Height ，则每幅图像可以用按行或者按列堆成一个d 维向量表示。

令111()()N T T t i i t t i S x m x m N N==--=ΦΦ∑，其中1(,,)t N x m x m Φ=-- 。

特征脸方法(KL 变换方法)是从PCA 方法导出的。

PCA 方法就是要寻找正交归一的变换矩阵12(,,,)d L L W u u u R ⨯=∈ ，1T WW =，使得判别准则()()T t J W tr W S W =达到最大，即arg max ()T t WW tr W S W =。

也就是在T y W x =的正交变换后使得总体散度矩阵y T t t S W S W =的迹最大，即各个样本分离得尽量远，将样本的方差尽量的保留下来，和原样本的均方误差尽量小。

可以证明求得12(,,,)L W u u u = 就是对应于矩阵t S 的前L 个最大的特征值的特征向量。

即12(,,,)L W u u u = 的各列向量是下面特征方程的解：t i i i S u u λ=显然，变换后的y T t t S W S W =是一个对角阵，对角线上的元素为相应的特征值，即1()dt i i tr S λ==∑，也就是说每个特征值都代表了相应的特征向量保留总体散度（方差）的能力。

Karhunen-Loeve 变换的应用摘要：本文对Karhunen-Loeve 变换的原理进行了说明，重点分析了K-L 变换的性质，结合K-L 变换的性质，对K-L 变换的具体应用进行了展示。

利用K-L 变换在人脸识别、遥感图像特征提取、地震波噪声抑制、数字图像压缩、语音信号增强中的具体利用，深入总结了K-L 变换在模式识别、噪声抑制和数据压缩领域的重要性。

关键字： Karhunen-Loeve 变换 K-L 变换 K-L 展开1、 Karhunen-Loeve 变换定义1.1Karhunen-Loeve 变换的提出在模式识别和图像处理等现实问题中，需要解决的一个主要的问题就是降维，通常我们选择的特征彼此相关，而在识别这些特征时，数据量大且效率低下。

如果我们能减少特征的数量，即减少特征空间的维数，那么我们将以更少的存储和计算复杂度获得更好的准确性。

于是我们需要一种合理的综合性方法，使得原本相关的特征转化为彼此不相关，并在特征量的个数减少的同时，尽量不损失或者稍损失原特征中所包含的信息。

Karhunen-Loeve 变换也常称为主成分变换(PCA)或霍特林变换，就可以简化大维数的数据集合，而且它的协方差矩阵除对角线以外的元素都是零，消除了数据之间的相关性。

所以可以用于信息压缩、图像处理、模式识别等应用中。

Karhunen-Loeve 变换，是以矢量信号X 的协方差矩阵Ф的归一化正交特征矢量q 所构成的正交矩阵Q ，来对该矢量信号X 做正交变换Y=QX ，则称此变换为K-L 变换（K-LT 或KLT ），K-LT 是Karhuner-Loeve Transform 的简称，有的文献资料也写作KLT 。

可见，要实现KLT ，首先要从信号求出其协方差矩阵Ф，再由Ф求出正交矩阵Q 。

Ф的求法与自相关矩阵求法类似。

1.2Karhunen-Loeve 展开及其性质设零均值平稳随机过程u(n)构成的M 维随机向量为u(n),相应的相关矩阵为R ，则向量u(n)可以表示为R 的归一化特征向量M 21q ,q ,q 的线性组合，即iMi i q c n u ∑==1)(，此式称为u(n)的Karhunen-Loeve 展开式，展开式的系数i c 是由内积 )(c i n u q Hi =M ,1,2,i =定义的随机变量，且有{}0E =i c ，{}⎩⎨⎧≠==l i li c c i li ,0,E *λ。

K_L变换，SVD变换在图像处理中实际应用讨论一．K_L变换在人脸识别中的实际应用模式识别中一个最基本的问题就是特征抽取, 抽取有效的鉴别特征是解决识别问题的关键。

K-L 变换是特征抽取的最有效的方法之一。

基于K -L 变换的特征抽取有四种最为经典和常用的技术, 其中主分量分析( PCA, 或称以总体散布矩阵为产生矩阵的K-L 变换) 被广泛地应用在人脸等图像识别领域。

主分量分析在处理图像识别问题时, 其弱点也很明显: 首先, 将图像矩阵转化为图像向量后, 造成图像向量的维数一般较高, 比如, 图像的分辨率为100×100, 所得图像向量的维数高达10000, 在如此高维的图像向量上完成PCA 是非常耗时的。

尽管在PCA 中利用奇异值分解定理可加速总体散布矩阵的本征向量的求解速度, 但整个特征抽取过程所耗费的计算量相当可观; 其次, 主分量分析的训练是非监督的，即PCA 无法利用训练样本的类别信息。

另一种经典的K-L 变换方法, 即包含在类平均向量中判别信息的最优压缩技术, 其识别性能优于非监督的主分量分析方法。

受Hong代数特征抽取思想的启发, 本研究提出了一种直接基于图像矩阵的包含在类平均图像中判别信息的最优压缩方法——广义K-L 变换, 它克服了传统主分量分析的弱点。

在ORL 标准人脸库上的试验结果表明, 广义K-L变换方法不仅在识别性能上优于Eigenfaces 方法和Fisherfaces 方法, 且特征抽取的速度提高了近19 倍。

1. 1基本思想定义1: 设X 表示n 维列向量, 将m×n 的图像矩阵A 通过以下线性变换直接投影到X 上。

Y= A X ( 1)得到一个m 维列向量Y, X 为投影轴, Y 为图像A的投影特征向量。

设有C 个已知的模式类别, ni 表示第i 类的训练样本数, N 表示各类的训练样本总数。

第i 类的第j 个训练样本图像为m×n 的矩阵Aj ( i=1, 2, ⋯, C; j = 1, 2, ⋯, ni) 。

“模式识别(三).PDF”课件课后上机选做作业参考解答(武大计算机学院袁志勇, Email: yuanzywhu@) 上机题目：两类问题，已知四个训练样本ω1={(0,0)T,(0,1)T}；ω2={(1,0)T,(1,1)T}使用感知器固定增量法求判别函数。

设w1=(1,1,1)Tρk=1试编写程序上机运行(使用MATLAB、 C/C++、C#、JA V A、DELPHI等语言中任意一种编写均可)，写出判别函数，并给出程序运行的相关运行图表。

这里采用MATLAB编写感知器固定增量算法程序。

一、感知器固定增量法的MATLAB函数编写感知器固定增量法的具体内容请参考“模式识别(三).PDF”课件中的算法描述，可将该算法编写一个可以调用的自定义MATLAB函数：% perceptronclassify.m%% Caculate the optimal W by Perceptron%% W1-3x1 vector, initial weight vector% Pk-scalar, learning rate% W -3x1 vector, optimal weight vector% iters - scalar, the number of iterations%% Created: May 17, 2010function [W iters] = perceptronclassify(W1,Pk)x1 = [0 0 1]';x2 = [0 1 1]';x3 = [1 0 1]';x4 = [1 1 1]';% the training sampleWk = W1;FLAG = 0;% iteration flagesiters = 0;if Wk'*x1 <= 0Wk =Wk + x1;FLAG = 1;endif Wk'*x2 <= 0Wk =Wk + x2;FLAG = 1;endif Wk'*x3 >= 0Wk=Wk-x3;FLAG = 1; endif Wk'*x4 >= 0Wk =Wk -x4; FLAG = 1; enditers = iters + 1; while (FLAG) FLAG = 0; if Wk'*x1 <= 0Wk = Wk + x1; FLAG = 1; endif Wk'*x2 <= 0Wk = Wk + x2; FLAG = 1; endif Wk'*x3 >= 0 Wk = Wk - x3; FLAG = 1; endif Wk'*x4 >= 0 Wk = Wk - x4; FLAG = 1; enditers = iters + 1; endW = Wk;二、程序运行程序输入：初始权向量1W , 固定增量大小k ρ 程序输出：权向量最优解W , 程序迭代次数iters 在MATLAB 7.X 命令行窗口中的运行情况： １、初始化1[111]T W = 初始化W 1窗口界面截图如下：２、初始化1kρ=初始化Pk 窗口界面截图如下：３、在MATLAB 窗口中调用自定义的perceptronclassify 函数由于perceptronclassify.m 下自定义的函数文件，在调用该函数前需要事先[Set path…]设置该函数文件所在的路径，然后才能在命令行窗口中调用。

模式识别作业报告组员：2011302265 孔素瑶2011302268 马征2011302273 周昳慧一、实验要求用FAMALE.TXT 和MALE.TXT 的数据作为本次实验使用的样本集，利用K-L 变换对该样本集进行变换，与过去用Fisher 线性判别方法或其它方法得到的分类面进行比较，从而加深对所学内容的理解和感性认识。

二、具体做法1． 不考虑类别信息对整个样本集进行K-L 变换（即PCA ），并将计算出的新特征方向表示在二维平面上，考察投影到特征值最大的方向后男女样本的分布情况并用该主成分进行分类。

2． 利用类平均向量提取判别信息，选取最好的投影方向，考察投影后样本的分布情况并用该投影方向进行分类。

3． 将上述投影和分类情况与以前做的各种分类情况比较，考察各自的特点和相互关系。

三、实验原理设n 维随机向量T n x x x x ),,(21⋯=，其均值向量][x E u =，相关矩阵][Tx xx E R =，协方差矩阵]))([(Tx u x u x E C --=，x 经正交变换后产生向量T n y y y y ),,(21⋯=。

设有标准正交变换矩阵)),,((21n t t t T T ⋯=，（即I T T T=）Tn T n y y y x t t t x T y ),,(),,(2121⋯=⋯==，x t y Ti i = (1,2,)i n =∑=-===ni i i Tt y y T y T x 11（称为 x 的K-L 展开式）取前m 项为x 的估计值1ˆmi i i xy t ==∑ 1m n ≤<其均方误差为∑∑+=+=∧∧==--=nm i ii nm i iTy y E yE x x x x E m 1'12][][)]()[()(ξ∑∑∑+=+=+====nm i i x inm i inm i ii t R ttx x E t y y E m 1'1''1')(][)(ξ在I T T ='的约束条件下,要使均方误差min )]()[()(1''→=--=∑+=∧∧nm i i x it R tx x x x E m ξ为此设定准则函数∑∑+=+=--=nm i i Ti i nm i i x Ti t t t R t J 11)1(λ由0iJt ∂=∂可得()0x i i R I t λ-= 1,...,i m n =+ 即x i i i R t t λ= 1,...,i m n =+表明: λi 是x R 的特征值，而i t 是相应的特征向量。

模式识别作业
4. 已知一组数据的协方差矩阵为[1 1
2
1
2
1]试问：
(1) 协方差矩阵中各元素的含义是什么？
(2) K-L变换的最佳准则是什么？
(3) 为什么说经K-L变换后消除了各分量之间的相关性？
（1）协方差矩阵的x ij表示i号特征和j号特征之间的协方差的值，其大小表示两者的线性相关性。

（2）最佳准则为对一组数据进行按一组正交基分解，在只取相同数量分量的条件下，以均方误差计算截尾误差最小。

（3）因为经过K-L变换，主分量分解后，协方差矩阵成为对角矩阵，因而各主分量间相关性消除。

5.若有下列两类样本集：
ω1: x1=(0,0,0)T, x2=(1,0,0)T,x3=(1,0,1)T,x4=(1,1,0)T
ω2:x5=(0,0,1)T,x6=(0,1,0)T,x7=(0,1,1)T,x8=(1,1,1)T
要求用K-L变换法，分别把特征空间维数降到d=2和d=1。

试编写满足要求的。

模式识别作业三一、编程要求：编程实现KL变换，并对TM六波段图像进行演算。

KL变换的思想是：从n维特征选取m维特征，删去的n-m维特征不一定就是无用的信息，如何在信息损失最小的情况下选取特征，在理论上就显得更严密些。

通常采用正交变换，得到新的经变换的模式，以保证信息损失最小情况下获得有利于分类的特征。

二、编程思想：将6副图象依次输入获得灰度值存在一个6*size(size为一副图象的像素数)的二维数组中，计算每个波段的灰度均值，然后计算协方差矩阵，得出特征值和特征向量，再重新计算要输出的处理的图象的灰度。

三、核心代码：void CKLchangeView::Oninput() //TM图象数据读入{AfxMessageBox("注意:请一次将6副图象都读入!");//读图象的方法int flag=0;while(flag<6){CFileDialog dlg(TRUE,"bmp",".bmp");if(dlg.DoModal()==IDOK){m_NewFileName[flag]=dlg.GetFileName();for(int j=0;j<5;j++){if(flag==j||m_NewFileName[flag]!=m_NewFileName[j]) continue;else{AfxMessageBox("已打开过的图象!请按方法重新打开!");flag=0;for(int i=0;i<6;i++) m_NewFileName[i]="";return;}}NewFileName=m_NewFileName[flag];flag++;}//判断是否是已经打开的图象，若已经打开，则要重新按方法打开else //不读图象，或是只是要看图象{for(int i=0;i<6;i++) m_NewFileName[i]="";return;}CDC *pDC = GetWindowDC();CFile f(NewFileName,0);m_pDib.Read(&f);//读图象if(m_pDib.m_nColorTableEntries!=256)//判断是否是灰度图象{AfxMessageBox("请打开灰度图象!");return;}CKLchangeDoc* pDoc = GetDocument();pDoc->pDib=&m_pDib;OnDraw(pDC);//调用OnDraw函数绘图sizeImage = m_pDib.GetDimensions();sizeSaveImage = m_pDib.GetDibSaveDim();a=sizeSaveImage.cx * sizeSaveImage.cy;//象素数if(flag==1) tm=CMatrix(6,a);//开第一副图象时为矩阵分配内存//将每幅图象象素灰度值放入矩阵中if(flag==1)for(int i=0;i<a;i++){tm.m_pData[i]=*(m_pDib.m_lpImage+sizeSaveImage.cx*(i/sizeImage.cx)+i%sizeImage.cx );}if(flag==2)for(int i=0;i<a;i++){tm.m_pData[a+i]=*(m_pDib.m_lpImage+sizeSaveImage.cx*(i/sizeImage.cx)+i%sizeImage. cx);}if(flag==3)for(int i=0;i<a;i++){tm.m_pData[2*a+i]=*(m_pDib.m_lpImage+sizeSaveImage.cx*(i/sizeImage.cx)+i%sizeImag e.cx);}if(flag==4)for(int i=0;i<a;i++){tm.m_pData[3*a+i]=*(m_pDib.m_lpImage+sizeSaveImage.cx*(i/sizeImage.cx)+i%sizeImag e.cx);}if(flag==5)for(int i=0;i<a;i++){tm.m_pData[4*a+i]=*(m_pDib.m_lpImage+sizeSaveImage.cx*(i/sizeImage.cx)+i%sizeImag e.cx);}if(flag==6)for(int i=0;i<a;i++){tm.m_pData[5*a+i]=*(m_pDib.m_lpImage+sizeSaveImage.cx*(i/sizeImage.cx)+i%sizeImag e.cx);}}}void CKLchangeView::OnKl() //KL变换计算{int i,j;xfea.Init(6,1);xfeat.Init(1,6);temp=CMatrix(6,6);feature=new double[6];fea.Init(6,6);double tmave[6];//求每幅图象的灰度均值for (i=0;i<6;i++){tmave[i]=0;for(j=0;j<a;j++){tmave[i]+=tm.m_pData[i*a+j];}tmave[i]/=(double)a;}//求方差矩阵Rfor(j=0;j<36;j++){temp.m_pData[j]=0.0;for(i=0;i<a;i++){temp.m_pData[j]+=(tm.m_pData[(j/6)*a+i%6]-tmave[j/6])*(tm.m_pData[(j%6)*a+i%6]-tmave[ j%6]);}temp.m_pData[j]/=(double)a;}temp.JacobiEigenv2(feature, fea, 0.000001);//求特征值和特征向量double t=0.0;double featuretemp[6];//临时变量for(i=0;i<6;i++)featuretemp[i]=feature[i];//给临时数组赋值for(i=0;i<6;i++)//"冒泡法"将特征值排序for(j=0;j<5;j++)if(feature[j]>feature[j+1]){t=feature[j];feature[j]=feature[j+1];feature[j+1]=t;}for(i=0;i<6;i++)//求出特征值排序后对应在特征值数组中的位置for(j=0;j<6;j++)if(feature[i]==featuretemp[j]){index[i]=j;}}void CKLchangeView::OnFirst() //第一主分量{if(NewFileName=="") {AfxMessageBox("请先读入图象数据!");return;} if(fea.m_pData[0]==0) OnKl();result.Init(1,a);for(int i=0;i<a;i++){for(int j=0;j<6;j++)result.m_pData[i]+=fea.m_pData[index[5]+j*6]*tm.m_pData[j*a+i];}Draw();}void CKLchangeView::OnFive() //第五主分量{if(NewFileName=="") {AfxMessageBox("请先读入图象数据!");return;} if(fea.m_pData[0]==0) OnKl();result.Init(1,a);for(int i=0;i<a;i++){for(int j=0;j<6;j++)result.m_pData[i]+=fea.m_pData[index[1]+j*6]*tm.m_pData[j*a+i];}Draw();}void CKLchangeView::OnFour() //第四主分量{if(NewFileName=="") {AfxMessageBox("请先读入图象数据!");return;} if(fea.m_pData[0]==0) OnKl();result.Init(1,a);for(int i=0;i<a;i++){for(int j=0;j<6;j++)result.m_pData[i]+=fea.m_pData[index[2]+j*6]*tm.m_pData[j*a+i];}Draw();}void CKLchangeView::OnSecond() //第二主分量{if(NewFileName=="") {AfxMessageBox("请先读入图象数据!");return;} if(fea.m_pData[0]==0) OnKl();result.Init(1,a);for(int i=0;i<a;i++){for(int j=0;j<6;j++)result.m_pData[i]+=fea.m_pData[index[4]+j*6]*tm.m_pData[j*a+i];}Draw();}void CKLchangeView::OnSix() //第六主分量{if(NewFileName=="") {AfxMessageBox("请先读入图象数据!");return;}if(fea.m_pData[0]==0) OnKl();result.Init(1,a);for(int i=0;i<a;i++){for(int j=0;j<6;j++)result.m_pData[i]+=fea.m_pData[index[0]+j*6]*tm.m_pData[j*a+i];}Draw();}void CKLchangeView::OnThird() //第三主分量{if(NewFileName=="") {AfxMessageBox("请先读入图象数据!");return;}if(fea.m_pData[0]==0) OnKl();result.Init(1,a);for(int i=0;i<a;i++){for(int j=0;j<6;j++)result.m_pData[i]+=fea.m_pData[index[3]+j*6]*tm.m_pData[j*a+i];}Draw();}void CKLchangeView::Draw() //图象输出{int i;for(i=0;i<a;i++)//处理越界的灰度值{if(result.m_pData[i]<0) result.m_pData[i]=0;if(result.m_pData[i]>255) result.m_pData[i]=255;}for (i = 0; i <a; i ++)//重新给图象每个象素赋灰度值{ unsigned char *pData=m_pDib.m_lpImage+sizeSaveImage.cx*(i/sizeImage.cx)+i%sizeImage.cx;*pData=(BYTE)(long)(result.m_pData[i]+0.5);}CDC *pDC = GetWindowDC();CKLchangeDoc* pDoc = GetDocument();pDoc->pDib=&m_pDib;OnDraw(pDC);//将处理后图象绘出}。