第8章 矩阵特征值计算.ppt

第8章矩阵特征值计算 8 1特征值性质和估计8 2幂法及反幂法8 3正交变换与矩阵分解8 4QR方法 8 1特征值性质和估计 工程技术中有多种振动问题 如桥梁或建筑物的振动 机械零件 飞机机翼的振动 及一些稳定性分析和相关分析在数学上都可转化为求矩阵特征值与特征向量的问题 8 1 1特征值问题及性质 求A的特征值问题等价于求A的特征方程 设矩阵A Rn n 特征值问题是求 C和非零向量x Rn 使 的根 Ax x 其中 是矩阵A的特征值 x是矩阵A属于特征值 的特征向量 定理1设 为A Rn n的特征值 且Ax x x 0 则有 2 为A I的特征值 且 A I x x 1 c 为的cA特征值 c 0为常数 且 cA x c x 3 k为Ak的特征值 且Akx kx 4 若A为非奇异矩阵 则 1为A 1的特征值 且A 1x 1x 特征值的基本性质 定理2 1 设矩阵A Rn n可对角化 即存在非奇异矩阵P使 的充分必要条件是A具有n个线性无关的特征向量 2 如果矩阵A Rn n有m个 m n 不同的特征值 1 2 m 则对应的特征向量x1 x2 xm线性无关 定理3设A Rn n为对称矩阵 则 3 存在一个正交矩阵P使的 且 1 2 n为A的特征值 而P u1 u2 un 的列向量uj为A的对应于 j的单位特征向量 1 A的特征值均为实数 2 A有n个线性无关的特征向量 定义设A Rn n为对称矩阵 对于任一非零向量x 0 称 为对应于向量x的瑞利 Rayleigh 商 定理4设A Rn n为对称矩阵 其特征值依次记为 1 2 n 则 1 对任何非零x Rn 2 证明只证 1 由于A为实对称矩阵 可将 1 2 n对应的特征向量x1 x2 xn正交规范化 即有 xi xj ij 设x 0为Rn中任一向量 则 于是 瑞利商必位于 n和 1之间 1 成立 8 1 2特征值估计与扰动 定义1设n阶矩阵A aij 令 1 2 集合称为复平面上以aii为圆心 以ri为半径的n阶矩阵A的n个格什戈林 Gerschgorin 圆盘 定理5 Gerschgorin圆盘定理 特别地 如果A的一个圆盘Di是与其它圆盘分离 即孤立圆盘 则Di中精确地包含A的一个特征值 1 设n阶矩阵A aij 则A的每一个特征值必属于下面某个圆盘之中 2 如果A有m个圆盘组成一个连通的并集S 且S与余下n m个圆盘是分离的 则S内恰包含A的m个特征值 或者说A的特征值都在n个圆盘的并集中 证明只就 1 给出证明 Ax x 其中x x1 x2 xn T 0 或 记 考虑Ax x的第k个方程 即 于是 即 设 为A的特征值 A的每一个特征值必位于A的一个圆盘中 并且相应的特征值 一定位于第k个圆盘中 其中k是对应特征向量x绝对值最大的分量的下标 利用相似矩阵性质 有时可以获得A的特征值进一步的估计 即适当选取非奇异对角阵 并做相似变换 适当选取可使某些圆盘半径及连通性发生变化 例1估计矩阵的特征值范围 解矩阵A的3个圆盘为 由圆盘定理可知 A的3个特征值位于3个圆盘的并集中 由于D1是孤立圆盘 所以D1内恰好包含A的一个特征值 1 为实特征值 即 A的其它两个特征值 2 3包含在D2 D3的并集中 取对角阵 做相似变换 矩阵A1的3个圆盘为 3个圆盘都是孤立圆盘 所以 每一个圆盘都包含A的一个特征值 为实特征值 且有估计 定理6 Bauer Fike定理 设 是A I Rn n的一个特征值 且P 1AP D diag 1 2 n 则有 其中 p为矩阵的p范数 p 1 2 证明 A 时显然成立 故只考虑 A 这时D I非奇异 设x是A I对应于 的特征向量 由 A I I x 0左乘P 1可得 而对角矩阵 D I 1的范数为 所以有 这时总有 A 中的一个 取到m值 cond P 是特征值扰动的放大系数 但将A对角化的相似变化矩阵不是唯一的 所以取cond P 的下确界 称为特征值问题的条件数 只要v A 不很大 矩阵微小扰动只带来特征值的微小扰动 但是v A 难以计算 有时只对一个P 用cond P 代替v A 特征值问题的条件数和解线性方程组时的矩阵是两个不同的概念 例 二阶矩阵A diag 1 10 10 有v A 1 但解线性方程组的矩阵条件数cond A 1010 计算矩阵A的特征值 当n 2 3时 可按行列式展开的办法求特征方程p 0的根 但当n较大时 如果按展开行列式的办法 首先求出p 的系数 再求p 的根 工作量就非常大 用这种办法求矩阵特征值是不切实际的 由此需要研究A的特征值及特征向量的数值方法 介绍 计算机上 两类方法 1 幂法及反幂法 迭代法 2 正交相似变换的方法 变化法 8 2幂法及反幂法 幂法与反幂法都是求实矩阵的特征值和特征向量的向量迭代法 幂法是计算矩阵的按模最大的特征值和相应特征向量的一种向量迭代法 特别适用于大型稀疏矩阵 反幂法是计算非奇异 可逆 矩阵按模最小的特征值和相应特征向量的一种向量迭代法 特别是计算海森伯格阵或三对角阵的对应一个给定近似特征值的特征向量的有效方法之一 8 2 1幂法 又称乘幂法 现讨论求 1及x1的方法 设实矩阵A aij 有一个完全的特征向量组 即A有n个线性无关的特征向量 设矩阵A的特征值为 1 2 n 相应的特征向量为x1 x2 xn 已知A的主特征值 1是实根 且满足条件 幂法的基本思想是 任取非零的初始向量v0 由矩阵A构造一向量序列 vk 称为迭代向量 由假设 v0可唯一表示为 于是 其中 由假设故从而 为 1的特征向量 所以当k充分大时 有 则vk为矩阵A的对应特征值 1的一个近似特征向量 用 vk i表示vk的第i个分量 则当k充分大时 有 即为A的主特征值 1的近似值 由于 这种由已知非零向量v0及矩阵A的乘幂Ak构造向量序列 vk 以计算A的主特征值 1及相应特征向量的方法就称为 乘 幂法 幂法的思想 由矩阵A的乘幂Ak与非零向量v0相乘来构造向量序列 vk Akv0 从而计算主特征值 1及其对应的特征向量 的收敛速度由比值 来确定 r越小收敛越快 但当r 1时收敛可能很慢 定理7设A Rn n有n个线性无关的特征向量 主特征值 1满足 1 2 n 则对任何非零向量v0 又设A有n个线性无关的特征向量 1对应的r个线性无关的特征向量为x1 x2 xr 则 如果A的主特征值为实的r重根 即 1 2 r 且 r r 1 n 为 1的特征向量 这说明当A的主特征值是实的重根时 定理7的结论还是正确的 应用幂法计算A的主特征值 1及其对应的特征向量时 如果 1 1 或 1 1 迭代向量vk的各个不等于零的分量将随k 而趋向于无穷 或趋向于零 这样在计算机实现时就可能 溢出 为克服这个缺点 就需要将迭代向量加以规范化 设有一向量v 0 将其规范化得向量为 其中max v 表示v的绝对值最大的第一个分量 下标最小 即如果有 则max v vq 且q为所有绝对值最大的分量中的最小下标 由于 取一规范化初始向量v0 0 构造规范化向量序列为 规范化向量序列收敛到主特征值对应的特征向量 收敛速度由比值r 2 1 确定 同理 定理8设A Rn n有n个线性无关的特征向量 主特征值 1满足 1 2 n 则对任意非零初始向量v0 u0 a1 0 有幂法计算公式为 uk vk 则有 例2用幂法计算矩阵 的主特征值和相应的特征向量 解取v0 u0 1 1 1 T 则 计算结果如下表 这个结果是用8位浮点数字进行运算得到的 uk的分量值是舍入值 于是得到 及相应的特征向量 0 7482 0 6497 1 T 1和相应的特征向量的真值 8位数字 为 例用幂法计算矩阵 的主特征值与其对应的特征向量 解取v0 u0 0 0 1 T 则 直到k 8时的计算结果见下表 从而 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的特征值的某种分布 它是十分有效的 例3设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的特征值都是实数 满足 2 10 式 且 2 n能初步估计出来 我们就能确定P 的近似值 例4用原点平移加速法求例2中矩阵A的主特征值与其对应的特征向量 对B应用幂法 仍取v0 1 1 1 T 则 解取p 0 75 做平移变换B A pI 则 计算结果如下表 由此得B的主特征值为 1 1 7865914 A的主特征值 1为 1 1 0 75 2 5365914 原点位移的加速方法 是一个矩阵变换方法 这种变换容易计算 又不破坏矩阵A的稀疏性 但p的选择依赖对A的特征值分布的大致了解 及相应的特征向量为x1 0 7482 0 6497 1 T 与例2结果比较 上述结果比例2迭代15次还好 若迭代15次 1 1 7865258 相应的 1 2 5365258 就是真值 设A Rn n为对称矩阵 称 为向量x的瑞利商 其中 x x xTx为内积 由定理4知道 实对称矩阵A的特征值 1及 n可用瑞利商的极限值表示 下面我们将瑞利商应用到用幂法计算实对称矩阵A的主特征值的加速上来 2 瑞利商 Rayleigh 加速 定理9设A Rn n为对称矩阵 特征值满足 对应的特征向量xi满足 xi xj ij 单位正交向量 应用幂法公式 2 9 计算A的主特征值 1 则规范化向量uk的瑞利商给出 1的较好的近似值为 由此可见 R uk 比 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 此时反幂法的迭代公式为 反幂法的迭代公式为 任取v0 u0 0 构造 对给定的误差 当 k k 1 时 得 显然 反幂法的收敛速度取决于比值 比值越小 收敛越快 定理10设A Rn n为非奇异矩阵 且有n个线性无关的特征向量 其对应的特征值满足 1 2 n 2 n 0 则对任意非零初始向量u0 an 0 由反幂法计算公式构造的向量序列 vk uk 满足 1 2 收敛速度的比值r n n 1 在反幂法中也可以用原点平移法加速迭代过程 或求其它特征值与其对应的特征向量 如果矩阵 A pI 1存在 显然其特征值为 对应的特征向量仍然是x1 x2 xn 现对矩阵 A pI 1应用幂法 得到反幂法的迭代公式 如果p是A的特征值 j的一个近似值 且设 j与其它特征值是分离的 即 就是说1 j p 是矩阵 A pI 1的主特征值 可用反幂法 2 12 计算特征值 j及特征向量 设A Rn n有n个线性无关的特征向量x1 x2 xn 则 其中 同理可得下面的定理 定理11设A Rn n有n个线性无关的特征向量 矩阵A的特征值及对应的特征向量分别记为 i及xi i 1 2 n 而p为 j的近似值 A pI 1存在 且 1 2 则对任意非零初始向量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 进行迭代 反幂法计算公式 1 分解计算P A pI LU 且保留L U及P信息 2 反幂迭代法 1 解Uv1 1 1 1 T求v1 2 解k 2 3 解Lyk Puk 1求yk 解Uvk yk求vk k max vk 计算uk vk k 例5用反幂法求 的对应于计算特征值 1 2679 精确特征值为 的特征向量 用5位浮点数进行计算 解用部分选主元的三角分解将A pI 其中p 1 2679 分解为 P A pI LU 其中 由Uv1 1 1 1 T 得 由LUv2 Pu1 得 3对应的特征向量是 由此看出u2是x3的相当好的近似 特征值 3的真值为 8 3正交变换与矩阵分解 8 3 1豪斯霍尔德变换 定义设向量w Rn 且wTw 1 称矩阵 为初等反射矩阵 这个矩阵也称为豪斯霍尔德变换 如果记w w1 w2 wn 则有 定理12设有初等反射矩阵 其中H I 2wwT 其中wTw 1 则 1 H是对称矩阵 即HT H 2 H是正交矩阵 即H 1 H 3 设A为对称矩阵 那么A1 H 1AH HAH即亦是对称矩阵 证明只证H的正交性 其中显然 HTH H2 I 2wwT I 2wwT I 4wwT 4w wTw wT I 设向量u 0 则显然 是一个初等反射矩阵 初等反射矩阵的几何意义 考虑以为w法向量且过原点O的超平面S wTx 0 任取向量v Rn 则v x y 其中x S y ST 于是 Hx I 2wwT x x 2wwTx x 从而 Hv x y v 其中v 为v关于平面S的镜面反射 初等反射矩阵在计算上的意义是它能用来约化矩阵 设向量Hx 0 可选择一初等反射矩阵H 使Hx y Hy I 2wwT y y 2wwTy y 定理13设x y为两个不相等的n维向量 x 2 y 2 则存在一个初等反射矩阵H 使Hx y 证明令 则得到一个初等反射矩阵 而且 因为 所以 w是使Hx y成立的唯一长度等于1的向量 不计符号 定理14 约化定理 设x x1 x2 xn T 0 则存在初等反射矩阵H 使Hx e1 其中 证明记y e1 设x y 则有 x 2 y 2 由定理13 存在变换H I 2wwT 使Hx y e1 其中 其中记u x1 x2 xn T 如果 和x1异号 那么计算x1 时有效数字可能损失 我们取 和x1有相同的符号 即取 在计算 时 可能上溢或下溢 为了避免溢出 将x规范化 记u x e1 u1 u2 un T 于是 则有H 使H x e1 其中 例6设x 3 5 1 1 T 则 x 2 6 取 6 可直接验证Hx e1 6 0 0 0 T 8 3 2吉文斯变换 设向量x y R2 则变换 或 是平面上向量的一个旋转变换 其中 为正交变换 R2中变换 其中x x1 x2 xn T y y1 y2 yn T 而 称为Rn中平面 xi xj 的旋转变换 也称为吉文斯变换 P P i j P i j 称为平面旋转矩阵 显然 P P i j 具有性质 1 P与单位阵I只是在 i i i j j i j j 位置元素不一样 其它相同 2 P为正交矩阵 P 1 PT 3 P i j A 左乘 只需计算第i行与第j行元素 即对A aij m n有 其中c cos s sin 4 AP i j 右乘 只需计算第i列与第j列元素 利用平面旋转变换 可得向量x中的指定元素变为零 定理15 约化定理 设x x1 xi xj xn T 其中不全为零 则可选择平面旋转阵P P i j 使 其中 于是 由c s的取法得 证明取c cos xi x i s sin xj x i 由P i j x x x 1 x i x j x n T 利用矩阵乘法 显然有 8 3 3矩阵的QR分解与舒尔分解 定理16设A Rn n非奇异 则存在正交矩阵P 使PA R 其中R为上三角矩阵 证明1 用吉文斯变换给出构造P的方法 1 第1步约化 由设有j j 1 2 n 使aj1 0 则可选择吉文斯变换P 1 j 将aj1处的元素化为零 即若aj1 0 j 2 3 n 则存在P 1 j 使得 可简记为P1A A 2 其中P1 P 1 n P 1 2 2 第k步约化 设上述过程已完成第1步至第k 1步 于是有 由设有j n j k 使ajk k 0 j k 1 n 则可选择吉文斯变换P k j j k 1 n 使 PkA k P k n P k k 1 A k PkPk 1 P1A A k 1 其中Pk P k n P k k 1 3 继续上述约化过程 最后则有 Pn 1 P2P1A A n R 上三角形矩阵 令P Pn 1 P2P1 它是一个正交矩阵 有PA R 2 用豪斯霍尔德变换构造出正交矩阵P 记A 0 A 它的第一列记为 不妨设a1 0 0 找到矩阵 于是 其中 一般地 设 其中D j 1 为j 1阶方阵 其对角线以下元素均为0 矩阵 A j 1 为n j 1阶方阵 设其第一列为a1 j 1 可选择n j 1的豪斯霍尔德矩阵变换 使 根据 Hj构造n n阶的变换矩阵Hj为 于是有 它和A j 1 有类似的形式 只是D j 为j阶方阵 其对角线以下元素是0 这样经过n 1步运算得到 Hn 1 H2H1A A n 1 R 上三角形矩阵 其中P Hn 1 H2H1是一个正交矩阵 从而有PA R 定理17 QR分解定理 设A Rn n为非奇异矩阵 则存在正交矩阵Q与上三角矩阵R 使 A QR 且当R的对角元素为正时 分解是唯一的 证明从定理16可知 只要令Q PT就有A QR 下面证明分解的唯一性 设有两种分解 A Q1R1 Q2R2 其中Q1 Q2为正交矩阵 R1 R2为对角元素均为正的上三角矩阵 则 由假设及对称正定矩阵ATA的楚列斯基 平方根法 分解的唯一性 则得R1 R2 从而可得Q1 Q2 证毕 定理16保证了A可分解为A QR 若A非奇异 则R也非奇异 如果不规定R的对角元素为正 则分解不是唯一的 一般按吉文斯变换或豪斯霍尔德变换方法作出的分解A QR R的对角元不一定是正的 设上三角矩阵R rij 只要令 则为正交矩阵 为对角元是 rii 的上三角矩阵 这样便是符合定理17的唯一QR分解 例7用豪斯霍尔德变换作矩阵A的QR分解 解按 3 2 式找豪斯霍尔德矩阵H1 R3 3 使 则有 再找 H2 R2 2 使 H2 1 44949 3 44949 T 0 T 得 这是一个上三角矩阵 但对角元素皆为负数 只要令D I 则有R H2H1A是对角元素为正的上三角矩阵 取 则得A QR 除了QR分解 矩阵的Schur分解也是重要的工具 它解决矩阵A Rn n可约化到什么程度的问题 对复矩阵A Cn n 则存在酉矩阵U 使UHAU为一个上三角矩阵R 其对角线元素就是A的特征值 A URUH称为A的舒尔分解 对于实矩阵A 其特征值可能有复数 A不能用正交相似变换约化为上三角矩阵 但它可约化为以下形式 其中对角块Rii i 1 2 m 为一阶或二阶方阵 且每个一阶Rii是A的实特征值 每个二阶对角块Rii的两个特征值是A的两个共轭复特征值 定理18 实Schur分解 设A Rn n 则存在正交矩阵Q使 A QRQT称为A的实舒尔分解 可以通过逐次正交变换使A趋于实舒尔型矩阵 以求A的特征值 8 3 4用正交相似变换约化一般实矩阵为上海森伯格矩阵 设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 重复上述过程 则有 定理19 豪斯霍尔德约化矩阵为上海森伯格阵 设A Rn n则存在初等反射矩阵U1 U2 Un 2使 为上海森伯格矩阵 总结上述结论 有 本算法约需要5n3 3次乘法运算 如果要把U0也算出来还需要增加2n3 3次乘法 例8用豪斯霍尔德方法将矩阵 约化为上海森伯格阵 解选取初等反射阵R1使 其中c1 2 4 T 1 计算R1 则有 2 约化计算 令 则得到上海森伯格阵为 如果A是对称的 则H U0TAU0也对称 这时H是一个对称三对角矩阵 定理20 豪斯霍尔德约化对称矩阵为对称三对角矩阵 设A Rn n为对称矩阵 则存在初等反射矩阵U1 U2 Un 2使 为对称三对角矩阵 证明由定理19 存在初等反射矩阵U1 U2 Un 2使为上海森伯格阵 且An 1亦是对称矩阵 因此 An 1为对称三对角矩阵 由上面讨论可知 当A为对称阵时 由Ak Ak 1 AkUkAk一步约化计算中只需计算Rk及RkA22 k Rk 又由于A的对称性 故只需计算RkA22 k Rk的对角线以下元素 注意到 引进记号 则有 对对称阵A用初等反射阵正交相似约化为对角三对角阵大约需要2n3