四元数保结构算法及其在彩色人脸识别中的应用

1、硕士学位论文(20161)四元数保结构算法及其在彩色人脸识别中的应用Quaternion structure-preserving method and itsapplication in color face recognition作者马茹茹导师贾志刚副教授江苏师范大学数学与统计学院学 号：103201 1321028 姓名：马茹茹论文题目：四元数保结构算法及其在彩色人脸识别中的应用中文关键词：四元数实表示矩阵；四元数特征值问题；保结构Jacobi方法； 彩色人脸识别.江苏师范大学学位论文版权使用授权书本学位论文作者完全了解江苏币范大学有权保留并向国家有关部门或机构送交论文的 复印件和电子版

2、，允许论文被查阅和借阅。本人授权江苏4币范大学可以将本学位论文的全部或 部分内容编入有关数据库进行检索和传播，可以采用影印、缩印或扫描等手段保存、汇编、出版本 学位论文。导师K【&年5月/彳日保密的学位论文在解密后适用本授权书。作者签名：乌办公勿伤年e月分日学位论文原创性声明本人郑重声明：所呈交的学位论文，是我个人在导师指导下进行的研究工作及取得的研究成果。 论文中除了特别加以标注和致谢的地方外，不包含任何其他个人或集体已经公开发表或撰写过的研 究成果。其他同志对本研究的启发和所做的贡献均已在论文中作了明确的声明并表示了谢意。本人学位论文与资料若有不实，愿意承担一切相关的法律责任。学位论文作者

3、签名：X年r月目日中图分类号：0241.6单位代码：1032()UDC：51密级：无硕士学位论文四元数保结构算法及其在彩色人脸识别中的应用Quaternion structure-preserving method and itsapplication in color face recognition作者马茹茹导 师贾志刚副教授申请学位理学硕士培养单位江苏师范大学学科专业计算数学研究方向科学与工程计算答辩委员会主席魏木生教授评阅人致谢时光如梭，日子总是在经意与不经意间匆匆流过.转眼间，我的三年硕士学 习生涯即将结束.回首这三年的每一天，仿佛所有的一切都历历在目.很幸运，在 人生最美的年华来到

4、师大,遇到了可亲、可爱、可敬的你们.值此论文完成之际, 我真诚地想向所有曾引导过、鼓励过、帮助过我的师长、朋友、同学表达我的 感激之情.首先，我要感谢我的导师贾志刚老师.谢谢老师毫无保留地传授给我们 专业知识,谢谢老师引领我们走进学术的殿堂，谢谢老师以其宽阔的胸襟包容我 们犯过的许许多多的小错误！由于研一期间老师在英国做访问学者，第一年的学 习都是通过网络授课给我们，这种严谨的教学态度以及对学生的负责之心让我们 感动！老师高尚的品质，严谨的治学态度，渊博的知识，求实的工作作风对我产生 深刻的影响，这些必将令我受益终生！同时，我要衷心地感谢庞宏奎老师、胡秀 玲老师、凌思涛老师和李建波老师,谢谢他

5、们在讨论班上的悉心教导！感谢孙振 威，马龙翠,车增艳，庞烁烁等亲爱的同学们对我学习和生活上的帮助!感谢一直 陪伴我的小伙伴-吴杰同学，谢谢你和我一同经历所有的时光.感谢师大数学与统 计学院为我们提供优良的实验环境，保证本论文的顺利完成.特别地，我要感谢我 的父母和家人，谢谢他们含辛茹苦地把我养育成人,谢谢他们一直以来对我的理 解和支持,谢谢他们在我最失落的时候一次又一次的鼓励!家,将永远是我心灵停 靠的港湾!最后,衷心感谢各位专家、教授在百忙之中评阅和参加本论文答辩！四元数保结构算法及其在彩色人脸识别中的应用中文摘要四元数理论是威廉卢云-哈密尔顿于19世纪40-60年代创立，是复数在四 维空间

6、的不可交换延伸，在理论研究和科学与工程计算中都起着重要作用.由于 四元数乘法的不可交换性，四元数运算面临计算速度慢、难于实现大规模计算等 困难.实表示方法是当今求解四元数矩阵计算问题常用的算法之一，但是实表示 矩阵的维数是四元数矩阵维数的4倍,这使得存储量和运算量过大.本文针对四元 数特征值问题提出一类保结构方法一一保结构Jacobi算法,算法的基本思想是: 利用四元数实表示矩阵的JRS-对称性,构造保结构的Jacobi旋转变换，每步变换 将一个8 x 8的子矩阵约化为对角矩阵，通过循环迭代逐步将整个矩阵的非对角 元素化为零，得到一个仍具有JRS对称性的对角矩阵.该算法充分利用实运算实 现He

