普通矩阵特征值的QR算法

普通矩阵特征值的QR 算法
摘 要
求矩阵的特征值有多种不同的办法，本文主要介绍用QR 法求矩阵的特征
值，QR 法是目前求中等大小矩阵全部特征值的最有效方法之一，使用于求实矩阵或复矩阵的特征值，它和雅可比法类似，也是一种变换迭代法。

关键词：QR 分解 迭代序列 特征值 Matlab
一 、QR 方法的理论：
对任意一个非奇异矩阵（可逆矩阵）A ，可以把它分解成一个正交阵Q 和一个上三角阵的乘积，称为对矩阵A 的QR 分解，即A=QR 。

如果规定R 的对角元取正实数，这种分解是唯一的。

若A 是奇异的，则A 有零特征值。

任取一个不等于A 的特征值的实数μ，则A-μI 是非奇异的。

只要求出A-μI 的特征值和特征向量就容易求出矩阵A 的特征值和特征向量，所以假设A 是非奇异的，不是一般性。

设A=A 1 ，对A 1 作QR 分解，得A 1 = Q 1R 1，，交换该乘积的次序，得A 2 = R 1Q 1=，由于Q 1正交矩阵，A 1到A 2的变换为正交相似变换，于是A 1和A 2就有相同的特征值。

