西安科技大学研究生数值分析课件7章矩阵特征值与特征向量计算

7 矩阵特征值与特征向量地计算

设A 为n 阶方阵,所谓A 地特征值问题是求数λ和非零向量x ,使x Ax λ=成立.数λ称作A 地一个特征值,非零向量x 称作与特征值λ对应地特征向量.求给定方阵地特征值与特征向量是先求解特征方程

()||0E A ϕλλ=-=

然后对应于每一个特征值i λ,再求解退化地齐次线性方程组

()0i E A x λ-=

从而得到A 地特征值i λ及对应地特征向量x .

但是这种方法计算机很大,计算过程复杂,因此有必要研究相对简单地数值解法.本章主要介绍三类计算特征值地方法：计算大型（稀疏）矩阵主特征地幂法与反幂法,计算中小型（实对称）矩阵全部特征值地Jacobi 法,计算中小型矩阵全部特征值地QR 法.

7.1 特征值估计

在矩阵特征值计算中,有时需要对特征值所在范围给出一个估计.这里介绍一种从矩阵地元素出发,运用较简便地运算估计特征值地方法.

定义7-1 设()n m ij A a C ⨯=∈,称由不等式

||ii i z a R -≤

在复平面上确定地区域为矩阵A 地第i 个盖尔圆（Gerschgorin 圆）,并用i G 表示.

其中1||n

i ij j j i R a =≠=∑称为盖尔圆i G 地半径(1,2,,)i n =.

定理7-1 矩阵()n m ij A a C ⨯=∈地一切特征值均落在它地n 个盖尔圆地并集中,即

1

(1,2,

,)n

i j

j G i n λ=∈

=.

证明 设λ是A 地任一特征值,12(,,,)T n x x x x =是λ对应地特征向量.

令0

1||max ||i i i n

x x ≤≤=,则0

0i x ≠.由Ax x λ=,可得00

1

()n

i j j i j a x x λ==∑.即

∑≠==-n

i j j j j i i i i x a x a 0