7、rmitian四元数矩阵的对角化，具有计算精度高、计算速度快、存储量利运 算量小的特点.针对四元数实表示阵我们设计了最大元保结构Jacobi算法和行循 环保结构Jacobi算法，并给出了收敛性分析.数值实验表明,保结构Jacobi算法与 经典的Jacobi方法相比在运算时间、存储量及计算精度上都有明显的优势.最后, 我们利用保结构Jacobi算法计算彩色人脸的特征脸，并且利用保结构投影处理了 带不同表情、模糊、噪声等干扰因素的彩色人脸照片的识别问题.关键词：四元数实表示矩阵;四元数特征值问题;保结构Jacobi方法;彩色 人脸识别.Quaternion structure-preserving

8、 method and itsapplication in color face recognitionAbstractThe theoiy of quaternion was carved by William Rowan Hamilton in the 1840s to 1860s, is a non-commutative extension in four-dimensional space of complex field, plays an important role in theoretical research and scientific engineering compu

9、ting problem. Since the multiplication of quaternion is non-commutative, quaternion operations faced difficulties in slow calculation speed, is hard to compute with large scale quaternion matrix. The real representation method is one of the most commonly algorithms to solve quaternion matrix problem

10、, however, the dimension of real representation matrix is four times as quaternion matrix, it makes storage space and computation so expensive. In this paper, for quatcrnionic eigenvalue problem, we proposed a kind of structure-preserving methodstruclure-preservdng Jacobi method, The basic ideaof th

11、e algorithm is: lake full advantage of JRS-symmetry of real representation matrix, tectonic structure-preserving Jacobi rotation, transform a 8 x 8 submatrix into a diagonal form each step, make off-diagonal elements of whole matrix to zeros step by step while every iteration, get a diagonal matrix

12、with the same JRS-symmetry finally. The algorithm takes full use of real operation to realize diagonalization in Hermitian quaternion matrix, with high accuracy, fast speed, small amount of storage space and computation. For quaternion real representation matrix, we proposed maximum element structur

13、e-preserving Jacobi method and row cyclic structure-preserving Jacobi method, and presented the convergency analysis. Numerical experiments illustrate that, structure-preserving Jacobi method is belter than classical Jacobi method in operation time、storage space and compulation intensive readings. F

14、inally, we apply structurepreserving Jacobi method to compute color faces eigcnfaces, and deal with different expressionss blurrings、noises etcetera factors by using structure-preserving projection for color faces recognition.Key Words: Quatemion real representation matrix; Quaternion eigenvalue pro

15、blem; Structure-preserving Jacobi method; Color face recognition.目录 TOC o 1-5 h z 摘要1 HYPERLINK l bookmark34 o Current Document 目录TV HYPERLINK l bookmark37 o Current Document 第一章引言1 HYPERLINK l bookmark40 o Current Document 1.1哈密尔顿系统的辛算法1 HYPERLINK l bookmark43 o Current Document 1.2四元数保结构算法4 HYPERL

16、INK l bookmark52 o Current Document 1.2.1正交JRS辛矩阵6 HYPERLINK l bookmark55 o Current Document 1.2.2保结构三对角化 6 HYPERLINK l bookmark61 o Current Document 1.2.3右特征值问题7 HYPERLINK l bookmark64 o Current Document 1.3人脸识别简介8 HYPERLINK l bookmark67 o Current Document 1.3.1人脸识别的发展 8 HYPERLINK l bookmark70 o Cu

17、rrent Document 1.3.2人脸识别方法9 HYPERLINK l bookmark80 o Current Document 1.3.3彩色人脸识别简介12 HYPERLINK l bookmark83 o Current Document 1.4本文的结构安排13 HYPERLINK l bookmark86 o Current Document 第二章 四元数代数及四元数彩色图像模型 15 HYPERLINK l bookmark89 o Current Document 2.1四元数代数15 HYPERLINK l bookmark92 o Current Document

18、 2.1.1四元数的定义和性质15 HYPERLINK l bookmark95 o Current Document 2.1.2四元数的运算 17 HYPERLINK l bookmark98 o Current Document 2.2四元数矩阵的特征值及特征向量19 HYPERLINK l bookmark104 o Current Document 2.3彩色图像的四元数表示21 HYPERLINK l bookmark110 o Current Document 第三章四元数保结构Jacobi方法25 HYPERLINK l bookmark113 o Current Documen

19、t 3.1经典Jacobi方法 25 HYPERLINK l bookmark125 o Current Document 3.2保结构Jacobi旋转变换27 HYPERLINK l bookmark128 o Current Document 3.3 8x8的对称实表示矩阵分解3() HYPERLINK l bookmark143 o Current Document 3.4最大元保结构Jacobi方法33 HYPERLINK l bookmark148 o Current Document 3.5行循环保结构Jacobi方法37 HYPERLINK l bookmark151 o Cur