一般的令A 1=A ，对k=1，2，3，…..
⎩⎨
⎧==+)
()(1迭代定义分解k
k k k k k Q R A QR R Q A
这样，可得到一个迭代序列{A k }，这就是QR 方法的基本过程。

二、QR 方法的实际计算步骤
Householder A Hessenberg B -----→=用阵作正交相似变换
上第阵一步............*
::::***⎛⎫
⎪ ⎪ ⎪* ⎪** ⎪ ⎪ ⎪ ⎪*⎝
⎭
Householder 变换：如果 v 给出为单位向量而 I 是单位矩阵，则描述上述线性变换的是 豪斯霍尔德矩阵 （v * 表示向量 v 的共轭转置）H=I -2VV*
1k k k
Given k k k B Q R B B R Q +=⎧-------→→⎨=⎩用变换产生迭代序列第二步12
*
**n λλλ⎡⎤⎢⎥⎢⎥⎢⎥⎢
⎥⎣
⎦
Householder A B -----→=用阵
作正交相似变换
（对称阵）三对角阵***
**
⎛⎫ ⎪ ⎪
⎪*
⎪*
*
⎪ ⎪* ⎪ ⎪*⎝
⎭
三、化一般矩阵为Hessenberg 阵
称形如
1112
11121
22212323331n n n n n nn nn h h h h h h h h h h h H h h ---⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦
的矩阵为上海森堡（Hessenberg ）阵。

如果此对角线元 (i=2,3,….n)全不为零，则该矩阵为不可约的上Hessenberg 矩阵 。

用Householder 变换将一般矩阵A 相似变换为 Hessenberg 矩阵。

首先，选取Householder 矩阵，使得经相似变换后的矩阵的第一列中有尽可能多的零元素。

为此，应取为如下形式
111000
0H H ⎡⎤
⎢⎥⎢
⎥=⎢⎥⎢⎥⎢⎥⎣
⎦
其中1H 为n-1阶Householder 矩阵。

于是有 111
21111
1221T
a a H H AH H a H A H ⎡⎤
=⎢
⎥⎢⎥⎣⎦
121311212131(,,
,),(,,
,),T T n n a a a a a a a a == 222222n n nn a a A a a ⎡⎤
⎢
⎥=⎢
⎥⎢⎥⎣⎦
只要取1H 使得111(,0,,0)T H a σ=，就会使得变换后的矩阵11H AH 的第一列出
现n-2个零元。

同理，可构造如下列形式Householder 矩阵
22
1000010000
00H H ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦
使得2112***
*********
**H H AH H ⎡⎤
⎢⎥⎢
⎥⎢⎥=⎢
⎥⎢⎥⎢⎥⎣⎦
如此进行n-2次，可以构造n-2个Householder 矩阵12,,,H H 2n H -，使得
221122 n n H H H AH H H H --=。

其中H 为上Hessenberg 矩阵。

特别地，当A
为实对称矩阵，则经过上述正交变换后，H 变为三对角阵。

用Householder 方法对矩阵A 作正交相似变换，使A 相似于上Hessenberg 阵，算法如下：
1,2,...,2k n =-
（1）
计算
(1)(1)11
122
1
11
1111(1)112,1
()()(()),
()
(0,...,0,,,...,)
k k T
k k n
k k k ik i k k k k k k k k k k k k nk H I U U sign a a a U a a a ρσ
ρσσσ++++++=+++++++++=-
==+=+∑，，，
（2）
计算1k H A A +→
(1)1
1
(1)
,1,,1
1()21,,n
k j l lj l k k k ij ij j i j k k n
t u a i k n
a a t u ρ+=+++=+==+=-∑（）（）
（3） 计算
1k AH A +→
(1)1
1(1)
1,...,1
(1)(2)1,...,n
k i il
l
l k k k ij ij i j i n
t a u
j k n
a a t u ρ+=+++==
=+=-∑
四、上Hessenberg 阵的QR 分解
对上Hess enberg 阵只需要将其次对角线上的元素约化为零，用Given 变换比
用Householder 变换更节省计算量。

1、 平面旋转阵（Givens 变换阵） 定义 n 阶方阵
θθθθ⎡⎤⎢⎥⎢⎥⎢⎥⎢
⎥
--⎢
⎥⎢⎥⎢⎥
=⎢
⎥
⎢⎥⎢
⎥
---⎢
⎥⎢⎥⎢⎥⎢⎥⎢⎥
>⎣
⎦,11cos sin 11sin cos 11(
)
i j i R j j i 称为平面旋转阵，或称为Givens 变换阵。

平面旋转阵R i,j 的性质：
（1） -==1,,,,,T T
i j i j i j i j R R I R R 平面旋转阵是非对称的正交阵。

（2） ,T i j R 也是平面旋转阵。

（3） ,i j R det()=1 平面旋转阵R i,j 的作用：
（1） 将向量(),,...,x x x x T
12n =的第j 个分量约化为零。

,i j R 左乘向量(),,...,x x x x T
12n =只改变X 的第i 个分量和第j 个分量。

若令=,i j y R x ，有
θθθθ=+⎧⎪
=-+⎨⎪==≠⎩cos sin sin cos 1,...,;,i i j j i j k
k y x x y x x y x k n k i j
, 调整θ，可将j y 约化为零。

令=0j y ，得θ=tan j i
x x
所以，取cos i x
C r θ==
=
，sin j x x S r
θ===
于是0i i j j y Cx Sx r y =+==
=
(),,...,,,,...,,0,,...,i j R x x x r x x x x T
1i-1i+1j-1j+1n =
（2） 将向量(),,...,x x x x T
12n =的第i+1个分量到第n 个分量约化为零。

(
)+=,1,...,,,0,,...,,i i R x x x r x x r T
1i-1i+2n =(
)++=,2,1,...,,,0,0,,...,,
i i i i R R x x x r x x r T
1i-1i+3n =
)++=+
+,,2,12,...,,,0, 0
i n
i i i i n
R R R x x x r r x
T
1i-1=
（3） 用,i j R 对矩阵A 作变换得到的结论
,i j R 左乘A 只改变A 的第i ，j 行。

,i j R T 右乘A 只改变A 的第i ，j 列。

,,i j i j R AR T 只改变A 的第i ，j 行和第i ，j 列。

2、用Givens 变换对上Hessenberg 阵作QR 分解
对上Hessenberg 阵(1)(1)
(1)11121(1)(1)(1)2122
2(1)(1)1
n
n nn nn b b b b b b B b b -⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣
⎦
通常用n-1个Givens 变换阵可将它化成上三角矩阵，从而得到B 的QR 分解式。