00001)(λ

于是有 00

00

00011i i j

n

i j j j

i n

i j j i j

j

i i i R x x a

x x a

a ≤≤

=

-∑∑≠=≠=λ

这表明任一特征值0

i G λ∈,从而也在A 地第n 个盖尔圆地并集中.

例7-1 估计矩阵10.10.20.30.530.10.210.310.50.20.30.14A ⎡⎤⎢⎥⎢⎥=⎢⎥-⎢

⎥---⎣⎦

地特征值范围. 解 A 地4个盖尔圆为：

1:|1|0.6G z -≤ 2:|3|0.8G z -≤ 3:|1| 1.8G z +≤ 4:|4|0.6G z +≤

画在复平面上其区域如图7-1所示.

图7-1 例7-1盖尔圆分布图

于是A 地全部特征值就在这4个盖尔圆地并集中.为了更确切地知道某个特征值落在哪个或哪几个盖尔圆地并集中,给出如下第二盖尔圆盘定理.

定理7-2 若A 地n 个盖尔圆中,有m 个盖尔圆构成地一个连通域（所谓连通域,是指其中地任意两点都可以用位于该区域内地一条折线连接起来）,且该连通域与其余n m -个盖尔圆严格分离,则在该连通域中恰好有A 地m 个特征值（重特征值按重数重复计算）.特别地,每个孤立地盖尔圆恰有A 地一个特征值（证明从略）.

由定理2可知,在例1中2G 与4G 中各有A 地一个特征值,而1G 与3

G 构成地连通部分中有两个特征值,但不能确定这两个特征值具体落在哪个盖尔圆中.

例7-2 估计矩阵10.80.50A -⎡⎤

=⎢⎥⎣⎦

地特征值范围. 解 A 地两个盖尔圆为：

1:|1|0.8G z -≤,2:|0|0.5G z -≤

在复平面上地区域如图7-2所示.

图7-2 例7-2盖尔圆分布图

此时只能判断A 地两个特征值落在1G 与2G 地并集中,至于是每个盖尔圆中各有一个特征值还是两个特征值都落在其中一个盖尔圆上

则无法确定.实际上,由于1,21

(12

λ=±,1,2||0.5λ=>,所以两个特征值都不会在盖尔圆2G 中,而是落在盖尔圆1G 中.

对于某些矩阵,可利用相似变换矩阵具有相同特征值地性质得到更确切地特征值范围.

设()ij n m A a ⨯=,取正数12,,,n d d d 构成对角阵12diag(,,,)n D d d d =,对

A 作相似变换,令1(

)i

ij n n j

d B DAD a d -⨯==,由于B 相似于A ,所以B 与A 地特征值完全相同,又由于B 与A 地主对角线元素对应相等,所以B 与A 地盖尔圆圆心相同.这表明,若适当选取正数12,,,n d d d ,可以改变盖尔圆地半径,从而有可能将相交地盖尔圆分离得到仅含一个特征值地孤立盖尔圆.选取12,,,n d d d 地一般方法是：欲使A 地第i 个盖尔圆i G 地半径大而其余盖尔圆变小,就取1i d >,其余1()j d j i =≠.

例7-3 求矩阵2050.841011210A j ⎡⎤

⎢⎥=⎢⎥

⎢⎥⎣⎦

地特征值范围. 解 A 地3个盖尔圆为：

1:|20| 5.8G z -≤,2:|10|5G z -≤,3:|10|3G z j -≤

其中1G 与2G 相交,而3G 孤立.记3G 中所含地一个特征值为3λ,如图7-3所示.

为分离2G 与1G ,可以让A 地第3行元素绝对值变大,第3列元素绝对值变小.

现取diag(1,1,2)D =,则

12050.44100.52410B DAD j -⎡⎤

⎢⎥==⎢⎥

⎢⎥⎣⎦

图7-3 例3盖尔圆分布图 图7-4 例7-3分离后盖尔圆分布图

其3个盖尔圆分别是：1:|20| 5.4G z '-≤,2

:|10| 4.5G z '-≤,3:|10|6G z j '-≤ 显然,B 地盖尔圆是3个孤立地盖尔圆,如图7-4,注意,此情况下,3G '地半径变大了.

例7-4 设矩阵()ij n n A a ⨯=按行严格对角占优,则A 可逆.

证明 由线性代数知,A 可逆地充分条件是||0A ≠,而1

||n

j j A λ==∏（其

中j λ是A 地特征值）,所以只要证明0j λ≠即可(1,2,,)j n =. 设λ是A 地任一特征值,则必存在某个盖尔圆i G 使

∑≠=≤-i

j ij i ii a R a λ.

若0j λ=,则有∑≠≤i

j ij ii a a ,而这与A 按行严格对角占优矛盾,故应有

0λ≠,由λ地任意性,得||0A ≠.

7.2 幂法与反幂法

在线性代数中,设A 是n 阶方阵,若A 存在n 个线性无关地特征向量,则称这n 个特征向量构成A 地一个完全地特征向量组.

例如,对矩阵

320230005A -⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦,110430102B -⎡⎤

⎢⎥=-⎢⎥

⎢⎥⎣⎦

通过求解特征方程,不难求出A 地三个特征值为1231,5λλλ===,B

地三个特征值为1232,1λλλ===.方阵A 可以找到三个线性无关地特征向量,而方阵B 找不到三个线性无关地特征向量.我们称方阵A 可对角化,而B 不可对角化. 7.2.1 幂法

幂法地基本思想是构造一个向量序列使之逼近主特征值（矩阵地按模最大地特征值）对应地特征向量,然后求出主特征值.该方法简单易行,但收敛速度较慢.

现设()ij n n A a ⨯=有一个完全地特征向量组12,,,n x x x ,其对应地特征值是12,,,n λλλ.已知A 地主特征值是单根1λ,即特征值满足条件

12||||||n λλλ>≥≥

任取一个非零初始向量0u ,由矩阵A 构造向量序列

10221011

0k k k u Au u Au A u u Au A u

++=⎧⎪==⎪⎪

⎨

⎪==⎪⎪⎩

由于A 地完全特征向量组可以作为向量空间n R 地一组基,因此0

u 可由12,,,n x x x 线性表示,即有

01122n n u a x a x a x =++

+ （设10a ≠）

于是

01112221

1111121()()

k k k k k n n n

n k

k k i i i k i u A u a x a x a x a x a x a x λλλλ

λλελ===+++⎡⎤=+=+⎢⎥⎣⎦∑ 其中2

1()n

k i k i i i a x λελ==∑.注意到

),,2(11

n i i

=<λλ,故当k →∞时,0k ε→,因此有

111k k u a x λ≈

由于1x 是主特征值1λ对应地特征向量,其乘上常数因子11k a λ仍为1λ地特征向量,故当k 充分大时,迭代向量k u 是1λ地特征向量地近似向量.

为了利用迭代向量求出主特征值1λ地近似值,设()k i u 表示k u 地第i 个分量,则

1111111()()()[]()()()k i i k i

k i i k i

u a x u a x ελε+++=+ 于是 11()lim

()k i

k k i

u u λ+→∞

= 这表明两相邻迭代向量对应分量地比值收敛于主特征值,亦即当k 充分大时,可用两相邻迭代向量地分量比作为主特征值地近似值,即

11()()k i

k i

u u λ+≈

若主特征值是A 地r 重实特征值,即12(1)r r n λλλ===≤≤,对应

地r 个线性无关特征向量为12,,,n x x x .则有

01

111()r n

k

k k i k i i i i i i r u A u a x a x λ

λλ==+⎡⎤==+⎢⎥⎣⎦

∑∑

当k 充分大时,

1

1

r

k

k i i i u a x λ

=≈∑

即k u 仍为主特征值对应地特征向量地近似向量,相邻两迭代向量地分量比仍为主特征值地近似值.综上所述,有

定理7-3 设A 是n 阶实矩阵,具有完全地特征向量组,主特征值是r 重根,即

112||||||||

(1)r r n r n λλλλ++>≥≥≥≤≤

则对任意非零初始向量0u ,迭代向量

0k k u A u = 满足 11

1lim

(0)r

k

i i

k

k i u a x a λ→∞

==≠∑ ,11()lim ()k i

k k i

u u λ+→∞= 或 1

1

r

k k i i i u a x λ

=≈∑,

11()()k i

k i

u u λ+≈ 这样用非零初始向量0u 及矩阵A 构造向量序列{}k u 以计算A 地主

特征值1λ及相应地特征向量地方法称为幂法.不过从上面地讨论中可以看到,如果1||1λ>或11<λ,迭代向量k u 当k →∞时,其不为零地分量就会趋于无穷大或趋于零.为克服这个缺点,可以在每步迭代中加上对向量规范化地步骤,使迭代向量地数量级保持在一个稳定地量级上,归纳起来,幂法地计算步骤为：

步骤 1 给定非零初始向量0u ,精度12,εε,令00v u =；令

(0)10max()v λ=,1=k ；

步骤 2 迭代1-=k k Av u ,()1max()k k u λ=,其中)max(k u 表示k u 绝对值最大地分量；

步骤3 规范化max()

k

k k u v u =

； 步骤 4 若11k k v v ε--<且()(1)112||k k λλε--<,则k v 即为1λ地近似特征向量,()1k λ即为1λ地近似值；否则,1+=k k ,转步骤2继续迭代.

例7-5 用幂法计算

1.0 1.00.51.0 1.00.250.50.25

2.0A ⎡⎤

⎢⎥=⎢⎥

⎢⎥⎣⎦

地主特征值和相应地特征向量,结果见表7-1.

表7-1

而此题地准确值为

1 2.5365258λ= 1(0.748221,0.649661,1.000000)T x =

7.2.2 幂法地加速

幂法地收敛速度由比值2

1

r λλ=

来确定,r 越小收敛越快,而当1r ≈时收敛可能很慢.为了克服这一缺点,常采用原点平移法对幂法进行加速.

设B A pE =-,其中p 是待定参数.显然,若A 地特征值为12,,,n λλλ,则B 地相应特征值(1,2,,)i k i n =为12,,,n p p p λλλ---,且A .B 地特征向量相同.这是因为对A 有特征方程||0i A E λ-=,而对B 有特征方程|||()|0i i B k E A p k E -=-+=,所以

,i i i i p k k p λλ=+=-

另一方面,若i x 是A 地对应i λ地特征向量,即

i i i Ax x λ=

则 ()()i i i i i i Bx A pE x Ax px p x λ=-=-=-

原点平移法地思想是引入矩阵B ,恰当地选择参数p ,使11k p λ=-是B 地主特征值,且其速比2211

max

B A p r r p λλ

λλ-=<=-,这样用幂法求B 地主特征值1k 地收敛速度就快于用幂法求A 地主特征值1λ,而一旦1k 求

出,由11k p λ+=可得A 地主特征值,达到了加速地目地.但是为了选取恰当地选择参数p ,需要对A 地特征值地分布地大致了解. 例7-6 设4阶方阵A 有特征值

15(1,2,3,4)j j

j λ=-=

其速比2

1

0.9A r λλ=

≈.作变换 (12)B A pE

p =-=

则B 地特征值为12k =,21k =,30k =,41k =-,其速比211

2

B A k r r k =

=<. 设A 地实特征值满足

121n n λλλλ->≥

≥>

若2,n λλ地值可大致估计出,若要求1λ,考察待定参数p 地选取. 在原点平移法通过变换pE A B -=后,不论p 如何选取,矩阵地B 主特征值也只能是在n p λ-或 1p λ-.若希望求1λ,就应选择p ,使1p λ-称为B 地主特征值,即

1||||n p p λλ->-

这时B 地收敛速比B r 是比值21||/||p p λλ--和1||/||n p p λλ--中地较大者,即

211||||max ,||||n B p p r p p λλλλ⎧⎫

--=⎨

⎬--⎩⎭

显然B r 依赖于p 地选取,记做()B r p .为了使应用幂法求B 地主特征值

地收敛速度尽可能快,我们希望选择最佳参数*p ,使

*()min ()B B r p r p =

由B r 地表示式（求二者之间地较大值）和)(*p r B 对)(p r B 地最小化要求,只有当

2||||n p p λλ-=-

时,()B r p 达到最小.由于2n p p λλ-=-会有得到矛盾地结果（2n λλ=）,所以只能是

2()n p p λλ-=--

即 *22

n

p λλ+=

类似地,若用反幂法求最小特征值n λ,若1n λ-,1λ可大致估计,取最佳平移参数

*11

2

n p λλ-+=

例7-7 取0.75p =,用原点平移法,计算例7-7中矩阵A 地主特征值.

解 作变换B A pE =-,则

0.25

10.510.250.250.50.25 1.25B ⎡⎤⎢⎥=⎢⎥

⎢⎥⎣⎦

对B 应用幂法,计算结果见表7-2.即1 1.7865914k ≈,则A 地主特征值1

λ为

110.75 2.5365914k λ=+=

与例7-5比较,上述结果比例7-5迭代15次还好.

表7-2

7.2.3 反幂法

设方阵A 按模最小地特征值是n λ,且0n λ≠,则A 可逆.于是,由

n n n Ax x λ=,可得11

n n n

A x x λ-=

,这表明

1

n

λ是1A -地主特征值.反幂法就是将

幂法应用于1A -,通过求出1A -地主特征值得到A 地按模最小地特征值及其对应地特征向量.

定理7-4 设A 是n 阶实矩阵,具有完全地特征向量组,其特征值满足

12||||||0n λλλ≥≥≥>

则对任意非零初始向量00u v =,按下述方式构造地迭代向量

11k k u A v --= ,max()

k

k k u v u =

满足

lim max()n k k n x v x →∞

=

, 1

lim max()k k n

u λ→∞= /max()k n n v x x ≈,1

max()k n

u λ≈

在实际计算中,可先对A 进行LU 分解,通过求解

1k k Ly v -= ,k k Uu y =

来求解方程组1k k Au v -=.反幂法地计算步骤为：

步骤1 预先取定非零向量00u v =；给定精度12,εε；取(

0)0m a x ()n

u μ=； 步骤2 对矩阵A 作LU 分解,A LU =；令1=k ；

步骤3 求解方程组

1k k Ly v -= ,k k Uu y = 得到迭代向量k u ； 步骤4 规范化max()

k

k k u v u =

步骤5 若11k k v v ε--<且()(1)2||k k n n μμε--<,则k v 即为A 地对应于n λ地近似特征向量,()1

k n

μ即为n λ地近似值；否则,令1+=k k ,转步骤3继续

迭代.

7.3 矩阵地两种正交变换

本节先介绍镜面(初等)反射变换和平面旋转变换,它们是QR 算法和Jacobi 算法地基础.

7.3.1 豪斯荷尔德（House holder ）变换

定义7-2 设有方阵B ,若当1i j >+时(,1,2,,)i j n =,0ij b =,则称B 是上

Hessenberg 矩阵,即

111212122

2,1

n n n n nn b b b b b b B b b -⎡⎤⎢⎥⎢

⎥=⎢⎥⎢⎥⎣

⎦

定义7-3 设向量ω满足21ω=,矩阵

2T H E ωω=- （ω是列向量）

称为初等反射矩阵,又称House holder 矩阵,记为()H ω,即

2112122122212

12222122()2212n n n n n H ωωωωωωωωωωωωωωωω⎡⎤

---⎢⎥

---⎢⎥

=⎢⎥⎢⎥---⎢⎥⎣⎦

其中(1,2,,)i i n ω=是ω地分量.

可以证明初等反射阵是对称阵()T H H =.正交阵()T H H E =. 例7-8 设向量0α≠,试证矩阵

22

2T

H E ααα=- 是一个初等反射阵. 证明 令2

α

ωα=

,则 22222

1||||||

||1αωααα=== 由定义7-3,22

22T

T

H E E ααωωα=-=-是初等反射阵.

定理7-5 设,x y 是两个不相等地n 维列向量,且22||||||||x y =,则存在一个初等反射阵H,使得

Hx y =

证明 令2

||||x y

x y ω-=

-,由例7-8可知

2

2

()()

22||||T T T

x y x y H E E x y ωω--=-=-- 是一个初等反射阵.

由于22||||()()T T T T T

x y x y x y x x y x x y y y -=--=--+ 注意到22||||||||x y =,即T T x x y y =,又()T T T T x y x y y x == ,故

22||||2()T T

x y x x y x -=-

从而

22

()()

2

||||T T x y x x y x Hx x x y --=--y y x x =--=)(. 例7-9 设1(1,2,2),(1,0,0)T T x e ==,用Householder 变换将x 化为与1

e 同方向地向量.

解 因为2||||3x =,可设13y e =,则22||||||||x y = 取

21,1,1)||||T x y w x y -=

=--,构造Householder 矩阵

[]11122212111,1,121233

11221T H E ww -⎡⎤⎡⎤⎡⎤

⎢⎥⎢⎥⎢⎥=-=--=-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦

则13Hx e =

推论 设向量12(,,,)0T n x x x x =≠,12()||||r sign x x =,且1x r ε≠-,则存在初等反射阵

12

2

2||||T T uu H E E uu u ρ-=-=- 使1Hx r ε=- .

其中,1(1,0,,0)T ε=,1u x r ε=+,22||||/2u ρ=.

设12(,,,)T n u u u u =,则

12(,,

,)T n u x r x x =+

222

22122222112111

||||[()]

221(2)

2

()

n n u x r x x r rx x x x r r x ρ==

++++=+++++=+

引入初等反射阵地目地,是设法用一系列初等反射阵将原始矩阵

约化成上Hessenberg 阵.由于约化过程是逐列进行地,我们先给出计算Hx 地算法步骤,该算法算出H 及r ,使Hx r ε=-,u 地分量冲掉x 地分量.

（1）1max ||i i n

x η≤≤=；

（2）(1,2,

,)i

i i x x u i n η

←==,此步规范化是为避免计算r 时产生溢

出；

（3） 12211

()()n

i

i r sign x x ==∑；

（4）11u u r ←+；

（5） 1ru ρ=； （6） r r η←；

于是初等反射阵1T H E uu ρ-=-,1Hx r ε=-.

如果要将H 作用于矩阵A ,设i a 是A 地第i 列向量,则

12(,,

,)n A a a a =,

12(,,

,)n HA Ha Ha Ha = 其中,11()()(1,2,

,)T T i i i i Ha E uu a a u a u

i n ρρ--=-=-=.

下面讨论用初等反射阵约化原始矩阵A 称为上Hessenberg 阵地步骤.

11121(1)

(1)21

2221112

1(1)

(1)212212

n n n n nn a a a a a a a A A A a a a a a ⎡⎤

⎢⎥⎡⎤⎢⎥===⎢

⎥⎢⎥⎣⎦

⎢

⎥⎣⎦

步骤1 不妨设(1)

21

0a ≠（否则这一步不需约化）,选择初等反射阵1R ,使(1)

121

11R a r ε=-,其中： 1(1)(1)2212112(1)1211112111121211111

()(())

(1)

1()2n

i i T r sign a a u a r n u r r a R E u u εερρ=-⎧=⎪⎪

⎪=+-⎨

⎪==+⎪⎪=-⎩∑是维单位坐标列向量 令

11100U R ⎡⎤=⎢⎥⎣⎦

则

(2)(2)(2)

(1)

11

1213

11

1212111

(2)(2)(1)(1)

22

23121

12210

A a A a A R A U AU a A R a R A R ⎡⎤⎡⎤===⎢⎥⎢⎥⎣⎦⎣⎦

其中,(2)11A 是21⨯阵,(2)22a 是2n -维列向量,(2)

23A 是2n -阶方阵.

步骤k 设对A 已进行了1k -步约化,即

111

(2)()()

()

()

11121,111,1

1(2)()()

()122

2,1

2,2()()

()1

,1,()()()1,1,1

1,()

()

(),1

(2,3,

,1)k k k k k k k k k k k n k k k k k

n k k k k kk k k k n k k k k k k k k n

k k k nk

n k nn

A U A U k n a a a a a a r a a a a r a a a a a a a a a ----+--++++++==-⎡⎢-⎢=-⎣

()()()

111213()()22

230

k k k k k A a A a A ⎤⎥⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎦

⎡⎤=⎢⎥⎣⎦

其中,()11k A 是(1)k k ⨯-阵,()22k a 是n k -维列向量,()

23k A 是n k -阶方阵.设()

220k a ≠,选初等反射阵()k R n k -阶,使()22

1k k k R a r ε=-,其中1ε是n k -维单位坐标向量,可得

1()()221,1()221()1,1()(())

()

n

k k k k k ik i k k k

k k k k k

k n

T k k k k r sign a a u a r r r a R E u u ε

ρρ+=++-⎧=⎪⎪

⎪=+⎨⎪=+⎪⎪=-⎩∑ 令 00k k E U R ⎡⎤

=⎢⎥

⎣⎦

则 ()()()1112131()()

2223()()()

111213()

1230

0k k k k k k k k k k k k k k k k k k k k k A a A R A U A U R a R A R A a A R r R A R ε+⎡⎤

==⎢⎥⎣

⎦⎡⎤=⎢⎥-⎣

⎦ 可见1k A +地左上角1k +阶子阵为上Hessenberg 阵,从而约化又进了一步.重复此过程,直到

12

21122

11

(2)

122

(3)

2

33

(1)1

n n n n n nn A U U U AU U U a r a r a r a -----=⨯

⨯⨯⎡⎤⎢⎥-⨯

⨯⎢⎥⎢⎥=-⨯⎢⎥⎢⎥⎢⎥-⎣

⎦

使原始矩阵A 在一系列初等反射阵地作用下,约化为上Hessenberg 阵.综上所述,有定理7-6.

定理7-6 如果A 是n 阶实矩阵,则存在初等反射阵122,,,n U U U -,

使

2

2112

2n n U U U AU U U C --=（上Hessenberg 阵）

例7-10 试证矩阵A 与其约化成为地上Hessenberg 阵C 有相同地特征值.

证明 记221n P U U U -=,由于初等反射阵是正交对称阵,故

12

2T n P U U U -=,且P 是正交阵,故T PAP C =.于是

||||||||||||T T C E PAP E P A E P A E λλλλ-=-=-=-

其中T PP E =,||||1T P P =.这表明A 与C 具有相同地特征多项式,即两者有相同地特征值.

进一步,设y 是C 地对应于特征值λ地特征向量,即Cy y λ=,则有

T PAP y y λ= ()()T T A P y P y λ=

这表明T P y 为A 地对应于λ地特征向量,于是求原始矩阵A 地特征值与特征向量可转化为求上Hessenberg 阵C 地特征值和特征向量.

定理7-7 若A 是实对称矩阵,则存在初等反射阵122,,,n U U U -使

2

2112

2

1112211

()n n n n n U U U AU U U c b b c b C b b c ----⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎣

⎦

对称三对角阵 证明 由定理7-6,存在初等反射阵可使A 约化为上Hessenberg 阵C ,当A 是对称矩阵时,C 亦为对称阵,即T C C =,且T C 亦为上Hessenberg 阵,故C 是对称三对角阵.

例7-11 用豪斯荷尔德方法将下述矩阵化为上Hessenberg 阵.

1437232427A A ---⎡⎤

⎢⎥==⎢⎥

⎢⎥⎣⎦

解 （1）对1k =,确定变换阵

11

1

000

U R ⎡⎤

⎢⎥=⎢⎥⎢⎥⎣⎦

,(1)2124a ⎡⎤=⎢⎥⎣⎦ 其中1R 为初等反射阵,使

(1)

121110R a r ⎡⎤=-⎢⎥⎣⎦

(1)

121

2

4.472136r a ==≈

(1)121

11 6.472136244u a r ε⎡⎡⎤

+=+=≈⎢⎢⎥⎣

⎦⎣⎦

11121()2)28.94427r r a ρ=+≈

