矩阵特征值问题计算.ppt

第8章矩阵特征问题的计算 8 1引言8 2幂法及反幂法8 3豪斯霍尔德方法8 4QR方法 8 1引言 工程技术中有多种振动问题 如桥梁或建筑物的振动 机械机件 飞机机翼的振动 及一些稳定性分析和相关分析在数学上都可转化为求矩阵特征值与特征向量的问题 下面先复习一些矩阵的特征值和特征向量的基础知识 定义1 已知n阶矩阵A aij 则 称为A的特征多项式 一般有n个根 实的或复的 复根按重数计算 称为A的特征值 用 A 表示A的所有特征值的集合 A的特征方程 设 为A的特征值 相应的齐次方程组 注 当A为实矩阵时 0为实系数n次代数方程 其复根是共轭成对出现 的非零解x称为矩阵A的对应于 的特征向量 例1求A的特征值及特征向量 其中 解矩阵A的特征方程为 求得矩阵A的特征值为 对应于各特征值矩阵A的特征向量分别为 定理1设 为A Rn n的特征值 且Ax x x 0 则有 p为A pI的特征值 即 A pI x p x c 为的cA特征值 c 0为常数 下面叙述有关特征值的一些结论 k为Ak的特征值 即Akx kx 设A为非奇异矩阵 那么 0 且 1为A 1的特征值 即A 1x 1x 定理2设 i i 1 2 n 为n阶矩阵A aij 的特征值 则有 称为A的迹 定理3设A Rn n 则有 定理4设A为分块上三角矩阵 即 其中每个对角块Aii均为方阵 则 定理5设A与B为相似矩阵 即存在非奇异矩阵P使B P 1AP 则 定理5说明 一个矩阵A经过相似变换 其特征值不变 一个亏损矩阵是一个没有足够特征向量的矩阵 亏损矩阵在理论上和计算上都存在困难 1 A与B有相同的特征值 2 如果y是B的特征向量 则Py是A的特征向量 定义2如果实矩阵A有一个重数为k的特征值 且对应于 的A的线性无关的特征向量个数 k 则A称为亏损矩阵 定理6 1 A Rn n可对角化 即存在非奇异矩阵P使 的充分必要条件是A具有n个线性无关的特征向量 2 如果A Rn n有m个 m n 不同的特征值 1 2 m 则对应的特征向量x1 x2 xm线性无关 定理7 对称矩阵的正交约化 设A Rn n为对称矩阵 则 存在一个正交矩阵P使的 且 1 2 n为A的特征值 而P u1 u2 un 列向量uj为A的对应于 j的单位特征向量 A的特征值均为实数 A有n个线性无关的特征向量 定义3设n阶矩阵A aij 令 下面讨论矩阵特征值界的估计 集合称为复平面上以aii为圆心 以ri为半径的n阶矩阵A的n个Gerschgorin圆盘 定理8 Gerschgorin圆盘定理 特别地 如果A的一个圆盘Di是与其它圆盘分离 即孤立圆盘 则Di中精确地包含A的一个特征值 设n阶矩阵A aij 则A的每一个特征值必属于下面某个圆盘之中 如果A有m个圆盘组成一个连通的并集S 且S与余下n m个圆盘是分离的 则S内恰包含A的m个特征值 或者说A的特征值都在n个圆盘的并集中 证明只就 给出证明 设 为A的特征值 即 Ax x 其中x x1 x2 xn T 0 或 记 考虑Ax x的第k个方程 即 于是 即 这说明 A的每一个特征值必位于A的一个圆盘中 并且相应的特征值 一定位于第k个圆盘中 其中k是对应特征向量x绝对值最大的分量的下标 利用相似矩阵性质 有时可以获得A的特征值进一步的估计 即适当选取非奇异对角阵 并做相似变换 适当选取可使某些圆盘半径及连通性发生变化 例2估计矩阵A的特征值范围 其中 解矩阵A的3个圆盘为 由定理8 可知A的3个特征值位于3个圆盘的并集中 由于D1是孤立圆盘 所以D1内恰好包含A的一个特征值 1 为实特征值 即 A的其它两个特征值 2 3包含在D2 D3的并集中 现在取对角阵 做相似变换 矩阵A1的3个圆盘为 显然 3个圆盘都是孤立圆盘 所以 每一个圆盘都包含A的一个特征值 为实特征值 且有估计 当A为实矩阵 如果限制用正交相似变换 由于A有复的特征值 A不能用正交相似变换约化为上三角阵 用正交相似变换能约化到什么程度呢 定理9 Schur定理 设A Rn n 则存在酉矩阵U使 其中rii i 1 2 n 为A的特征值 下面给出理论上有关通过酉相似变换及正交变换可以约化一般矩阵A到什么程度的问题 其中Rii i 1 2 m 为一阶或二阶方阵 且每个一阶Rii是A的实特征值 每个二阶对角块Rii的两个特征值是A的两个共轭复特征值 定理10 实Schur分解 设A Rn n 则存在正交矩阵Q使 定义4设A Rn n为对称矩阵 对于任一非零向量x 称 我们转向实Schur型的实际计算 为对应于向量x的瑞利 Rayleigh 商 定理11设A Rn n为对称矩阵 其特征值次序记为 1 2 n 则 1 对任何非零x Rn 2 3 证明只证1 关于2 3自己作练习 由于A为实对称矩阵 可将 1 2 n对应的特征向量x1 x2 xn正交规范化 则有 xi xj ij 设xi 0为Rn中任一向量 则有展开式 于是 从而1成立 结论1说明瑞利商必位于 n和 1之间 关于计算矩阵A的特征值问题 当n 2 3时 我们还可按行列式展开的办法求 0的根 但当n较大时 如果按展开行列式的办法 首先求出 的系数 再求 的根 工作量就非常大 用这种办法求矩阵的特征值是不切实际的 由此需要研究求A的特征值即特征向量的数值解法 本章将介绍一些计算机上常用的两类方法 一类是幂法及反幂法 迭代法 另一类是正交相似变换的方法 变换法 幂法与反幂法都是求实矩阵的特征值和特征向量的向量迭代法 所不同的是幂法是计算矩阵的主特征值 矩阵按模最大的特征值称为主特征值 其模就是该矩阵的谱半径 和相应特征向量的一种向量迭代法 而反幂法则是计算非奇异 可逆 矩阵按模最小的特征值和相应特征向量的一种向量迭代法 下面分别介绍幂法与反幂法 8 2幂法及反幂法 现讨论求 1及x1的方法 设实矩阵A aij 有一个完全的特征向量组 即A有n个线性无关的特征向量 设矩阵A的特征值为 1 2 n 相应的特征向量为x1 x2 xn 已知A的主特征值 1是实根 且满足条件 8 2 1幂法 又称乘幂法 幂法的基本思想是 任取非零的初始向量v0 由矩阵A构造一向量序列 vk 称为迭代向量 由假设 v0可唯一表示为 于是 其中 由假设故从而 为 1的特征向量 所以当k充分大时 有 即为矩阵A的对应特征值 1的一个近似特征向量 用 vk i表示vk的第i个分量 则当k充分大时 有 即为A的主特征值 1的近似值 由于 这种由已知非零向量v0及矩阵A的乘幂Ak构造向量序列 vk 以计算A的主特征值 1 2 7 及相应特征向量 2 5 的方法就称为幂法 迭代公式实质上是由矩阵A的乘幂Ak与非零向量v0相乘来构造向量序列 vk Akv0 从而计算主特征值 1及其对应的特征向量 这就是幂法的思想 的收敛速度由比值 来确定 r越小收敛越快 但当r 1时收敛可能很慢 定理12设A Rn n有n个线性无关的特征向量 主特征值 1满足条件 1 2 n 则对任何非零向量v0 a1 0 幂法的算式成立 又设A有n个线性无关的特征向量 1对应的r个线性无关的特征向量为x1 x2 xr 则由 2 2 式有 如果A的主特征值为实的重根 即 1 2 r 且 r r 1 n 为A的特征向量 这说明当A的主特征值是实的重根时 定理5的结论还是正确的 应用幂法计算A的主特征值 1及其对应的特征向量时 如果 1 1 或 1 1 迭代向量vk的各个不等于零的分量将随k 而趋向于无穷 或趋向于零 这样在计算机实现时就可能 溢出 为克服这个缺点 就需要将迭代向量加以规范化 设有一向量v 0 将其规范化得向量为 其中max v 表示v的绝对值最大的分量 即如果有 则max v vq 且q为所有绝对值最大的分量中的最小下标 在定理12的条件下幂法可这样进行 任取一初始向量v0 0 a1 0 构造规范化向量序列为 由 2 3 式 收敛速度由比值r 2 1 确定 总结上述结论 有 同理 可得到 定理13设A Rn n有n个线性无关的特征向量 主特征值 1满足 1 2 n 则对任意非零初始向量v0 u0 a1 0 有幂法计算公式为 例1用幂法计算矩阵 的主特征值与其对应的特征向量 解取v0 u0 0 0 1 T 则 直到k 8时的计算结果见下表 从而 见书p303 例3 8 2 2幂法的加速方法 1 原点平移法 由前面讨论知道 应用幂法计算A的主特征值的收敛速度主要由比值r 2 1 来决定 但当r接近于1时 收敛可能很慢 这时 一个补救办法是采用加速收敛的方法 其中p为参数 设A的特征值为 i 则对矩阵B的特征值为 i p 而且A B的特征向量相同 引进矩阵B A pI 如果要计算A的主特征值 1 只要选择合适的数p 使 1 p为矩阵B A pI的主特征值 且 那么 对矩阵B A pI应用幂法求其主特征值 1 p 收敛速度将会加快 这种通过求B A pI的主特征值和特征向量 而得到A的主特征值和特征向量的方法叫原点平移法 对于A的特征值的某种分布 它是十分有效的 例4设A R4 4有特征值 比值r 2 1 0 9 做变换 B A 12I p 12 则B的特征值为 应用幂法计算B的主特征值 1的收敛速度的比值为 虽然常常能够选择有利的p值 使幂法得到加速 但设计一个自动选择适当参数p的过程是困难的 下面考虑当A的特征值是实数时 怎样选择p使采用幂法计算 1得到加速 且使收敛速度的比值 设A的特征值都是实数 且满足 则对实数p 使矩阵A pI的主特征值为 1 p或 n p时 当我们计算 1及x1时 首先应选取p使 显然 当 2 p n p 时 即P 2 n 2 P 时 为最小值 这时收敛速度的比值为 当希望计算 n时 应选取p 1 n 1 2 P 使得应用幂法计算 n得到加速 当A的特征值都是实数 满足 且 1 2能初步估计出来 我们就能确定P 的近似值 例2用原点平移加速法求例1中矩阵A的主特征值与其对应的特征向量 对B应用幂法 仍取v0 0 0 1 T 则 解取p 2 5 做平移变换B A pI 则 迭代5步的计算结果见下表 可得到B的主特征值 1 13 5000 主特征向量v1 0 5 1 0 0 7500 T 因此 A的主特征值为 1 1 p 11 0000 主特征向量仍为x1 0 5 1 0 7500 T 原点位移的加速方法 是一个矩阵变换方法 这种变换容易计算 又不破坏矩阵A的稀疏性 但p的选择依赖对A的特征值分布的大致了解 见书p306 例5 设A Rn n为对称矩阵 称 为向量x的瑞利商 其中 x x xTx为内积 由定理11知道 实对称矩阵A的特征值 1及 n可用瑞利商的极限值表示 下面我们将瑞利商应用到用幂法计算实对称矩阵A的主特征值的加速上来 2 瑞利商 Rayleigh 加速 定理14设A Rn n为对称矩阵 特征值满足 对应的特征向量xi满足 xi xj ij 单位正交向量 应用幂法公式 2 9 计算A的主特征值 1 则规范化向量uk的瑞利商给出 1的较好的近似值为 由此可见 R xk 比 k更快的收敛于 1 证明由 2 8 式及 得 幂法的瑞利商加速迭代公式可以写为 其中A为n阶实对称矩阵 对给定的误差限 当 k k 1 时 取近似值 8 2 3反幂法 反幂法是用于求非奇异矩阵A的按模最小的特征值和对应特征向量的方法 而结合原点平移法的反幂法则可以求矩阵A的任何一个具有先了解的特征值和对应的特征向量 设矩阵A非奇异 其特征值 i i 1 2 n 满足 其相应的特征向量x1 x2 xn线性无关 则A 1的特征值为1 i 对应的特征向量仍为xi i 1 2 n 此时 A 1的特征值满足 因此 对A 1应用幂法 可求出其主特征值1 n k和特征向量xn uk 从而求得A的按模最小特征值 n 1 k和对应的特征向量xn uk 这种求A 1的方法就称为反幂法 为了避免求A 1 可通过解线性方程组Avk uk 1得到vk 采用LU分解法 即先对A进行LU分解A LU 此时反幂法的迭代公式为 反幂法的迭代公式为 对给定的误差 当 k k 1 时 得 显然 反幂法的收敛速度取决于比值 比值越小 收敛越快 定理15设A Rn n为非奇异矩阵 且有n个线性无关的特征向量 其对应的特征值满足 1 2 n 2 n 0 则对任意非零初始向量u0 an 0 由反幂法计算公式构造的向量序列 vk uk 满足 在反幂法中也可以用原点平移法加速迭代过程 或求其它特征值与其对应的特征向量 如果矩阵 A pI 1存在 显然其特征值为 对应的特征向量仍然是x1 x2 xn 现对矩阵 A pI 1应用幂法 得到反幂法的迭代公式 如果p是A的特征值 j的一个近似值 且设 j与其它特征值是分离的 即 就是说1 j p 是矩阵 A pI 1的主特征值 可用反幂法 2 12 计算特征值及特征向量 设A Rn n有n个线性无关的特征向量x1 x2 xn 则 其中 同理可得 定理16设A Rn n有n个线性无关的特征向量 矩阵A的特征值及对应的特征向量分别记为 i及xi i 1 2 n 而p为 j的近似值 A pI 1存在 且 则对任意非零初始向量u0 aj 0 由反幂法计算公式 2 12 构造的向量序列 vk uk 满足 且收敛速度为 由该定理知 对A pI 其中p j 应用反幂法 可用来计算特征向量xj 只要选择p是 j的一个较好的近似且特征值分离情况较好 一般r很小 常常只要迭代一二次就可完成特征向量的计算 反幂法迭代公式中的vk是通过解方程组 求得的 为了节省工作量 可以先将A pI进行三角分解 于是求vk相对于解两个三角形方程组 实验表明 按下述方法选择u0是较好的 选u0使 用回代求解 2 13 即得v1 然后再按公式 2 12 进行迭代 反幂法计算公式见书p311 前面已提到 见书p311 例6 8 3豪斯霍尔德方法 1 用初等反射矩阵作正交相似变换约化一般实矩阵A为上海森伯格矩阵 8 3 1引言 本节讨论两个问题 2 用初等反射矩阵作正交相似变换约化对称矩阵A为对称三对角矩阵 于是 求原矩阵特征值问题 就转化为求上海森伯格矩阵或对称三对角矩阵的特征值问题 8 3 2用正交相似变换约化一般实矩阵为上海森伯格矩阵 设A Rn n 下面来说明 可选择初等反射矩阵U1 U2 Un 2使A经正交相似变换约化为一个上海森伯格阵 1 设 其中c1 a21 an1 T Rn 1 不妨设c1 0 否则这一步不需要约化 于是 可选择初等反射阵使 其中 令 则 其中 2 第k步约化 重复上述过程 设对A已完成第1步 第k 1步正交相似变换 即有 或 且 其中为k阶上海森伯格阵 设ck 0 于是可选择初等反射阵Rk使其中 Rk计算公式为 令 则 其中为k 1阶上海森伯格阵 第k步约化只需计算及 当A为对称矩阵时 只需要计算 3 重复上述过程 则有 定理17 豪斯霍尔德约化矩阵为上海森伯格阵 设A Rn n则存在初等反射矩阵U1 U2 Un 2使 为上海森伯格矩阵 总结上述结论 有 算法1 豪斯霍尔德约化矩阵为上海森伯格阵 设A Rn n 本算法计算U0TAU0 H 上海森伯格型 其中U0 U1U2 Un 2为初等反射阵的乘积 1 U0 I 2 对于k 1 2 n 2 1 计算初等反射阵Rk使 本算法约需要5n3 3次乘法运算 要明显形成U0还需要附加2n3 3次乘法 2 约化计算 3 U0 U0Uk 例7用豪斯霍尔德方法将矩阵 约化为上海森伯格阵 解选取初等反射阵R1使 其中c1 2 4 T 1 计算R1 则有 2 约化计算 则得到上海森伯格阵为 8 3 3用正交相似变换约化对称矩阵为对称三对角矩阵 定理18 豪斯霍尔德约化对称矩阵为对称三对角矩阵 设A Rn n为对称矩阵 则存在初等反射矩阵U1 U2 Un 2使 为对称三对角矩阵 证明由定理17 存在初等反射矩阵U1 U2 Un 2使为上海森伯格阵 且An 1亦是对称阵 因此 An 1为对称三对角阵 由上面讨论可知 当A为对称阵时 由Ak Ak 1 AkUkAk一步约化计算中只需计算Rk及RkA22 k Rk 又由于A的对称性 故只需计算RkA22 k Rk的对角线以下元素 注意到 引进记号 则有 对对称阵A用初等反射阵正交相似约化为对角三对角阵大约需要2n3 3次乘法 用正交矩阵进行相似约化有一些特点 如构造的 Uk容易求逆 且Uk的元素数量级不大 这个算法是十分稳定的 算法2见书p318 8 4QR方法 8 4 1QR算法 Rutishauser 1958 利用矩阵的三角分解提出了计算矩阵特征值的LR算法 Francis 1961 1962 利用矩阵的QR分解建立了计算矩阵特征值的QR算法 QR方法是一种变换方法 是计算一般矩阵 中小型矩阵 全部特征值问题的最有效方法之一 1 上海森伯格矩阵的全部特征值问题 2 计算对称三对角矩阵的全部特征值问题 目前QR方法主要用来计算 且QR方法具有收敛快 算法稳定等特点 对于一般矩阵A Rn n 或对称矩阵 则首先用豪斯霍尔德方法将A化为上海森伯格阵B 或对称三对角阵 然后再用QR方法计算B的全部特征值 设A Rn n 对A进行QR分解 即 A QR 其中R为上三角阵 Q为正交阵 于是可得到一新矩阵 B RQ QTAQ 显然 B是由A经过正交相似变换得到 因此B与A的特征值相同 再对B进行QR分解 又可得一新矩阵