特征值和Jacobi方法

第 9章 矩 阵 特 征 值 问 题 的 数 值 方 法 9.1 特 征 值 与 特 征 向 量 引 言工 程 实 践 中 有 多 种 振 动 问 题 ， 如 桥 梁 或 建筑 物 的 振 动 ， 机 械 机 件 、 飞 机 机 翼 的 振动 ， 工 程 实 践 中 有 多 种 振 动 问 题 ， 如 桥梁 或 建 筑 物 的 振 动 ， 机 械 机 件 、 飞 机 机翼 的 振 动 ， 及 一 些 稳 定 性 分 析 和 相 关 分析 可 转 化 为 求 矩 阵 特 征 值 与 特 征 向 量 的问 题 。 London, England: Millennium (Wobbly) Bridge (1998-2002, Norman Foster and Partners and Arup Associates) I decide that I have to write something today, otherwise I would not know how to speak English here. This is a very quick story about a bridge. London launched three major construction projects to celebrate the arrival of the Millennium. After all, Greenwich (pronounced green-ich) is supposed to be (supposed to be?!) where the prime meridian lies, and the place where the Millennium officially starts in the world. The three projects are the Millennium Dome in North Greenwich, so far the largest single roofed structure in the world, London Eye right across Westminster, which becomes so far the largest observation wheel in the world, and the Millennium Bridge that links Southeast London with St. Pauls Cathedral, which is currentlywell.not swinging any more, it is said. The bridge was designed by Imperial College, a college of my former university. On the very first day that the bridge was open to public, there were simply so many people going there to walk from the south bank to St. Pauls that the weight completely exceeded the architects expectation. The slender steel truss bridge began to vibrate with a million people on there. The opening ceremony ended up in an embarrassing vertigo. Millennium left Londoners a happy adage about swinging bridge, meaning fancy technology that looks good but functions in a funny fashion. Am I using too many Fs here? Or is it simply because my tongue starts to swing in the same direction when I am writing about this wobbly bridge? Next time you visit London, I strongly recommend this place. After all, with a little swing, this is a shortcut to dash into St. Pauls directly from the southeast! xxG T 1ex TG : G oogle Matrix, “the worlds largest matrix computation”. 4,300,000,000 x: PageRank（ 网 页 级 别 ） vector “The $25,000,000,000 Eigenvector”搜 索 引 擎 9.1 特 征 值 与 特 征 向 量 设 A是 n阶 矩 阵 ， x是 非 零 列 向 量 . 如 果 有数 存 在 ， 满 足 ， (1)那 么 ， 称 x是 矩 阵 A关 于 特 征 值 的 特 征 向量 . Ax x 如 果 把 (1)式 右 端 写 为 ， 那 么 (1)式 又 可 写为 : Ix( ) 0I A x | | 0I A 即 11 1 0( ) | | .n nnf I A a a a 记它 是 关 于 参 数 的 n次 多 项 式 ， 称 为 矩 阵 A的 特征 多 项 式 ， 其 中 a0=(-1)n A . (2) 显 然 ， 当 是 A的 一 个 特 征 值 时 ， 它 必 然是 的 根 . 反 之 ， 如 果 是 的 根 ，那 么 齐 次 方 程 组 (2)有 非 零 解 向 量 x， 使 (1)式成 立 . 从 而 ， 是 A的 一 个 特 征 值 . A的 特 征 值 也 称 为 A的 特 征 根 . ( ) 0f ( ) 0f 矩 阵 特 征 值 和 特 征 向 量 有 如 下 主 要 性 质 ： 定 理 9.1.1 n阶 矩 阵 A是 降 秩 矩 阵 的 充 分 必要 条 件 是 A有 零 特 征 值 . 定 理 9.1.2 设 矩 阵 A与 矩 阵 B相 似 ， 那 么 它 们 有 相 同 的 特 征 值 . 定 理 9.1.3 n阶 矩 阵 A与 AT有 相 同 的 特 征 值 . 定 理 9.1.4 设 i j是 n阶 矩 阵 A的 两 个 互 异特 征 值 ， x、 y分 别 是 其 相 应 的 右 特 征 向 量 和 左特 征 向 量 ， 那 么 ， xTy=0 . 9.2 Hermite矩 阵 特 征 值 问 题 设 A为 n阶 矩 阵 ， 其 共 轭 转 置 矩 阵 记 为 AH. 如果 A=AH， 那 么 ， A称 为 Hermite矩 阵 . 9.2.1 Hermite矩 阵 的 有 关 性 质 设 是 Hermite矩 阵 A的 n个 特 征值 .有 以 下 性 质 : 全 是 实 数 .1 2, ,., n 1 2, ,., n 有 相 应 的 n个 线 性 无 关 的 特 征向 量 ， 它 们 可 以 化 为 一 组 标 准 酉 交 的 特 征向 量 组 ,即 1 2, ,., n 1 2, ,., nu u u Hi ju u ij 是 酉 空 间 中 的 一 组 标 准 酉 交 基 . 1 2, ,., nu u u 记 U=( ),它 是 一 个 酉 阵 ， 即UHU=UUH=I， 那 么即 A与 以 为 对 角 元 的 对 角 阵 相 似 .1 2, ,., nu u u 1H nU AU D 1 2, ,., n A为 正 定 矩 阵 的 充 分 必 要 条 件 是 全 为 正 数 . 1 2, ,., n 定 理 9.2.1 设 是 Hermite矩 阵 A的 n个 特 征 值 ， 那 么证 : 1 2, ,., n 2 1max ii nA 2 1n iF iA 2 2 22 ( ) ( ) ( ( )HA A A A A 由 2 1max ii nA 因 此2 2 21( ) ( ) nH iF iA tr A A tr A 又 由 2 1n iF iA 得 设 x是 一 个 非 零 向 量 ， A是 Hermite矩 阵 ，称 为 矩 阵 A关 于 向 量 x的 Rayleigh商 ，记 为 R(x). HHx Axx x定 理 9.2.2 如 果 A的 n个 特 征 值 为 其 相 应 的 标 准 酉 交 的 特 征 向 量 为 那 么 有 1 2 . n 1 2, ,., nu u u1 ( ) nR x 定 理 9.2.3 设 A是 Hermite矩 阵 ， 那 么 10 0min ( ) min ( )k n kk kx C x x C xR x R x 且 且或 9.2.2 极 值 定 理 定 理 9.2.4(极 值 定 理 ) 设 Hermite矩 阵 的 n个特 征 值 为 ， 其 相 应 的 标 准 酉交 特 征 向 量 为 . 用 Ck表 示 酉 空间 Cn中 任 意 的 k维 子 空 间 ， 那 么1 2 . n 1 2, ,., nu u u 1 10 0max min ( )min max ( )kk n k n kk x C xCk C x C xR xR x 且 且或 9.2.3 Hermite矩 阵 特 征 值 问 题 的 性 态 矩 阵 特 征 值 问 题 与 求 解 线 性 方 程 组 问 题 一样 ， 都 存 在 当 矩 阵 A的 原 始 数 据 有 小 变 化 (小 扰动 )时 ， 引 起 特 征 值 问 题 的 变 化 有 大 有 小 的 问 题 ，如 果 引 起 的 变 化 小 ， 称 该 特 征 值 问 题 是 良 态 的 . 反 之 ， 称 为 病 态 的 . 矩 阵 特 征 值 问 题 的 性 态 是 很 复 杂 的 ， 通 常分 别 就 单 个 特 征 值 或 整 体 特 征 值 给 出 条 件 数 进行 分 析 . 对 于 Hermite矩 阵 ， 由 于 其 特 征 值 问 题的 特 殊 性 质 ， 其 特 征 值 都 是 良 态 的 .下 面 先 证 明Hermite矩 阵 特 征 值 的 扰 动 定 理 . 定 理 9.2.5 设 矩 阵 A， E， A+E都 是 n阶 Hermite矩 阵 ， 其 特 征 值 分 别 为 那 么 ,证 设 矩 阵 A关 于 特 征 值 1， 2， ， n 的标 准 酉 交 特 征 向 量 为 u1， u2， ， un， 是 由 ui， ui+1， ， un生 成 的 n-i+1维 子 空间 . 对 中 任 意 非 零 向 量 x,由 极 值 定 理 ,有1 2 n 1 2 n 1 2 n 1i n i i 1n iC 1n iC 1 1 100 0( )maxmax maxn in i n i Hi Hx C x H HH Hx C x x C xx A E xx xx Ax x Exx x x x 且且 且 由 定 理 9.2.3，又 由 定 理 9.2.2， 对 任 意 x0， 有从 而 有另 一 方 面 , A=(A+E)-E. 记 为 矩 阵 -E的 特 征 值 ， 那 么 ，重 复 上 面 的 过 程 ,可 得从 而 有 1 0maxn i H iHx C x x Axx x 且11 0maxn i H nHx C x x Exx x 且 1i i 1 2 n 1i n i 1i i i i n 定 理 9.2.5通 常 又 称 为 Hermite矩 阵 特 征 值的 扰 动 定 理 定 理 9.2.6 设 矩 阵 A和 A=A+E都 是 n阶 Hermite矩 阵 ， 其 特 征 值 分 别 为 和 ， 那 么 1 2 n 1 2 n 22 2i iE E 9.3 Jacobi方 法 理 论 上 ， 实 对 称 矩 阵 A正 交 相 似 于 以 A的特 征 值 为 对 角 元 的 对 角 阵 . 问 题 是 如 何构 造 这 样 的 正 交 矩 阵 呢 ? Jacobi方 法 就 是通 过 构 造 特 殊 的 正 交 矩 阵 序 列 ， 通 过 相似 变 换 使 A的 非 对 角 线 元 素 逐 次 零 化 来 实现 对 角 化 的 . 9.3.1 平 面 旋 转 矩 阵 与 相 似 约 化先 看 一 个 简 单 的 例 子 . 设 是 二 阶 实 对 称 矩 阵 ， 即 a21=a12，其 特 征 值 为 1， 2. 令 使 得 记 容 易 验 证 BT=B, 且11 1221 22a aA a a cos sinsin cosR 1 1TR AR 11 1221 22T b bB R AR b b 2 211 11 12 22 2 212 21 22 11 122 222 11 12 22cos 2 sin cos sin( )sin cos (cos sin )sin 2 sin cos cosb a a ab b a a ab a a a 解 之 得 :当 时当 时 可 选 取 11 22a a 1211 2221arctan2 aa a 11 22a a 4 为 使 RTAR为 对 角 阵 ， 要 求 b12=b21=0 一 般 的 n阶 平 面 旋 转 矩 阵 21arctan2 pq pp qqaa a 4 9.3.2 经 典 的 Jacobi方 法 设 A是 实 对 称 矩 阵 ， 记 A1=A.Jacobi方法 的 基 本 思 想 是 用 迭 代 格 式 Ak+1=QTkAkQk , k=1,2, 构 造 一 个 相 似 矩 阵 序 列 ， 使 Ak 收敛 于 一 个 对 角 阵 . 其 中 Qk为 平 面 旋 转 矩 阵 ，其 旋 转 角 k由 使 Ak的 绝 对 值 最 大 元a(k)pq=a(k)qp=0 或 按 列 依 次 使 A的 非 对 角 元 零 化 来 确 定 . 定 理 9.3.1 设 A是 n阶 实 对 称 矩 阵 ， 那么 由 Jacobi方 法 产 生 的 相 似 矩 阵 序 列 Ak的 非 对 角 元 收 敛 于 0. 也 就 是 说 ， Ak 收 敛于 以 A的 特 征 值 为 对 角 元 的 对 角 阵 . 记 其 中 Ek是 Ak除 主 对角 元 外 的 矩 阵 .由 平 面 旋 转 矩 阵 的 性 质 中 ,对 于 ,有因 此 , ( )( )kk ii kA diag a E 1 Tk pq k pqA R A R ,i p q( 1) 2 ( 1) 2 ( ) 2 ( ) 2( ) ( ) ( ) ( )k k k kip iq ip iqa a a a 2 2 ( ) 21 2( )kk k pqF FE E a 又 由 假 设 ,因 此 ,这 样 ,便 有从 而 ,当 2 2( ) ( ) ( ) ( )1 , , 1max ( 1)nk k k kpq ij pq pqi j n i j i j i ja a a n n a 且 且且 22 ( )( 1) kk pqFE n n a 2 2 21 12 22 2(1 ) (1 )kk kF F FE E En n n n 1, 0k Fk E 9.3.3 实 用 的 Jacobi方 法 循 环 Jacobi方 法 必 须 一 次 又 一 次 扫 描 ， 才 能 使 Ak收 敛 于 对 角 阵 ， 计 算 量 很 大 . 在 实 际 计 算 中 ， 往往 用 一 些 特 殊 方 法 来 控 制 扫 描 次 数 ， 减 少 计 算 量 . 下 面 介 绍 一 种 应 用 最 为 广 泛 的 特 殊 循 环 Jacobi方法 阈 Jacobi方 法 . 阈 Jacobi方 法 首 先 确 定 一 个 阈值 ， 在 对 非 对 角 元 零 化 的 一 次 扫 描 中 ， 只 对 其 中绝 对 值 超 过 阈 值 的 非 对 角 元 进 行 零 化 . 当 所 有 非对 角 元 素 的 绝 对 值 都 不 超 过 阈 值 后 ， 将 阈 值 减 少 ， 再 重 复 下 一 轮 扫 描 ， 直 至 阈 值 充 分 小 为 止 . 减 少阈 值 的 方 法 通 常 是 先 固 定 一 个 正 整 数 Mn， 扫 描一 次 后 ， 让 . 而 阈 值 的 下 界 是 根 据 实 际 问 题的 精 度 要 求 选 定 的 . M 9.3.4 用 Jacobi方 法 计 算 特 征 向 量 假 定 经 过 k次 迭 代 得 到 Ak+1=RTkRT1AR1Rk， (15) 这 时 Ak+1是 满 足 精 度 要 求 的 一 个 近 似 的 对 角 阵 . 如 果记 Qk=R1R2Rk=Qk-1Rk， (16) 那 么 ， Qk是 一 个 正 交 矩 阵 ， 且 (15)式 又 可 表示 为 Ak+1=QTkAQk.当 Ak+1的 非 对 角 元 素 充 分小 ， Qk的 第 j列 qj可 以 看 成 是 近 似 特 征 值a(k+1)jj相 应 的 特 征 向 量 了 . 在 实 际 计 算 中 ， 可 以 按 (16)式 在 迭 代 过程 中 形 成 Qk， 把 Qk看 成 是 Qk-1右 乘 一 个 平 面旋 转 矩 阵 得 到 . 不 妨 记 Q0=I， Qk的 元 素 按 下式 计 算 ：( ) ( 1) ( 1) ( ) ( 1) ( 1)( ) ( 1)cos sin ,sin cos , , ,1,2,k k kip ip k ip kk k kiq ip k iq kk kij ijq q qq q qq q j p qi n 9.4 对 分 法 理 论 上 ， 一 个 实 对 称 矩 阵 正 交 相 似 于一 个 以 其 特 征 值 为 对 角 元 的 对 角 阵 . 但 是 ，经 典 的 结 果 告 诉 我 们 ， 一 个 大 于 4次 的 多项 式 方 程 不 可 能 用 有 限 次 四 则 运 算 求 根 . 因 此 ， 我 们 不 可 能 期 望 只 用 有 限 次 相 似变 换 将 一 个 实 对 称 矩 阵 约 化 为 一 个 对 角阵 .下 面 先 介 绍 将 一 个 实 对 称 矩 阵 相 似 约化 为 实 对 称 三 对 角 矩 阵 的 方 法 ， 再 讨 论求 其 特 征 值 的 对 分 法 . 9.4.1 相 似 约 化 为 实 对 称 三 对 角 矩阵 将 一 个 实 对 称 矩 阵 正 交 相 似 约 化 为 一 个 实 对 称三 对 角 矩 阵 的 算 法 ， 可 归 纳 如 下 ： 记 A(1)=A,对 k=1,2, ,n-2 按 (4)式 、 (5)式 和 (8)式 计 算 ； 按 (9) (12)式 ， 计 算 A(k+1). ,k k ku 和 ( )1,( )2,( ) (4)kk k kkk kk knkaau a ( ) ( ) 21, 1( ) ( ) (5)nk kk k k iki ksign a a ( )1,( )(8)kk k k k ka 1 ( )1 ( 1) ( ) ,(9),(10)1 ,(11)2 ( ),(12)kk k kTk k k kk k k k k k T Tk k k k g A uu gy g uA A u y y u 9.4.2 Sturm序 列 的 性 质 设 实 对 称 三 对 角 矩 阵 为 其 中 i0 (i=1,2,， n-1) 其 特 征 矩 阵 为 T-I. 记 T-I的 第 i阶 主 子 式 为 1 11 2 2 1n naT a a 1 11 2 11( )i ii ia ap a 这 是 关 于 的 i次 多 项 式 ， 当 i=n时 ， pn()=T-I 是 矩 阵 T的 特 征 多 项 式 . 令 p0()1，则 有 p1()=1-， pi()=(i-)pi-1()-2i-1pi-2()，i=2,3， ， n.(15) 多 项 式 序 列 pi() (i=0,1,， n)称 为Sturm序 列 1 11 2 11( )i ii ia ap a 定 理 9.4.1 pi() (i=1,2,， n)的 根 都 是 实 根 . 证 由 (14)式 ， pi()是 i阶 实 对 称 矩 阵 的 特 征 多 项式 ， 因 此 ， pi() (i=1,2,， n)的 根 全 是 实 根 . 定 理 9.4.2定 理 9.4.2 设 是 pi()的 一 个 根 ， 那 么 p i-1( )pi+1( ) 0， 即 相 邻 的 两 个 多 项 式 无 公 共 根 ； pi-1( )pi+1( )0， 即 pi-1( )与 pi+1( )反 号 . lim ( ) 0,ip lim ( ) 0, lim ( ) 0,i ip p 当 i为 奇 数 当 i为 偶 数定 理 9.4.4 pi()的 根 都 是 单 根 ， 并 且 将 pi+1()的 根 严 格 隔离 . 9.4.3 同 号 数 和 它 的 应 用 定 义 1 设 p0()1， pi() (i=1,2,， n) 是 一 个 Sturm序 列 ， 称 相 邻 的 两 个 数 中 符号 一 致 的 数 目 为 同 号 数 ， 记 为 ai(). 若 某个 pi()=0， 规 定 与 pi-1()反 号 . 定 理 9.4.5 设 两 个 实 数 x 3 n ， 可 以 用 乘 幂法 计 算 2及 其 相 应 的 特 征 向 量 . 在 计 算 1和 v 后 ， 按 (15)式 形 成 n-1阶 矩 阵 B的 计 算 过 程称 为 收 缩 方 法 . 9.6 反 幂 法 反 幂 法 可 以 求 一 个 非 奇 异 矩 阵 A的 逆 矩 阵 A-1的 按 模 最 小 的 特 征 值 及 相 应 的 特 征 向 量 ， 又可 以 求 A的 一 个 近 似 特 征 值 相 应 的 特 征 向 量 . 9.6.1 求 按 模 最 小 特 征 值 及 相 应 特 征 向 量 的反 幂 法 ,又 称 为 反 迭 代 法 . 1( 1) 1 , 0,1,kk kk kTk k kTk kxz xLUx zz x kz z 9.6.2 求 近 似 特 征 值 的 特 征 向 量 的 反 幂 法 先 对 矩 阵 进 行 LU分 解 ， 记 那 么 ， (7) 下 面 介 绍 一 种 选 取 特 殊 的 初 始 向 量 x0的 反 幂法 半 迭 代 法 . lI A lI A LU 1 , 0,1,kk kk kk kxz xLy zUx y k 假 设 ,选 取 初 始 向 量 x0满 足 x0=1，这 时 z0=x0.对 照 (7)式 中 的 第 二 个 式 子 .可 把 z0看 成 满 足 Le=z0.(8) 这 里 ， e=(1,1， ， 1)T， 而 z0的 各 个 分 量 的取 值 多 少 是 无 关 重 要 的 .这 样 ， 在 第 一 个 迭代 步 的 计 算 中 ， 只 需 求 解 (7)式 中 的 上 三 角方 程 组 Ux1=e. “半 迭 代 法 ” 的 命 名 也 由 此 而得 . lI A 9.7 QR方 法 定 理 9.7.1 设 A是 n阶 矩 阵 ， 其 n个 特 征 值为 .那 么 存 在 一 个 酉 矩 阵 U， 使U AU是 以 为 对 角 元 的 上 三 角 矩 阵 . 1 2, , , n 1 2, , , n 9.7.1 两 个 基 本 定 理定 理 9.7.2 设 A是 n阶 实 矩 阵 ， 那 么 ， 存 在 一个 正 交 矩 阵 Q， 使 Q AQ为 一 个 准 上 三 角 矩 阵 ，它 的 对 角 元 是 A的 一 个 特 征 值 ， 对 角 元 上 的 二阶 块 矩 阵 的 两 个 特 征 值 是 A的 一 对 共 轭 复 特 征值 . 9.7.2 相 似 约 化 为 上 Hessenberg矩 阵 对 一 般 n阶 矩 阵 ， QR算 法 的 每 一 个 迭 代 步 需 要O(n )次 乘 法 运 算 .如 果 矩 阵 阶 数 稍 大 ， 这 个 算法 几 乎 没 有 实 际 的 应 用 价 值 .通 常 采 用 的 方 法 是 先 将 矩 阵 相 似 约 化 为 上Hessenberg形 式 的 矩 阵 ， 在 此 基 础 上 应 用 QR迭 代 .这 时 ， 一 个 QR迭 代 步 的 乘 法 运 算 次 数 只需 O(n ） 次 . 所 谓 上 Hessenberg矩 阵 是 指 一 个 n阶矩 阵 A， 如 果 当 ij+1时 ， aij=0， 称 A为上 Hessenberg矩 阵 .例 如 ： 一 个 5阶 的上 Hessenberg矩 阵 具 有 如 下 的 形 式 ： * * * * * * * * *0 * * * *0 0 * * * 0 0 0 * *A 下 面 介 绍 QR方 法 时 ， 都 假 设 矩 阵 A是 一 个上 Hessenberg矩 阵 .