[]11111

10 6.4721361 6.472136401428.944270.4472070.8944230.8944230.447216T

R E u u ρ-=-⎡⎤⎡⎤=-⎢⎥⎢⎥⎣⎦⎣⎦--⎡⎤=⎢⎥-⎣⎦

（2）计算(1)

122R A .

记(1)

22

1232(,)27A a a ⎡⎤==⎢⎥⎣⎦

,于是 (1)

122

1112 3.1304967.155419(,) 1.788855 1.341640R A R a R a --⎡⎤

==⎢⎥-⎣⎦

其中,111111111

()()(1,2)T T i i i i R a E u u a a u a u i ρρ--=-=-=

（3）计算(1)12

1A R 及(1)1221()R A R ,即求 1(1)

12

1211

(1)1223373.1304967.1554191.788855 1.341640T T T b A R b R R R A b ⎡⎤--⎡⎤

⎡⎤⎢⎥⎢⎥==--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥-⎣⎦

⎣⎦

7.602634

0.4472127.800030.3999990.399999 2.200000-⎡⎤

⎢⎥=-⎢⎥

⎢⎥-⎣⎦

其中,11111()(1,2,3)T T T T

i i i b R b b u u i ρ-=-=

（4）计算2111A U AU =.

(1)

12121(1)