20、rent Document 3.6数值算例39 HYPERLINK l bookmark157 o Current Document 第四章彩色人脸识别45 HYPERLINK l bookmark160 o Current Document 4.1四元数主成分分析法 45 HYPERLINK l bookmark166 o Current Document 4.2彩色特征脸的选取48 HYPERLINK l bookmark169 o Current Document 4.3带模糊与噪声人脸图像识别50 HYPERLINK l bookmark172 o Current Document 参

21、考文献55 HYPERLINK l bookmark198 o Current Document 作者简历59 HYPERLINK l bookmark13 o Current Document 学位论文独创性声明60 HYPERLINK l bookmark204 o Current Document 学位论文数据集61第一章引言本章简单介绍保结构算法和人脸识别发展及其方法，并列出本文整体的结构 安排.1.1哈密尔顿系统的辛算法在过去的几百年里,学者们对普通常微分方程数值解的理论研究已经是非常 成熟了，并且把相对应的许多方法编制成了商业软件投入使用.但是，研究人员 注意到一般的数值方法用在源

22、于理论物理、分子动力学、天文学及相关领域具 有一定结构的动力系统的时候并不总能够使人非常满意.后来学者们发现，能够 反映原问题几何结构特征的数值算法常常可以产生效果更好的定性行为,而且能 够进行稳定的长时间数值模拟工作.我们把这类能够有效保持原系统结构特征的 数值算法称为保结构算法，其中包括辛结构算法、保酉结构算法、四元数保结构 算法、保李-泊松结构算法、保首次积分算法、保泊松结构算法、保接触结构算 法、保体积算法等1.保结构算法的使用已经渗透到了自然科学的各个领域 中，比如流体动力学、天体力学、形态理论、等离子体物理、粒子加速器、分子 动力学、气象学和反应扩散方程等.考虑如下哈密尔顿系统(z

23、 =z=q)Tek,z(tQ)=免其中=表示动量坐标=知妍.扇,表示位置坐标, H(z)称为哈密尔顿函数.此处的矩阵J表示一个反对称矩阵，其定义如下J = Ik ,满足JT = J-1 = -J G 2kx2k.0Ik是一个k X k的单位矩阵.系统(1.1)具有两个重要的性质:能量守恒和辛结构.在欧氏空间W1中我们可定义下面的欧氏内积,即为一个对称非退化双线性函数：G 涉）：=?涉=其中，/为欧氏内积的度量矩阵.倘若我们把度量矩阵7换成一个反对称矩阵J,那 么就得到了一新的非退化双线性函数k可（3） ：= 5 =工（3料一 Q+戒）,e酣七i=l其满足E（。刃=所以它是反对称的.我们把称为辛

24、内积，具有 辛内积空间IT则称为辛空间.欧氏内积我们可以理解为其是对长度的度量，而辛 内积是对面积的度量.特别地，当* = 1时,，（C）则表示R2中由向量心分别 投影到。敞+该坐标平面上张成的平行四边形有向面积之和.在几何上,因为我 们要求R非退化，所以辛空间是偶数维的.b可以表示为笊=dz A Jdz,所以3也 称为辛2-形式.定义1.1 一般地，对于线性映射次：R2fc 一股2七倘若它能保持2一形式以也就 是满足中（吨,物）略,那么则称评为辛变换.线性映射冰一定程度上可以等同于一个矩阵&因此,相应地辛变换具有辛 矩阵的概念.我们称矩阵0是辛的，当且仅当可（3/，渺戒）=0（,心）或等价地

25、，疽=J.明显地,J自身就是一个标准的辛矩阵.对于非线性映射，我们利用其Jacobi矩阵来定义它的辛性.定义1.2如果非线性映射的Jacobi矩阵#(如)是辛矩阵,即等式，(pq)Tjg(p,q) = J 或可(mg)Cg(饱，)涉)=可(,涉)对所有的p,qe U及3 e略成立，则我们称可微映射g : /(c R2fe) 一 A是辛 变换.根据以上定义中的表达式,我们可以证明det成(p,q) = 1,这里隐含着辛映 射是保体积的.下述定理揭示了哈密尔顿系统的相流是辛变换的,并且我们可以 进一步证明，即对任一个光滑的一阶方程组;另=/(X),其相流是辛变换的充要条 件是该方程是局部的哈密尔顿

26、系统.定理L1假设，定义在区域U 上的哈密尔顿函数H(z)是二次连续可微 的，那么对任意给定的时刻,哈密尔顿系统的精确解流由为一个辛变换，也就是 满足令T保)=/换句话说,务可：=$(dz A J)dz = 0,也就是微分2-形式，为守恒的.定理1.2自治哈密尔顿系统(1.1)的能量沿着解的轨道可以保持不变，即H(次)=H.)=常量,上述两个定理给出了哈密尔顿系统的最基本同时也是非常重要的特性.就启 发我们来构造出恰当数值算法使得相对应的性质能够得以保持.然而，我国院士 冯康的一个学生葛忠和美国著名教授Marsden合作于1988年给出了一个否定性 的结论,即对于一般的哈密尔顿系统，构造出既保

