07学年论文 宫琳

本科学年论文(设计)题目：对于大规模矩阵矩阵特征值解法初步探索学院统计与数理学院专业信息与计算科学班级07级信息与计算科学二班学号20070534213姓名宫琳指导教师张爱英山东财政学院教务处制二Ｏ一一年摘要：现实生活中，很多问题不同系数都可以归结为大规模矩阵问题，这为处理解决问题提供了方法，但是大规模矩阵往往过于复杂，难于计算和处理，这就需要寻找一种可行的分类解决方法。

在矩阵解法中，矩阵特征值是矩阵内容的核心内容，本文将对于矩阵的特征值解法进行讨论。

对于乘幂法，反幂法，子空间代换法的讨论进行初探，并且对于其中的使用解决领域，精度等进行讨论。

关键字：矩阵特征值 乘幂法 反幂法 子空间代换法引言：在科学技术的应用领域中，许多问题都归为求解一个相应的特征系统。

特征系统里面的诸多因素通过矩阵的形式表现出来，如动力学系统和结构系统中的振动问题，求系统的频率与振型；物理学中的某些临界值的确定等等。

但是一般系统内部因素诸多，数值很多，大规模矩阵情况很复杂，计算繁琐，需要根据不同类型进行矩阵求解。

矩阵解法中利用矩阵特征值进行求解，是一个比较普遍的方法。

这在工程计算领域非常常见，如量子物理中的Kohn-Sham 方程求解化为哈密顿矩阵某些关键特征值的对策；财务风险中的决策树法以财务变量为判别点建立决策树模型，以最低误判成本为标准对样本公司进行分类，进行矩阵运算；会计中所有者权益变动表利用矩阵求解等。

本论文希望通过对大规模矩阵的主特征值解决的思想和方法进行探讨，例如乘幂法，反幂法和子空间代换法。

通过对不同类型的计算方法的结果和适用类型的对比和探讨，以及对于不同类型解决方法的优化等方面的研究，加深对矩阵实际的应用的理解，学以致用。

计算矩阵的时候会引入矩阵的核心内容—特征值λ，下面为相关介绍。