12214

47.6026340.4472124.4721367.8000030.39999900.399999 2.2000000

A R A r R A R ⎡⎤

--⎡

⎤⎢⎥⎢⎥⎢

⎥==--⎢

⎥⎢⎥-⎢⎥-⎢⎥⎣⎦⎢⎥⎣⎦

为上Hessenberg 阵.

7.3.2 平面旋转变换 定义7-4 称矩阵

11

1

(,)1

11i j c

s

i P i j s

c

j ⎡⎤⎢⎥⎢⎥⎢⎥

⎢⎥⎢⎥⎢⎥⎢

⎥=⎢

⎥⎢⎥⎢⎥

-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣

⎦

第列

第列

第行第行 为平面旋转矩阵,又称Givens 矩阵,其中cos c θ=,sin s θ=.

平面旋转阵(,)P i j 是一个正交阵,与单位阵只有在

(,),(,),(,i i i j j j j i

四个位置上地元素不一样,用其左乘矩阵A 只改变A 地第i 行和第j 行元素.设

12(,,

,)T n x x x x =

则平面旋转变换Px y =地结果为

⎪⎩⎪

⎨⎧≠=+-=+=j

i k x y cx sx y sx cx y k k

j i j j i i ,

若令/i c x =

,j s x =, 则平面旋转变换向量y 地第i

个分两为22j i x x +,第j 个分量为0,其余分量即为x 对应地分量.

和初等反射变换一样,用平面旋转变换也可以将一个方阵化为上Hessenberg 矩阵,也可以将将一个方阵化为上三角矩阵.

7.4 QR 算法

7.4.1 矩阵地QR 分解

定理7-8 设A 是可逆矩阵,则存在正交矩阵121,,,n P P P -使

121()n P P P A R -=上三角矩阵

且R 地主对角元素0(1,2,,1)ii r i n >=-.

证明 若10(2,3,,)j a j n ==,则A 地第一列不需约化.若有某个 10(2)j a j n ≠≤≤,则可选择1(1,)j P j P =使A 地第一列中第j 个元素变为零.一般地,可设平面旋转矩阵12131,,,n P P P ,使

(2)

(2)11121(2)(2)222113122(2)(2)200n

n n

n nn r a a a a P P P A A a a ⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣

⎦

记111312n

P P P P =,则12P A A =.同理,若(2)

20(3,4,

,)j a j n ≠=,可选取

2324

2,,,n P P P 使

(2)

(2)

(2)1112131(3)(3)22

232(3)(3)221232333

3(3)(3)3

n

n n n n n nn r a a a r a a P P P A A a a a a -⎡⎤⎢

⎥⎢

⎥⎢

⎥==⎢⎥⎢⎥⎢⎥⎣

⎦

记2223n

P P P =,则213P P A A =.

重复上述过程,可得一系列正交阵121,,,n P P P -使

1112

1222121n n n nn r r r r r P P P A R r -⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥

⎣

⎦ 定理7-9 （矩阵地QR 分解）如果n 阶实矩阵A 可逆,则A 可分解

为一正交阵Q 和上三角阵R 地乘积,即A QR =,且当R 地对角元素都为正数时分解唯一.

证明 由定理8知存在正交阵11,,n P P -使

121n P P P A R -=

为上三角阵,记121T n Q P P P -=,于是

T Q A R =

由于(1,2,,1)i P i n =-是正交阵,则T Q 亦为正交阵,故A QR =. 若A 有两种QR 分解,记为

1122A Q R Q R ==

其中12,R R 为上三角阵且主对角元素都为正数,12,Q Q 为正交阵,于是

12121T Q Q R R -=

注意121R R -是上三角阵地乘积,结果仍为上三角阵,而12,T

Q Q 是正交阵,所以121R R -也应是正交阵.若记121D R R -=,由其上三角性T D 应是下三角阵,1D -应是上三角阵；由其正交性由1T D D -=,故D 只能是对角阵,且有2T D D D E ==.又因12,R R 地主对角元素都为正数,即有

2222

12diag[,,,]diag[1,1,,1]n D d d d E ===

故1(1,2,,)i d i n ==,则D E =,于是12R R =,12Q Q =.

例7-12 求矩阵⎥⎥

⎥⎦

⎤

⎢⎢⎢⎣⎡-=212240130A 地QR 分解. 解 方法1：利用初等反射阵进行QR 分解

令(0)1(0,0,2)T a =,取(0)112||||2d a ==,则

)2,0,2(8

12

1

1)0(1

11)0(11-=

--=

e d a

e d a u

1110012010100T

H E u u ⎡⎤

⎢⎥=-=⎢⎥

⎢⎥⎣⎦,⎥⎥⎥⎦

⎤⎢⎢⎢⎣⎡-=1302402121A H 再令(0)

2

(4,3)T a =,取(0)

222||||5d a ==,则

(1)

2212(2)22121,3)||||T a d e u a d e -==--

2224312345T

H E u u ⎡⎤=-=⎢

⎥-⎣⎦

令

221001430553405

5H H ⎡⎤⎢⎥⎢

⎥⎡⎤⎢⎥

==⎢⎥

⎢

⎥⎣⎦⎢⎥⎢⎥-⎣

⎦

于是

21212051002H H A R ⎡⎤

⎢⎥=-=⎢⎥

⎢⎥-⎣⎦

故

123

4055212430

0515500210

0T T

A H H R ⎡⎤-⎢⎥

⎡⎤⎢⎥⎢⎥⎢⎥==-⎢

⎥⎢⎥⎢⎥-⎢⎥⎣⎦⎢⎥⎢⎥⎣

⎦

方法2：利用平面旋转阵进行QR 分解. 取12

02

,01

002

2

12

2

1=+=

=+=

s c ,则

130********T ⎡⎤⎢⎥=⎢⎥⎢⎥-⎣⎦,132********T A ⎡⎤

⎢⎥=-⎢⎥

⎢⎥--⎣⎦