27、能量又保辛结构的数值方法 是不可能的.Chartier等学者则利用B级数理论更进一步地证明了：对足够光滑 的任一个哈密尔顿系统,唯一的既保能量又保辛结构的数值方法.下面给出保能量算法和辛算法的定义,其叙述如下：定义1.3我们把数值方法备+i =抛(讯)用到哈密尔顿系统(1.1)中，如果相应 的数值解满足= H(Znh那么则称此方法为保能量算法.定义1.4如果我们把数值方法=如(.端)用到哈密尔顿系统(1.1)满足那么则称它为辛算法.虽然辛算法不能够严格地保能量，但是我们根据向后误差分析发现，它能够 保持一个扰动的哈密尔顿量.特别地，对于可积或者是近可积的哈密尔顿系统，辛 算法可以具有良好的长时

28、间的数值行为，例如对大多数首次积分能够长时间近似 保持、数值解的误差呈线性增长特性、指数在长时间内可以近似保持不变环面 等等.但是，传统意义上的非辛算法往往误差呈二次增长特性，并且其能量是耗散 的.经过大量的数值试验和理论推理证明，我们都发现辛算法具有独特的长时间 跟踪能力和计算稳定性，并且能够正确地反映原系统的定性行为特性，因此也就 还原了原系统的初始面貌.通过大量的国内外学者们的辛勤工作，辛算法亦或称为辛格式的构造理论巳 经变得相当丰富了,例如有组合分裂法、变分积分法、辛RK法、辛分块RK法、 生成函数法等.1.2四元数保结构算法四元数保结构算法是计算四元数实表示矩阵特征值和特征向量问题的

29、一种 新算法.我们知道，求解矩阵特征值问题有很多种方法.我们熟知的非对称特征值 问题的计算方法有:慕法、反慕法、QR方法等;对称特征值向题的计算方法有对 称QR方法、Jacobi方法、二分法、分而治之法等.在第二章的介绍中我们知道, 四元数实表示阵是实对称矩阵,所以就可以用对称特征值问题的计算方法来求解 该矩阵的特征值和特征向量.在本文中我们提出了四元数保结构Jacobi方法.具 体的算法过程会在第三章重点介绍.另外，这也是给我们今后的工作延伸的一个 方向.既然我们已经讨论了保结构Jacobi方法，并发现它在计算对称和JRS-对称 矩阵的特征值问题有较好的结果，那么我们可以继续研究保结构QR方

30、法、保结 构二分法等，总结归纳出各个算法的优缺点，以便将算法应用于实际生活中.一个丸x n Hermitian四元数矩阵的实表示矩阵有两种结构:对称和JRS-对 称,较于一个一般饥x 4n实矩阵的16疽个参数，它只依赖于2n2 - 2个参数.因 此，一个保结构的算法比一般方法更有效.现在我们定义对于任意一个四元数矩阵Q 瞄心”,若。=Q0 + Q12 + Q2J + 。3加QM印叩=0, 1,2,3）,则其实表示矩阵2Q1-Qi(1.2)Qi-Q2Q-2事实上，如果我们定义Q = 4 + Bj,其中刀=+ QE = Q2 + Q展G饥那么有Tq = U*A + B：j0000A-Bj00000

31、 +跖0000 A- Bj(1.3)这里我们定义的U为酉四元数矩阵:In-jin kinU = 1/2InjnkinLljhiUnInInUn kln其中In为单位矩阵.适合矩阵结构的一个理想方法满足:强向后稳定，即计算解是持有相同结构的扰动矩阵的精确解；可靠，在被考虑的矩阵类中能够计算所有的特征值问题；要求做。（静）次浮点型运算，最好是更少于一个有竞争力的一般方法.1.2.1正交JRS辛矩阵每一个正交的JRS辛矩阵IF有以下块状结构:Wo1的应3顷2形)昭_W1_匹1帝2-匹3_附2甲。W =,w0,wuw2,w3en.(1.4)两种类型的正交矩阵都有以上的形式，其中一个是一个4n x饥的广