设A 为n 阶方阵，n n ij R a A ⨯∈=)(，若)0(≠∈x R x n ，有数λ使Ax= λx则称λ为A 的特征值，x 为相应于λ的特征向量。

因此，特征问题的求解包括两方面： 1．求特征值λ，满足0)det()(=-=I A λλϕ 2．求特征向量)0(≠∈x R x n ，满足齐方程组0)(=-x I A λ称ϕ（λ）为A 的特征多项式，它是关于λ的n 次代数方程。

关于特征值的求解，一般用乘幂法，反幂法以及子空间代换法，但是不同方法使用的情况不同，精度不同，过程不同，因而对于不同的需求，计算机上机实现选择方式也不同，本文将对其进行阐述。

关于矩阵特征值的解法(一) 乘幂法与反幂法 在实际工程应用中，如大型结构的振动系统中，往往要计算振动系统的最低频率（或前几个最低频率）及相应的振型，相应的数学问题便为求解矩阵的按模最大或前几个按模最大特征值及相应的特征向量问题，或称为求主特征值问题。

1.1 乘幂法 介绍：乘幂法是用于求大型稀疏矩阵的主特征值的迭代方法，其特点是公式简单，易于上机实现。

基本思想：任取一个非零向量，由已知的矩阵A 的乘幂构造一个向量序列，通过不断迭代取得主特征值及相关特征向量。

乘幂法的公式推算为： 设A ∈R n ⨯n ，取初始向量x (0)∈R n ，令x (1) = Ax (0)，x (2) = Ax (1)，…，一般有)1()(-=k k Ax x（1）形成迭代向量序列{x (k)}。

由递推公式，有)0()1(2)2()()(x A x A Ax A x k k k k ====-- （2）这表明x (k )是用A 的k 次幂左乘x (0)得到的，因此称此方法为乘幂法， {x (k )}称为迭代序列。

例1 用规范化乘幂法计算矩阵A 的主特征值及相应特征向量⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---=20101350144A解：首先，通过矩阵的基本运算，不难得到A 的特征值λ1 = 6，λ2 = 3，λ3 = 2其次，运用乘幂法运算，对结果进行对比不妨取初始值x (0) = (1, 1, 1)T ，用规范化乘幂法公式计算1)max()0(=xT x x y )1,1,1()max()0()0()0(==T Ay x )1,8,10()0()1(==其它结果见表1（表中的向量均为转置向量）。

表1 k max(y (k )) x (k ) = y (k )/max(x (k )) x (k +1) = Ay (k )0 1 (1, 1, 1) (10, 8, 1) 1 10 (1, 0.8, 0.1) (7.2, 5.4, -0.8) 2 7.2 (1, 0.75, -0.111111) (6.5, 4.75, -1.222222)3 6.57 (1, 0.730769, -0.203704) (6.230766, 4.499997, -1.407408) 46.230766(1, 0.722222, -0.225880)(6.111108, 4.388886, -1.1451767)5 6.111108 (1, 0.718182, -0.237561) (6.054548, 4.336336, -1.475122)6 6.054548 (1, 0.716216, -0.243639) (6.027024, 4.310808, -1.487278)7 6.027024 (1, 0.715247, -0.246768) (6.013458, 4.298211, -1.483536)8 6.013458 (1, 0.714765, -0.248366) (6.00671, 4.291945, -1.496732)9 6.00671 (1, 0.714525, -0.249177) (6.00335, 4.28825, -1.496354) 10 6.00335 (1, 0.714405, -0.249586) (6.00167, 4.287265, -1.499172) 11 6.00167 (1, 0.714345, -0.239792) (6.00083, 4.286485, -1.499584)12 6.00083 (1, 0.714315, -0.249896)取max(x (12)) = 6.00083作为主特征值λ1的近似值，与真值λ1 = 6相比，有较好的近似程度，相应于λ1的特征向量的近似值取为y (2) = (1, 0.714315, -0.249896)T 。

1.2 乘幂法的加速及降价介绍：乘幂法计算A 的主特征值及相应的特征向量的时候，如果主特征值绝对值大于1，迭代向量的各个不等于零的分量都会随k →∞而趋于无穷，这样会造成计算机的“溢出”。

为了克服这个问题，利用向量的长度与方向无关这一性质，将迭代向量的长度规范化以改进幂法。

当λi ( i = 1, 2, …, n )为矩阵A ∈ R n ⨯n 的n 个特征值，且n λλλ≥≥> 21时，乘幂法的收敛速度由11λλ=r 决定，r << 1收敛得快。

因此为提高收敛速度或改善r ≈ 1的状况，可以采取原点移位的方法，改变原矩阵A 的状态。

取λ0（常数），用矩阵B = A - λ0I 来代替A 进行乘幂迭代。

设μi (i = 1, 2, …, n )为矩阵B 的特征值，则B 与A 特征值之间应有关系式：0λλμ-=i i (i = 1, 2, …, n )且若v i 是A 相应于λi 的特征向量，则v i 亦是μi 的特征值，即对i = 1, 2, …, ni i i i i i v v Av v I A Bv )()(000λλλλ-=-=-=因此，对任取x (0)∈R n ，关于矩阵B 的乘幂公式可为⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛--+-=⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛+=-==∑∑==j j nj j kj kj nj j k k k k v v v v x I A x B x 2010*******111)0(0)0()()()(λλλλααλλμμααμλ为加快收敛速度，适当选择参数λ0，使kj n j 01020max )(λλλλλω--=≤≤达到最小值。

如当λi (i = 1, 2, …, n )为实数，且λ1>λ2 ≥…≥λn 时，取)(212*0n λλλ+=则*0λ为ω(λ0)的极小值点。

这时122122122*01*02221212121λλλλλλλλλλλλλλλλλ<---=----=--n n n n原点移位法是一个矩阵变换过程，变换简单且不破坏原矩阵的稀疏性。

但由于预先不知道特征值的分布，所以应用起来有一定困难，通常对特征值的分布有一个大略估计，设定一个参数λ0进行试算，当所取λ0对迭代有明显加速效应以后再进行确定计算。

例2 计算A 的主特征值⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=0.225.05.025.00.10.15.00.10.1A解 先用规范化乘幂法计算，得表5.2表5.2k y (k ) = x (k )/max(x (k )) λ1 ≈ max (x (k ))0 (1, 1, 1)1 (0.9091, 0.8182, 1)T 2.75 19 (0.7482, 0.6497, 1)T 2.5365374 20(0.7482, 0.6497, 1)T2.5365323主特征值λ1及特征向量v 1为（8位有效数字） λ1 = 2.5362258 v 1 = (0.74822116, 0.64966116, 1)T 而用规范化乘幂法计算的相应近似值为： λ1 ≈ max (x (20)) = 2.5365323 v 1 ≈ y (20) = (0.7482, 0.6497, 1)T 如果采用原点位移的加速法求解，取λ0 = 0.75，矩阵b = A - λ0I⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=25.125.05.025.025.015.0125.0B对矩阵B 应用规范化乘幂法公式见表5.3。

表5.3k y (k ) = x (k )/max(x (k )) λ1 ≈ max (x (k ))0 (1, 1, 1)T9 (0.7483, 0.6497, 1)T 1.7866587 10(0.7483, 0.6497, 1)T1.7865914可见5365914.2011≈+=λμλ此结果与未加速的规范化乘幂法公式计算结果相比，收敛速度要快得多。

2、 反幂法 设A ∈R n ⨯n 可逆，则无零特征值，由)0(≠=x x Ax λ有x x A λ11=- 即若λ为矩阵A 的特征值，则λ1必为矩阵A -1的特征值，且特征向量相同。

如果A 的n 个特征值λi (i = 1, 2, …, n )为n λλλ≥≥≥ 21则A -1的n 个特征值更为11111λλλ≥≥≥- n n因此，若乘幂法可求A 的主特征值λ1，则用A -1做乘幂矩阵，由乘幂迭代格式)(1)1(k k x A x -+=便可求出A -1的按模最大特征值nλ1，取倒数，即为矩阵A 的按模最小特征值。