具体步骤为：
设(1)
210b ≠（否则进行下一步）
取旋转矩阵1
1
11
cos sin 00sin cos 00(1,2)0
01
1R ϕϕϕ
ϕ⎡⎤⎢⎥-⎢
⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣
⎦
则(2)(2)(2)1
12131(2)(2)(2)22232(2)(2)(2)1232
33
3(2)(2)1
(1,2) n
n n
nn nn r b b b b b b R B B b b b b b -⎡⎤⎢⎥⎢
⎥⎢
⎥==⎢⎥⎢⎥⎢⎥⎣
⎦
其中，(1)(1)
112111111cos , sin , b b r r r ϕϕ===设2320
（）
b ≠（否则进行下一步），再取旋转矩阵 222
2
1
0cos sin sin cos (3,2)1
1－R ϕϕϕϕ⎡⎤⎢⎥⎢
⎥⎢⎥=⎢
⎥⎢⎥⎢⎥⎢⎥⎣
⎦
则(3)
(3)(3)(3)11213111(3)(3)(3)2
23212(3)(3)(3)3331323(3)(3)(3)4341
4(3)(3)1
(3,2)n n
n n n n
n n nn nn r b b b b r b b b b b b R B B b b b b b -----⎡⎤⎢⎥⎢⎥⎢
⎥==⎢
⎥⎢⎥⎢⎥⎢⎥⎢⎥
⎣
⎦
其中(2)(2)
322222222
cos , sin , b b r r r ϕϕ===。

假设上述过程已进行了k-1步，有
1(1,)k k B R k k B -=+()
()
()()
111111
1()
()
()1
1111()()()1()()()111
1()()1
k k k k k k n n k k k k k k k n k n k k k kk kn kn k k k k k
k n k n
k k nn nn r b b b b r b b b b h h b b b b b --------++-+-⎡⎤⎢⎥⎢⎥⎢⎥⎢
⎥=⎢
⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣
⎦
设()10
k k k b +≠，取 1
1
(1,)cos sin sin cos 1
k k k
k
R k k ϕϕϕϕ⎡⎤⎢⎥⎢⎥⎢⎥⎢
⎥+=⎢
⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎣
⎦
其中()()1cos , sin k k kk k k k k k k
b b r r ϕϕ+==
， k r =
于是(1)(1)(1)
1
111
1(1)(1)1(1)(1)
(1)111111(1)(1)(1)21
21
2(1)(1)1
(1,)k k k k k n
k k k
kk kn
k k k k k k k k n k n k k k k k k n k n
k k nn nn r b b b r b b R k k B B b h b b h h b b ++++++++++++++-+++++++-+++-⎡⎤⎢⎥⎢⎥⎢⎥⎢
⎥+==⎢
⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣
⎦
因此，最多做n-1次旋转变换，即得
() (,1)(2,1)
(1,2)n H R n n R n n R B =---()()()1
12
131()()2
23
2()3
3n n n n n n n n n
n r b b b r b b R r b r ⎡⎤⎢⎥⎢
⎥⎢
⎥==⎢⎥⎢⎥⎢⎥⎣
⎦ 因为(,1),(2,3,
,)R i i i n -=均为正交矩阵，故2132
1T T
T
nn H R R R R QR -==其
中2132
1T T
T
nn Q R R R -=仍为正交矩阵。

可算出完成这一过程的运算量约为4,
比一般矩阵的QR 分解的运算量3()O n 少一个数量级。

可证明仍是上Hessenberg 阵，于是可按上述步骤一直迭代下去，这样的得到的QR 方法的运算量比基本QR 方法大为减少。

具体实例
现在我们就用上述所说的QR 方法求解矩阵532644445A -⎡⎤
⎢⎥=-⎢⎥⎢⎥-⎣⎦的全部特征值。