32、义辛Given 旋转：Gih-1000000000000COSQQ00COS 0200cos%00COS。3000Il-i000000000000A-l000000000 COS00cos a。00COS 任300COS000000Il-i000000000000A-i000000cos a i00cos a 300cos a。00cosa?000000000金1000000000000L_i000cosotg00COSQi00COS Q200cos a。000000000000h-i(1.5)其中 1 / 0t l/(r+ / + 丁 2)；elset = -l/(-T + a/1 + T2

33、)；endc = 1 / vl + 伐;s = tc;elsec = 1;s = 0;end算法3.2经典Jacobi方法给定对称阵A 股和容许界tol 0,本算法用覆盖A,其中V为正交阵且 offiyAV) tol选择(岔时使厨,讪=max详顶|q育|(c, s) = sym.schur2(4, j, k);v = vj&5end3.2保结构Jacob漩转变换设A = aij是4?i x 4实表示阵,形如式(1.2). Jacobi方法的思想就是逐步 地减小nn n。/(A)=(冈 -浦1 =( 篮)1，(3.1)。=1i=l j=l,详i即矩阵A的非对角元素的“范数” .注意A三a为不变量

34、,off(A)2为a减 去更新后矩阵的对角元素的模的平方和.实现此的基本工具就是如下的平面旋转 变换：定义3.1对一个8x8矩阵coSo$0勺0$2S-20()SS100$3S300c0$00S30_S1G(饱 q,。)=_$20-SoCo$30-5100_S10_S3CoSQ0S2-510 S30-SoCo宠00_$30Sl0_S2CoSo0Si0$20Co(3.2)其中诺+ + S： += 1.我们称为广义Jacobi旋转变换.定理3.1已知由定义3.1定义的G是正交JRS-辛阵，即它是正交阵且满足 GJnGT = J” GRnGT = Rn GSnGT = Sn.对于 JRS-对称实表示

35、阵 A = q,e 卧心4伫有GtAG仍是JRS-对称的，即在G的作用下不改变实表示阵的JRS-对 称性.证G是一个4n x 4n实正交JRS-辛矩阵，且刀是JRS-对称的那么，由第二 章第二节JRS-对称的定义我们可得Jn(GTAG)J = (JnGT)A(JnGTf=GTG(JnGT)A(JnGT)TGTG=GT(GJnGT)A(GJnGT)TG=GT(JnAJ)G=GtAG,Rn(GTAG)R = (RnGT)A(RnGT)T=GTG(RnGT)A(RnGT)TGTG=GT(GRnGT)A(GRnGT)TG=GT(RnAR)G=GtAG,Sn(GTAG)S = (5nGr)A(SnGT)

36、T=GTG(SnGT)A(SnGT)TGTG=GT(GSnGT)A(GSnGT)TG=GT(SnAS)G=GtAG.即G以G仍是JRS-对称的.四元数保结构Jacobi特征值求解的过程基本步骤：(1)选择满足1 Pqn的指标(0 g); (2)计算一组余弦正弦值(c0, Sq. su %第)使得33 2 0 B-BBR21 3 OB BBB-F 彪S0-BB12 13 BOW-B Go G2 Gi G3 G2 Go Gi GiG3GoG2-G3GiG2GorGoG2GiG3G2GoG3GiGiG3GoGa-G3GiG2GoAo A2 Ai 一由 A。 A3 _Ai _人3 Ao 一人3 Al

37、-A2(3.3)为对角阵.即31 =B2 = B3 =0 00 0其中，Oii 。00Q& =,& =。0 Qjja1 00,Bo =HP0bqq _,& =0,-3 =03Q，20_。，30cq %_s。 Co0 Si0 s20 S.3,Gi 一,G2 =,G3 一Si 0s2 0_ S3 0 _以及计算3 = GTAG并覆盖原矩阵A.注意到矩阵3和刀除了，和q这两行两列外均相等.而且有Frobenius范数 正交变换下不变,我们有4 扇 |2 + 8|qo|2 + 8|q+ 8|a2|2 + 8|a3|2 = 4|6PP|2 + 4bqq2这样4off(B)=仔临一工沥睨k=l4=MIIf

38、 一 E nakk + 4|Gpp+ 4|q”|2 - 4|物2 - 4|5W|2k=l=。/(人)-8|q，o, 8|(ii|2 - 8|。2, - 8厨|2在此意义下,每一次Jacobi循环后矩阵同更靠近对角形.在我们讨论怎样选取下标对(mg)前,让我们看看(mg)子问题的实际计算 过程.3.3 8x8的对称实表示矩阵分解定理3.2设a3o0a2oao OFOalopsao alo闻o如叩0 oo芯aooa2 020 如叩 Follo opd 为 Oa3o 7 2 1 3 ao以。ao Yo 妃ftoOa2oaloa3存在一个形如式(3.2)的广义Jacobi旋转变换矩阵G使得GTAG为具

39、有JRS-对 称性的对角矩阵.a3o-alol2oao听 OFOalOF 必ao.? 2 aloasoao%ao ofo做Clooa2 %o&o叱oalo o如0030 aoJFO00 如 Soa2oaloa3S30 70 轮 0约,vELo al(3.5)使得矩阵B的非对角元素为零.对角线元素为:= B33 =方55 = “77 诺姐 + ($3 + s： + s| + s)bjj2勺$0缶)2q)siS 2。0$2饥2cqS3 &3,B-22 B4 = 3(56 B8=CQjj +(S3 + + 崩 + 曷)加+2Co5(A)+ 2。0$1饥 + 2。0$2饥 + 2?0$3如gtag =

40、Sil00000000B-2200000000B3300000000S4400000000B.5500000000B6600000000&00000000&8上述证明也是构造正交JRS-辛矩阵的过程，给出了广义旋转Jacobi矩阵的方法.我们总结8 x 8实表示阵计算如下：算法3.3 8 x 8实表示阵Jacobi旋转过程给定一个4n x 4n对称及JRS-对称实表示阵A和整数p和q,满足1 M P q bqq,如知板,b3）Qi = sqrt（bl + 段 + 辱 + t =（饥0 -饥/）/（2。1）；t = signer）/（t + /l + T）；= 1 / Vl +12；3 =切2；

41、0 =。2；So =例（b（）/Qi）；si =。3也函1）;$2 =。3（切。1）；S3 =。3（枷/。1）；。0 = c0 So； So Co；G2 = 0 $2； s2 0；Gi = 0 Si； Si 0;（7.3 = 0 S3； S3 0;G = G（）G 0,这一算法用GTAG覆盖A,其中G是正交阵且off（GTAG） toliter 100iter = iter + 1;b荻=/!（）（，）;% = WoQj）；场=& （侦）；bi =刀1（槌）；勿=刀2（很）；b：3 = &（很）；G = RQ J R（ba, bjj,? 62）3）Bo = &（R,，,:）;& = &（很B

42、= Bq B2 Bi B3； B-2 Bq B3 B Bi B,313（）B2； B.3 Bi Bo； B = GB;血(皿:)=B(l,2,l:心&(很,:)=B(l?2,2n + 1: 3n); A2(i,j,:) = B(l,2,n + 1 : 2丸);&(底顶,:)=B(l,2,3n + 1 : 4n);Bo =力o(:很);务=Ai(:, i,j);B2 = 02(：,?打)屈=&(:，底顶);B =修0 B? By 氏；-屋 Bq Bl Bl 一 Bq B2； -B3 B- B2 Bq, B = BG;A)(：, )=3(1 ： W2);人 1(：杯)=3(1 项,2 x 2 +

43、1 : 3 x 2); &(:山顶)=3(1 顶,2 + 1 : 2 x 2); &(:/) = 8(1 项,3 x 2 + 1: 4 x 2);E = Aq A2 Al 人3； 山2 A)刀3 力1； 一人1 A3 Aq A2； A3 Al A-2 Aq； z, j = findmaxindE)i = intmod(i,n); j = intmod(j.n);M(k) = off A;论=& + 1;off A = offdiag(A0l A2, &);end while算法3.5求矩阵A模最大元的下标function: p. g = findinaxind(X)SA = sum(A(i).

44、2) sum(diag(A).2)m, n = size(A);T = triu(A, 1);T(：)；tmp, i dx = max( abs (7(:);/;, q = ind2sub(m： n, zt/x);对于最大元保结构Jacobi方法,我们有如下的收敛性定理：定理3.3存在4的特征值的一个排列人 1,人2, ,人,入 1,人2, ；人 1：人2, , ?入72，人1，人2，, , ,；人，使得lim Ak diagX A2? ?人”,入,人2, . ：舄?人 1；人2,Anj Ai,人2, ?人). (3.7) OO证我们首先证明随着迭代次数k的增加,4的非对角“范数off(Ak)

45、趋 于0.由前面的讨论我们知道。打时=。打料_1)2 4|或T)巳(3.8)这里的成一】)|2 = U=2就且或T)是&T的非对角元之中最大的元素.再注 意到。(4一1)2 饥伽 - 1)|唯-叩,(3.9)将(3.9)代入(3.8)式即有。打M。功妇=。”(&1)2(1侦4；1)1=。/(&-1) (1 - 页)，其中N = n(4n 1).由此即知lin以一名= 0.再证存在A的特征值的一个排列 人1，人2, , , , ?人九，人1，人2, , , ,；人n：入1；人2, , ?人n,人1，人2： , , ?人仇，使得lim = A,z 1,. ,n.(3.10)&TOC假定A的互不相同

