普通矩阵特征值的QR算法

H1 A22 H1 ⎥⎦
a1 = (a21 , a31 ,L, an1 )T , a2 = (a12 , a13 ,L, a1n )T ,
A22
=
⎡ ⎢ ⎢
a22 M
L L
a2n M
⎤ ⎥ ⎥
⎢⎣an2 L ann ⎥⎦
只要取 H 1 使得 H 1a1 = (σ1, 0,L, 0)T ，就会使得变换后的矩阵 H1 AH1 的第一列出
则
R(1, 2)B1
=
⎢ ⎢
⎢
b b (2)
(2)
32
33
O
L O
b(2) 3n M
⎥ ⎥ ⎥
=
B2
⎢⎣
b(2) nn−1
b(2) nn
⎥⎦
其中， cosϕ1
=
b(1) 11 r1
,
sinϕ1
=
b(1) 21 r1
,
r1 =
b(1) 11
+
b(1) 21
设
b（2） 32
≠
0
（否则进行下一步），再取旋转矩阵
具体步骤为：
设
b(1) 21
≠
0
（否则进行下一步）
⎡ cosϕ1 sinϕ1 0 L 0⎤
⎢ ⎢
−
sin
ϕ1
cosϕ1
0
L
0⎥⎥
取旋转矩阵 R(1, 2) = ⎢ 0
01
⎥
⎢ ⎢
O
⎥ ⎥
⎢⎣
1⎥⎦
⎡ ⎢
r1
⎢
b(2) 12
b(2) 22
b(2) 13
b(2) 23
L L
b(2) 1n
b(2) 2n
⎤ ⎥ ⎥
设 A=A1 ，对 A1 作 QR 分解，得 A1 = Q1R1，，交换该乘积的次序，得 A2 = R1Q1=， 由于 Q1 正交矩阵，A1 到 A2 的变换为正交相似变换，于是 A1 和 A2 就有相同的特 征值。一般的令 A1=A，对 k=1，2，3，…..
⎧ ⎨ ⎩
Ak = Qk Rk Ak +1 = RkQk
l = k +1
(2) j = k + 1, ..., n
aij = aij − ti u(jk+1)
四、上 Hessenberg 阵的 QR 分解
对上 Hessenberg 阵只需要将其次对角线上的元素约化为零，用 Given 变换比
用 Householder 变换更节省计算量。
1、 平面旋转阵（Givens 变换阵） 定义 n 阶方阵
R−1 i, j
=
RT i, j
平面旋转阵是非对称的正交阵。
（2） RiT, j 也是平面旋转阵。
（3） det(Ri, j )=1
平面旋转阵 Ri,j 的作用：
（1） 将向量 x= ( x1, x2 , ..., xn )T 的第 j 个分量约化为零。 Ri, j 左乘向量 x= ( x1 , x2 , ..., xn )T 只改变 X 的第 i 个分量和第 j 个分量。
⎡1
0⎤
⎢ ⎢
cosϕ2 sinϕ2
⎥ ⎥
⎢ R(3, 2) = ⎢
⎢
－sinϕ2 cosϕ2 1
⎥ ⎥ ⎥
⎢ ⎢
O
⎥ ⎥
⎣
1⎦
⎡ ⎢
r1
⎢
b(3) 12 r2
b(3) 13
b(3) 23
L
b(3) 1n−1
L
b(3) 2n−1
b(3) 1n
b(3) 2n
⎤ ⎥ ⎥
⎢ 则 R(3, 2)B2 = ⎢
（2） 计算 Hk+1 A → A
j = k, k + 1,L, n
∑ （1）t j
（2）i
=
1
n
(
ρk+1 l=k+1
= k + 1,L, n
ul(
k
a +1) lj
)
aij
=
aij
−
t
j
u( k +1) i
（3） 计算
AHk+1 → A
i = 1, ..., n
∑ (1)ti
=
1 ρk+1
n
a u(k+1) il l
首先，选取 Householder 矩阵，使得经相似变换后的矩阵的第一列中有尽可
能多的零元素。为此，应取为如下形式
⎡1 0 L 0⎤
H1
=
⎢⎢0
⎢ ⎢
M
⎥
⎥
H1
⎥ ⎥
⎢⎣0
⎥⎦
其中 H 1 为 n‐1 阶 Householder 矩阵。
于是有
H1 AH1
=
⎡ ⎢
a11
⎢⎣ H 1a1
a2T H 1
⎤ ⎥
b( k +1) k + 2k +1
LL
L h(k +1) k +1n−1
L h(k+1) k + 2n−1
OM b( k +1)
nn−1
b( k +1) 1n M
⎤ ⎥ ⎥
b( k +1) kn
b( k +1) k +1n
⎥ ⎥ ⎥
=
Bk +1
h(k +1) k+2n
⎥ ⎥
M⎥
b( k +1) nn
⎢ ⎢
h21
h22
L
h2n−1
h2n
⎥ ⎥
H=⎢ ⎢ ⎢
h32 h33 O
L O
h3n M
⎥ ⎥ ⎥
⎢⎣
hnn−1 hnn ⎥⎦
的矩阵为上海森堡（Hessenberg）阵。 如果此对角线元 (i=2,3,….n) 全 不 为
零，则该矩阵为不可约的上 Hessenberg 矩阵。
用 Householder 变换将一般矩阵 A 相似变换为 Hessenberg 矩阵。
算法如下：
k = 1, 2, ..., n − 2
（1） 计算
H k +1
=
I
−
1 ρk+1
U (k +1) (U (k +1) )T
n
1
∑ σ k+1 = sign(ak+1，k )( (aik )2 )2 ,
i = k +1
ρk+1 = σ k+1 (σ k+1 + ak+1，k )
U (k+1) = (0, ..., 0,σ k+1 + ak+1，k , ak+2,k , ..., ank )
T
x1, ..., xi-1, r, 0, ..., 0 ,
r = xi2 + L + xn2
（3） 用 Ri, j 对矩阵 A 作变换得到的结论
Ri, j 左乘 A 只改变 A 的第 i，j 行。
RiT, j 右乘 A 只改变 A 的第 i，j 列。
Ri, j ARiT, j 只改变 A 的第 i，j 行和第 i，j 列。
如此进行 n‐2 次，可以构造 n‐2 个 Householder 矩阵 H1, H2 ,L, Hn−2 ，使得
Hn−2 L H2 H1 AH1H2 L Hn−2 = H 。其中 H 为上 Hessenberg 矩阵。特别地，当 A
为实对称矩阵，则经过上述正交变换后，H 变为三对角阵。 用 Householder 方法对矩阵 A 作正交相似变换，使 A 相似于上 Hessenberg 阵，
→
⎢ ⎢ ⎢
⎣
* λ2
L O
*⎤
M
⎥ ⎥
*⎥
λ
n
⎥ ⎦
⎛∗ ∗
⎞
⎜ ⎜
*
O
O
⎟ ⎟
⎜
A（对称阵）− − − − − → 三对角阵B = ⎜
用Householder阵 作正交相似变换
⎜
O∗ * ∗ ∗O
⎟ ⎟ ⎟
⎜ ⎜⎜⎝
O
O *
∗⎟ ∗ ⎟⎟⎠
三、化一般矩阵为 Hessenberg 阵
称形如
⎡ h11 h12 L h1n−1 h1n ⎤
xi2
+
x
2 j
r
xi2
+
x
2 j
r
于是 yi = Cxi + Sx j = r =
xi2
+
x
2 j
,
yj = 0
( ) Ri, j x= x1 , ..., xi-1 , r, xi+1 , ..., xj-1 , 0, xj+1 , ..., xn T
（2） 将向量 x= ( x1, x2 , ..., xn )T 的第i+1 个分量到第 n 个分量约化为零。
若令 y = Ri, j x ，有
⎧ ⎪ ⎨
yi yj
= =
xi cosθ + x j sinθ − xi sinθ + x j cosθ
⎪⎩ yk = xk ,k = 1, ..., n;k ≠ i, j
调整θ ，可将 y j 约化为零。
令
yj
=
0 ，得 tanθ
=
xj xi
所以，取 C = cosθ = xi = xi ， S = sinθ = x j = x j
⎡ ⎢
r1
⎢
⎢
Bk
=
R(k
+ 1, k )Bk−1
=
⎢ ⎢
⎢
⎢
⎢
⎢⎣
b( k ) 1k −1
O rk −1
b( k ) 1k
b( k ) k −1k
b( k ) kk
b( k ) k +1k
b(k ) 1n−1
b(k ) k −1n−1
L
h(k ) kn−1
L b(k) k +1n−1
OM
b( k ) nn−1
普通矩阵特征值的 QR 算法
摘要
求矩阵的特征值有多种不同的办法，本文主要介绍用 QR 法求矩阵的特征 值，QR 法是目前求中等大小矩阵全部特征值的最有效方法之一，使用于求实矩 阵或复矩阵的特征值，它和雅可比法类似，也是一种变换迭代法。
关键词：QR 分解 迭代序列
特征值
Matlab
一 、QR 方法的理论：
现 n‐2 个零元。 同理，可构造如下列形式 Householder 矩阵
⎡1 0 0 L 0⎤
⎡* * * L *⎤
⎢⎢0 1 0 L 0⎥⎥
⎢⎢* * * L *⎥⎥
H2
=
⎢0 ⎢
0
⎢M M
H2
⎥ ⎥
使得
H2 H1 AH1H2
=
⎢ ⎢
⎥
⎢
* * L *⎥
M
O
M
⎥ ⎥
⎢⎣0 0
⎥⎦
⎢⎣
* * *⎥⎦
⎢
b(3) 33
b(3) 43
L L
b(3) 3n−1
b(3) 4n−1
b(3) 3n
b(3) 4n
⎥ ⎥ ⎥
=
B3
⎢ ⎢
OM
M
⎥ ⎥
⎢⎣
b(3) nn−1
b(3) nn
⎥⎦
其中 cosϕ2
=
b(2) 22 r2
,
sin ϕ 2
=
b(2) 32 r2
,
r2 =
假设上述过程已进行了 k‐1 步，有
(b2(22) )2 + (b3(22) )2 。
⎥⎦
因此，最多做 n‐1 次旋转变换，即得
⎡ ⎢
r1
⎢
b(n) 12 r2
b(n) 13
b(n) 23
L L
b(n) 1n
b(n) 2n
⎤ ⎥ ⎥
H (n)
=
R(n,
n
−
1)R(n
−
2,
n
−
1)L R(1, 2)B
=
⎢ ⎢
⎢
r3
L
b(n) 3n
⎥ ⎥
:⎟
:
⎟ ⎟
⎜ ⎜⎜⎝
O
O *
:⎟ ∗ ⎟⎟⎠
Householder变换：如果 v 给出为单位向量而 I 是单位矩阵，则描述上述线性变
换的是 豪斯霍尔德矩阵 （v * 表示向量 v 的共轭转置）H=I ‐2VV*
第二步
⎡λ 1 ⎢
B
−−−−−−−→
用Given变换产生迭代序列
⎧ ⎨ ⎩
Bk = Qk Rk Bk+1 = RkQk
b(k ) 1n
⎤ ⎥
⎥
b( k ) k −1n
h(k ) kn
⎥ ⎥ ⎥
b( k ) k +1n
⎥ ⎥
M⎥
b(k ) nn
⎥⎦
设 b(k) k +1k
≠
0 ，取
⎡1
⎢ ⎢
O
⎢
1
R(k
+
ຫໍສະໝຸດ Baidu
1, k)
=
⎢ ⎢
⎢
⎢
⎢
⎢⎣
cos ϕ k − sinϕk
sin ϕ k cos ϕ k
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ 1⎥ O ⎥⎦
对任意一个非奇异矩阵（可逆矩阵）A，可以把它分解成一个正交阵 Q 和一 个上三角阵的乘积，称为对矩阵 A 的 QR 分解，即 A=QR。如果规定 R 的对角元 取正实数，这种分解是唯一的。若 A 是奇异的，则 A 有零特征值。任取一个不 等于 A 的特征值的实数μ，则 A‐μI 是非奇异的。只要求出 A‐μI 的特征值和特 征向量就容易求出矩阵 A 的特征值和特征向量，所以假设 A 是非奇异的，不是 一般性。
(QR分解) (迭代定义)
这样，可得到一个迭代序列{Ak}，这就是 QR 方法的基本过程。
二、QR 方法的实际计算步骤
⎛ ∗ ... ... ... ... ∗⎞
第一步
⎜ ⎜
*
⎜
A− − − − − → 上Hessenberg阵B = ⎜
用Householder阵 作正交相似变换
⎜
O O
∗ ∗
∗
:
⎟ ⎟
2、用 Givens 变换对上 Hessenberg 阵作 QR 分解
对上
Hessenberg
阵
B
=
⎡⎢⎢bb12((1111)) ⎢
⎢
⎢⎣
b(1) 12
b(1) 22
O
L
O b(1)
nn−1
b(1) 1n
b(1) 2n
⎤ ⎥ ⎥
M b(1)
nn
⎥ ⎥ ⎥⎦
通常用 n‐1 个 Givens 变换阵可将它化成上三角矩阵，从而得到 B 的 QR 分解 式。
其中 cosϕk
=
b(k ) kk rk
,
sin ϕ k
=
b(k ) k +1k rk
，
rk
=
(bk(kk ) )2 + (bk(k+)1k )2 .
⎡ ⎢
r1
⎢
⎢
于是
R(k
+ 1, k)Bk
=
⎢ ⎢
⎢
⎢
⎢
⎢⎣
b( k +1) 1k
O rk
b( k +1) 1k +1
b( k +1) kk +1
b( k +1) k +1k +1
( ) Ri,i+1 x=
T
x1 , ..., xi-1, r, 0, xi+2 , ..., xn ,
r=
xi2
+
x2 i +1
( ) Ri,i+2 Ri,i+1 x=
x1 , ..., xi-1 , r, 0, 0, xi+3 , ..., xn
T
,
r=
xi2
+
x2 i +1
+
x2 i+2
( ) Ri,n L Ri,i+2 Ri,i+1 x=
⎡1
⎤
⎢ ⎢
O
⎥ ⎥
⎢
1
⎥
⎢ ⎢
cosθ
sinθ
⎥ ⎥
−−i
⎢
1
⎥
⎢ Ri, j = ⎢
O
⎥ ⎥
⎢ ⎢
1
⎥ ⎥
⎢
− sinθ
cosθ
⎥ −− j
⎢ ⎢
1
⎥ ⎥
⎢
O⎥
⎢ ⎣
1⎥⎦ ( j > i)
称为平面旋转阵，或称为 Givens 变换阵。
平面旋转阵 Ri,j 的性质：
（1）
RiT, j Ri, j = I ,