求特征值的Jacobi方法

山东科学 SHANDONG SCIENCE 第 24 卷 第 6 期 2011 年 12 月出版 Vol． 24 No． 6 Dec． 2011 收稿日期:2011-10-13 基金项目:山东省高等学校教学改革研究重点项目(2009044) ;山东省高等学校教学改革研究项目(2009262) 作者简介:于正文(1960 －) ，男，副教授，研究方向为科学计算。Email:yuzhengwen@ sdjzu． edu． cn DOI:10． 3976 / j． issn． 1002 － 4026． 2011． 06． 006 求特征值的 Jacobi方法 于正文1，尹庆莉2 (1．山东建筑大学理学院，山东 济南 250101;2．山东建筑大学实验与设备处，山东 济南 250101) 摘要:讨论了求实对称矩阵的特征值的经典 Jacobi方法，通过一系列的正交相似变换将实对称矩阵化为对角矩阵，从而 求出全部特征值和相应的特征向量。文中给出所有正交变换的计算公式，并用 MATLAB编程实现，为实际问题的计算提 供了简单实用的计算工具。 关键词:特征值;特征向量;Jacobi方法 中图分类号:O24 文献标识码:A 文章编号:1002-4026(2011)06-0019-03 Jacobi method for eigenvalues calculation YU Zheng-wen1，YIN Qing-li2 (1． School of Sciences，Shandong Jianzhu University，Jinan 250101，China; 2． Department of Laboratory and Facility Management，Shandong Jianzhu University，Jinan 250101，China) Abstract ∶ This paper addresses classical Jacobi method of eigenvalues calculation of a real symmetric matrix． It converts a real symmetric matrix into a diagonal matrix by a series of orthogonal similarity transformations，and then derives all eigenvalues and their corresponding eigenvectors． The paper presents the formulas of all orthogonal transformations， which are implemented by MATLAB programming． This provides a simple and practical calculation tool for the computation of practical problems． Key words ∶ eigenvalue;eigenvector;Jacobi method 特征值问题的解法有两类，即迭代法和变换法。前者算法简单易于编程，但仅适用于求模最大或最小特 征值及相应的特征向量;后者是将矩阵通过相似变换化矩阵为特殊形式，如将对称矩阵三对角化，以及将一 般矩阵化为上 Heisenberg矩阵，然后再求这两种特殊形式矩阵的特征值。作相似变换后所得矩阵的特征值 与原矩阵的特征值相同，但特征向量与矩阵的特征向量一般不同，需要记忆所作过的变换。因此，变换法的 算法和编程较复杂，但其优点是可求全部或部分特征值及相应的特征向量。经典的 Jacobi 方法是一种求实 对称矩阵的全部特征值及相应特征向量的方法，它是通过一系列的正交相似变换，将实对称矩阵化为对角矩 阵，并进而求出全部特征值及相应特征向量。本文将讨论经典 Jacobi方法的算法，并用 MATLAB编程实现。 1 经典的 Jacobi方法 设 A∈Rn × n是实对称矩阵，则正交矩阵 R∈Rn × n使得 山 东 科 学 2011 年 R － 1AR = diag(λ1，λ2，…，λn) ， (1) 其中{λ i} n i = 1是 A的 n个实特征值，而矩阵 R的第 j列 r(j)为 A的相应于特征值 λ j 的特征向量 u j(j = 1， 2，…，n)。 但是，由于通过实际计算不可能准确地给出正交矩阵 R，因而需要借助平面旋转矩阵，对 A进行一系列 的正交相似变换，以实现矩阵 A的对角化。 定义 称矩阵 Gij(θ)= 1  1 cos θ sin θ 1  1 － sin θ cos θ 1   1 第 i行 第 j行 (2) 为 n维空间 Rn 中(i，j)平面上的旋转矩阵或 Givens矩阵。 经典 Jacobi方法的基本思想是:对矩阵 A0 = A，构造一系列的平面旋转矩阵 G1，G2，…，Gk，并计算: Ak =GkAk － 1G T k (k = 1，2，…)。 (3) 则当 k→∞时，有 Ak→diag(λ1，λ2，…，λn)。注意到矩阵的 F-范数‖A‖F 是正交不变量，因此要使 A近似对 角化，就应使对角元的平方和尽可能的增大。假设已对矩阵 A 进行了 k － 1 次变换，得到矩阵 Ak － 1 = (aij (k － 1))n × n，则将平面旋转矩阵 Gk =Gpq(θ)代入(3)中，并比较其左右两端元素得 a(k)pp = a (k － 1) pp cos 2 θ + 2a(k － 1)pq cos θsin θ + a (k － 1) qq sin 2 θ a(k)pq = a (k － 1) pp sin 2 θ － 2a(k － 1)pq cos θsin θ + a (k － 1) qq cos 2 θ a(k)pq = 1 2 (a (k － 1) qq － a (k － 1) pp )sin2 θ + a (k － 1) pq cos2      θ (4) a(k)pi = a (k) ip = a (k － 1) ip cos θ + a (k － 1) iq sin θ (i≠p，q) a(k)qi = a (k) iq = － a (k － 1) ip sin θ + a (k － 1) iq cos θ (i≠p，q{ ) (5) a(k)ij = a (k － 1) ij (i，j≠p，q) (6) 变化了的 p，q行(列)的元素间满足关系式 a(k)[ ]pi 2 + a(k)[ ]qi 2 = a(k － 1)[ ]pi 2 + a(k － 1)[ ]qi 2 (i≠p，q)。 根据矩阵的 F-范数的正交不变性，有 a(k)[ ]pp 2 + a(k)[ ]qq 2 + 2 a(k)[ ]pq 2 = a(k － 1)[ ]pp 2 + a(k － 1)[ ]qq 2 + 2 a(k － 1)[ ]pq 2。 Ak 与 Ak － 1相比，其对角元的模平方和最多增加 2 a (k － 1)[ ]pq 2，而且只有当 θ满足 a(k)pq = 0 时才能达到。为了使 Ak 的对角元的模平方和尽可能地增加，在第 k步计算过程中，要选 Ak － 1的模最大的非对角元 a (k － 1) p0q0 ，并构造 适当的平面旋转矩阵 Gp0q0(θ)= Gk，使得 Ak = GkAk － 1G T k 的元(p0，q0) ，a (k) p0q0 = 1 2 (a (k － 1) q0q0 － a (k － 1) p0p0 )sin 2θ + a(k － 1)p0q0 cos 2θ = 0。如上选法可以保证 A的对角元的模平方和以最快速度增加。由于 A的元素 aij在某一步变 换化为 0 后，在以后的变换中有可能变为非 0。因此，一般很难在经有限次变换后将给定的矩阵 A化为对角 02 第 6 期 于正文，等:求特征值的 Jacobi方法 矩阵，但有如下的收敛性定理。 定理 1 按经典 Jacobi方法生成的矩阵序列{Ak}的非对角元素均收敛于 0，即limk→∞Ak = diag(λ1，λ2，…，λn)。 证明:见文献［1］。 2 算法 本算法的主要部分及计算公式为 (1)选定 Ak － 1 =(a (k － 1) ij )中的元素 a (k － 1) p0q0 使得 | a (k － 1) p0q0 | = max1≤i≤n 1≤j≤i － 1 | a(k － 1)ij |。 (2)确定旋转矩阵 Gk 使 θ 满足 tan 2θ = 2a(k － 1)p0q0 a(k － 1)p0p0 － a (k － 1) q0q0 (| θ | ≤ π4 )。当 a (k － 1) p0p0 = a (k － 1) q0q0 时，取 θ = π4 sign(a (k － 1) p0q0 )。否则 Gk 的 sin θ，cos θ分别取为 cos θ = 1 / 1 + tan 2槡 θ = 1 / 1 + t槡 2，sin θ = tcos θ。 其中 t = tan θ = sign (c)·(－ | c | + 1 + c槡 2) ，c = cot 2θ = a(k － 1)p0p0 － a (k － 1) q0q0 2a(k － 1)p0q0 。 (3)按式(4)～(6)计算 Ak。 (4)特征向量的计算 注意到 Ak 的计算格式，有 Ak = GkGk － 1…G1AG T 1G T 2…G T k，记 R T k = GkGk － 1…G1，并设 R0 = I，则有 ARk = RkAk，Rk = Rk － 1·G T k(7)。在计算 Ak 的同时，可按(7)式计算生成序列{Rk}，并且当 k→∞时，Rk 的各列就 是近似特征向量。Rk =(r k ij)的计算格式为: r(k)ip = r (k) ip cos θ + r (k) iq sin θ (i = 1，2，…，n) r(k)iq = － r (k － 1) ip sin θ + r (k － 1) iq cos θ (i = 1，2，…，n) r(k)ij = r (k － 1) ij (i = 1，2，…，n;j≠p，q { ) (8) 设 A∈Rn × n是实对称矩阵，则其全部特征值的近似值及相应近似特征向量可由 MATLAB程序求得。 限于篇幅，源程序略。 3 计算实例 例 求实对称矩阵 A = 10 8 12 － 9 7 8 － 7 0 11 5 12 0 － 6 9 12 － 9 11 9 － 3 5  7 5 12 5 － 9 的全部特征值及对应的特征向量。 解:在 MATLAB命令窗口输入 A，执行程序 Jacobi后，有下面的结果。 ＞＞ A =［10，8，12，－ 9，7;8，－ 7，0，11，5;12，0，－ 6，9，12;－ 9，11，9，－ 3，5;7，5，12，5，－ 9］ ＞＞［D，V］= Jacobi(A) D的对角线即为近似的特征值 V的列向量为相应的近似特征向量 D = 23． 7229 0 0 0 0 0 － 5． 5496 0 0 0 0 0 － 27． 1214 0 0 0 0 0 11． 0372 0 0 0 0 0 － 17． 0891 (下转第 100 页) 12 山 东 科 学 2011 年 ［3］韩延彬，李金屏． 一种基于局部特征的人脸识别方法［J］． 济南大学学报:自然科学版，2007，21(1) :56 － 59． ［4］唐资娜． 基于肤色和 Hausdorff距离的人脸检测系统［D］． 南京:南京航空航天大学，2007． ［5］龚理专，王巍．基于肤色信息和主分量分析的人脸实时检测系统［J］．计算技术与自动化，2005，1(24) :92 － 94． ［6］邢果． 彩色图像中人脸检测与识别方法的研究［D］． 郑州:中国人民解放军信息工程大学，2007． ［7］彭进业，余卞章，王大凯，等．多尺度对称变换及其应用于定位人脸特征点［J］．电子学报，2003，30(3) :363 － 366． ［8］严超，朱光大．人脸特征的定位与提取［J］．中国图编图形学报，1998，3(5) :375 － 380． (上接第 21 页) V = 0． 7174 0． 2172 － 0． 3703 － 0． 4960 － 0． 2347 0． 2917 0． 7459 0． 4861 0． 3021 0． 1762 0． 4828 － 0． 5524 0． 5547 0． 2015 － 0． 3367 0． 1145 0． 0736 － 0． 5224 0． 7596 － 0． 3627 0． 3926 － 0． 2931 － 0． 2143 0． 2124 0． 8179 4 结论 经典 Jacobi方法的计算工作量非常大，利用此程序容易计算实对称矩阵的全部特征值及对应的特征向 量，经过验证程序正确无误，给实际问题的应用和教学提供了简单实用的计算工具。 参考文献: ［1］翟瑞彩，谢伟松． 数值分析［M］． 天津:天津大学出版社，2000:87 － 95． ［2］威尔金森． 代数特征值问题［M］．北京:科学出版社，2001:275 － 293． ［3］HEATH M T． Scientific Computing:An Introductory Survey［M］． Second Edition，北京:清华大学出版社(影印版) ，2001:169 － 203． ［4］MATEJA J． Quadratic convergence of scaled matrices in Jacobi method［J］． Numer Math，2000，87 (1) :171 － 199． ［5］张智星． Matlab程序设计与应用［M］．北京:清华大学出版社，2002． ［6］HARRELLII E M． Eigenvectors and Eigenvalues［EB /OL］． (2007 － 12 － 10) ［2011 － 09 － 20］． http:/ /www． mathphysics． com / calc /eigen． html． 001