46、的特征值之间的最小距离为次即6 = min(|/ - A| : A, /z G A(A), A p.任取正数e满足e kQ有off(Ak) e , , ,入n,使得|舄一 Q?)* |山闾 一 D加 |2of/(&o) &必4, = 1,., 4几(3.11)这样的话只要能证明上式蕴含着|人/ 一g?+D| k0有|A.f - a a(2n+p)(2n+p)f tt(3n+p)(3n+p)时切 J a(n+q)(n+qY 戒潟伽+办Q僧%35)，所以只需证明(3.13)式对/ = P,gs + P预+们2“ + p： :2n + g, 3n + p; 3?. + g成立即可.由式(3.4)和式

47、(3.5)可知必+1)= q*o)+c2(2飕)+(礴。)-够)=礁。)+说2 成。)+札 1 F)礴。)=瘤)-飕。).同理可证噌+】)=必)+嶙。).(3-14)这样，对任何W我们有|必7)-月| = |必)-人出+人?-入+秘役叫2 I播-棚-1。*-L)一入计-川。(编)5 - e - e2e,(3.15)这里用到了 t 1.此外，由于M&E) = A(A)和。/7(棚F) &因此 q献1)必须与A的某个特征值之间的距离小于&这样，结合(3.15)式即知(3.13)式 对i = p成立.同样，从(3.14)式出发可以推出(3.13)式对i = q成立.从定理的证明我们可以看出，选择|Z