第一步：
将A 化成上Hessenberg 阵，取
(6,4)(6T T T u =+= 10 201T T uu I u u ⎡⎤-=-
⎢⎥⎣⎦
0.9160250.2773500.8320500.55470020.2773500.08397470.5547000.832050--⎡⎤⎡⎤
=⎢⎥⎢⎥
-⎣⎦⎣⎦
于是110000.8320500.55470000.5547000.832050 H ⎡⎤
⎢⎥=--⎢⎥⎢⎥-⎣⎦
115 1.386750 3.3282007.211102 1.2307688.15384000.153846 2.230767H H AH ⎡
⎤⎢⎥==---⎢⎥
⎢⎥-⎣⎦ H 即为与A 相似的上Hessenberg 阵。

将H 进行QR 分解，记
11, 8.774964B H r ===
111cos 50.56980. sin 0.821781r ϕϕ===- 取0.5698030.8217810(2,1)0.8217810.5698030001R -⎡⎤
⎢⎥=⎢⎥⎢⎥⎣⎦
于是18.774964 1.8015968.597089(2,1)00.438310 1.91103000.153846 2.230767R B ⎡⎤
⎢⎥=-⎢⎥
⎢⎥-⎣⎦ 再取
于是
1 (3,2)(2,1)R R B
18.774964 1.8015968.59708900.464526 2.54198200 1.471953R ⎡⎤
⎢⎥=-=⎢⎥
⎢⎥⎣⎦
10.5698030.7754030.272165(2,1)(3,2)0.8217810.5376430.18871200.3311890.943564T T Q R R ⎡⎤
⎢⎥==-⎢⎥
⎢⎥-⎣⎦ 第一次迭代得211 3.519482 4.92549110.8401170.381739 1.091627 2.31065300.487495 1.388883B R Q ⎡⎤
⎢⎥==--⎢⎥⎢⎥-⎣⎦
重复上述过程，迭代11次得12 2.992032 1.000385312.0133920.007496 2.004695 1.94197100.0003250.999895B -⎡⎤
⎢⎥=-⎢⎥
⎢⎥-⎣⎦
1232.992032, 2.004695,0.999895λλλ≈≈≈ 3,2,1精确值
从得出的结果我们可以看出其结果还是很接近精确值的，其实我们完全
可用Matlab 实现上述计算。

用海森伯格QR 基本算法求矩阵全部特征值Matlab 得到运算结果为：
源程序如下：
function l = hessqrtz(A,M)
%海森伯格QR 基本算法求矩阵全部特征值
2222222(0.438310)(0.153846)0.464526, cos 0.4383100.943564, sin 0.1538460.331189100 (3,2)00.9435640.33118900.3311890.943564r r r R ϕϕ=+-====-=-⎡⎤
⎢⎥=-⎢⎥
⎢⎥⎣⎦
%已知矩阵:A
%迭代步数:M
%求得的矩阵特征值:l
A=[5 -3 2;6 -4 4;4 -4 5];
M=11;
A = hess(A);
for(i=1:M)
[q,r]=qr(A);
A = r*q;
l = diag(A);
End
而用QR基本算法求矩阵全部特征值的运行结果如下：
其Matlab源程序如下：
function l = qrtz(A,M)
%QR基本算法求矩阵全部特征值
%已知矩阵:A
%迭代步数:M
%求得的矩阵特征值:l
A=[5 -3 2;6 -4 4;4 -4 5];
M=11;
for(i=1:M) %M为迭代步数
[q,r]=qr(A);
A = r*q;
l = diag(A);
end
由以上两种运算结果的对比可以看出结果基本一致，而且和精确值相当接近，由此可以得出
两种方法均可行。

参考文献
[1] 袁慰平、孙志忠、吴洪伟、闻震出.计算方法与实习.东南大学出版社，2005
[2] 李海涛、邓樱等．MATLAB程序设计教程[M]．北京：高等教育出版社，2004。