48、| 1对最大元保结构Jacobi算法的收敛 性起到了至关重要的作用，它保证了迭代产生的每一个对角元将目标一致地趋向 于矩阵A的某一固定的特征值.此外，这一定理的证明亦给出了该算法的收敛速 度的一个粗略估计：1ofJAkf (1 -页)切/(&咒(3.16)这表明最大元保结构Jacobi方法是线性收敛的.但是，实际上其渐近收敛速 度是二次的.更具体一点讲,我们可以证明，存在常数G使得。打(&+N) c -off(Ak)2,(3.17)对充分大的自然数N成立.关于这一结果的证明，可参看Golub和Van Loan合 著的矩阵计算的第八章所引用的参考文献.3.5行循环保结构Jacobi方法最大元保结

49、构Jacobi算法的麻烦之处是每次校正要做。(4冲次运算，但选 取最大元素(p,g)却需要。(16疽)次运算.解决此不平衡的一条途径是将变换的 顺序固定下来.一种合理的选择是逐行对每个非对角元素进行变换.例，若 =4,我们做如下循环:(饱 g) = (1,2),(1,3),(1,4),(2,3),(2,4),(3,4),(1,2),.这种排序格式称为按行循环，它导致如下算法：算法3.6行循环保结构Jacobi算法给定一四元数矩阵A e HIx”和容许界0,这一算法用G1 AG覆盖？1,其中 G 是正交阵且 off(GTAG) tolkiter toliter 100iter = iter +

50、1;for p = 1 : n - 1for g = p + 1 : nbpp = &(00)； bqq = &(g, q); bQ = A0(p, q);bi = & 3，g)；饥=2(p, q);b3 = &(R9)；G = RQ JR(bpp,bo,上,b2,】)3);Bq = AQ(p, q,:); B = &伽,q, :);B2 = A2(p,亦：)；& = &(?,亦：)；B = Bq B-2 Bl B3； B-2 Bq B3 B-y & Bq Bq； B B? Bq B = G 13;&(函,：)=B(l,2,l: n);&(叵亦:)=3(1,2,2 + 1: 3n);A2(p,

51、q,:) = B(l,2,?z + 1 : 2做&(饱亦:)=B(l,2,3n + 1 : 4n);Bq = A)(：, p, g); B=瓦(：血 g)；月2 =人 2(：血 g);& = &(:, p,g);B = Bq B2 Bi B3; -B2 Bq B3 Bi； By 氏 B()% -构 Bi B();B = BG血(:,0 q) = 3(1 项,1,2)3i(：,叵q) = 3(1 项,2 x 2 + 1 : 3 x 2);力2(： p,g) = B(1 : n,2 + 1 : 2 x 2);&(:, p,q) = 5(1 : n,3x 2+ 1 : 4x2); end forof

52、f A = offdiag(A0f A, A2, A3);end forend while行循环保结构Jacobi方法也是二次收敛(见Wilkinson(1962)和Van Kem- per(1966).然而，由于行循环保结构Jacobi算法不需要进行非对角线搜索，它比 最大元保结构Jacobi方法要快得多.在最大元保结构Jacobi算法和行循环保结构Jacobi算法中，除了需要确定保 结构Jacobi旋转变换外,还要在每一次循环后都要给出覆盖后矩阵的非对角元素 之和.通过判断非对角元素之和，我们可以得知最后矩阵是否趋于对角矩阵，从而 判断是否可以使算法停止.下面我们给出四元数矩阵的非对角元素

53、之和的具体算 法.算法3.7四元数矩阵A的非对角元素之和function: offA = offdiag, AbA2, &)a = 4* (norm(Ao, frof) + norm(Ai fro) + normA fro) + norm(A3/ n = size(AQ, 1);for ?： = 1 项offA = sqrt(a 4 * norm(diag(AQ)/ fro,)end for3.6数值算例本文所有数值算例运行处理器为：Intel(R)Core(TM)i5-3230M CPU 2.60GHz,安装内存(RAM): 4.00GB.例1.这个例子中，首先给出一个3 x 3的四元数矩阵

54、，我们希望说明算 法3.4最大元保结构Jacobi和算法3.5行循环保结构Jacobi对于求解该四元数矩阵 的实表示阵特征值问题的有效性.这里，我们选择求解特征值的精度为IO-】。最 高迭代次数为50.测试矩阵为：21 + 2* + 3顶 + 上 4 + 2 + 顶 + 5用Q= l-2i-3j-k53+i+j + k 4-2i-j-5k 3-i-j-k该四元数矩阵的实表示阵为:力()Ai 0：3-A-2&)&_ &A =-& A) &_刀3 -4 A Aq2 1 402 2& =1 5 3? A =-2 0 14 3 3-2 -1 0其中，0310152 =-301,& =-101-1-10

55、-5-10用最大元保结构Jacobi算法我们得到的数值结果为:A)A2-4oAi&-瓦A =_刀3Ao&_ &思其中,-5.3123-0.0000-0.0000& =0.00003.2660-0.0000-0.00000.000012.04630.0027-0.02790.19230.00000.01480.0000-0.17470.00010.0060-0.0100-0.00760.20130.00000.00790.0000-0.2308-0.0131-0.0085-0.1459-0.26080.25270.0000-0.01930.00000.4283-0.3982-0.2260A2 =

56、 10 14Ai = 10 14a3 = 10 15运行时间为0.0227秒.用行循环保结构Jacobi算法我们得到的数值结果为:刀0人2_旦2&人3-Aa2一 &3Ai&其中,-5.3123& =0.0000 0.0000-0.0000 12.0463 0.0000-0.0000-0.0171& = 10项0 2840-0.00560.0000 3.26600.00000.0000-0.2045 -0.0000-0.2328 -0.12370.03510.00000.0000A2 =io-150.4622-0.04470.0000-0.19310.0029-0.0832-0.41760.00

57、00-0.0000人3 =10150.3391-0.2373-0.00000.1109-0.01730.1366运行时间为0.0018秒.例2.这个例子我们想说明对于一幅彩色人脸图像如何用本文中的算法来实 现.选取JSNU彩色人脸库作为测试集，形成16384阶协方差矩阵.类似于例1用算 法3.5计算该四元数矩阵特征值问题.我们所得到的特征向量称为特征脸（图3.1）.图3.1部分彩色特征脸图例例3,针对一 512 X 512维的实表示矩阵九我们分别用最大元保结构Jacobi 算法、行循环保结构Jacobi算法和经典的Jacobi方法来计算该对称正定矩阵的 特征值问题,通过比较迭代次数、非对角线元

58、素的和以及运算时间我们可看出最 大元保结构Jacobi算法和行循环保结构Jacobi算法的优越性，矩阵 A 由 MATLAB 命令 W=rand(n) , A0= (W+Wf )/2, X=rand(n),Al=(X-Xf)/2, Y=rand(n), A2=(Y-Y，)/2, Z=rand(n),A3= (Z-Z)/2, A=A0 A2 Al A3；-A2 AO A3 -Al;-Al -A3 AOA2;-A3. Al -A2 A0,n=128产生，这是一个对称正定矩阵且是JRS- 对称矩阵.迭代精度为IO-%这里我们不设定最大迭代步数，以便比较两种算法 的有效迭代步数.下面三张图分别是经典的

59、Jacobi算法图3.2、最大元保结构 Jacobi算法图3.3和行循环保结构Jacobi算法图3.4所计算的非对角元素和的变化 图.图3.2经典Jacobi方法计算矩阵A的非对角元之和的变化图9000800070006000500040003000200010000图3.3最大元保结构Jacobi算法计算矩阵A的非对角元之和的变化图图3,4行循环保结构Jacobi算法计算矩阵A的非对角元之和的变化图观察三种算法计算矩阵A的非对角元素之和的变化图，我们可以看到, 经典的Jacobi方法所计算的非对角元之和变化幅度比较小,在经过1000步迭 代后的非对角元之和为2.9019 x 104,在经过3

60、个小时的运行后计算精度量级 仍为10气所以我们对经典的Jacobi算法设置了最大迭代步数1000,迭代时间 为59.8386秒.而最大元保结构Jacobi算法在收敛程度上要好于经典的Jacobi方 法,需要说明的是由于图32的横坐标迭代步数较大(大于4万步)，我们用MATLAB 命令set (gca,z xticklabel, x)重新刻画横坐标，运行时间为456.7655秒.而 行循环保结构Jacobi算法仅在第9步迭代后的计算精度为5.6284 x IO27,运行时 间为12.6782秒.第四章彩色人脸识别人脸识别是模式识别研究领域的重要课题，成功的识别算法在很大程度上依 赖于模式特征的提