-
第
9
章
矩阵特征值的数值解法
9.1
引言
矩阵特征值问题有广泛的应用背景
.
例如动力系统和结构系统中的振动问题、
电力系统
的静态稳定分
析上、工程设计中的某些临界值的确定等,都归结为矩阵特征值问题
.
数学中
诸如方阵的对角化及解微分方程组等问题,都要用到特征值的理论
.
本章介绍
n
阶实
矩阵
A
?
R
n
?
n
的特征值与特征向量的数值解法<
/p>
.
定义
9.1.1
已知
n
阶实
矩阵
A
?
(
a
ij
)
?
R<
/p>
n
?
n
,如果存
在常数
?
和非零向量
x
,使
Ax
?
?
x
或
p>
(
A
?
?
I
)
x
?
0
(9.1.1)
那么称
?
为
A
的
特征值
(eigenvalue)
,
< br>x
为
A
的相应于
?
的
特征向量
(eigenv
ector)
.
多项式
?
a
11
?
?
?
a
p
< br>n
(
?
)
?
det(
A
?
?
I
)
?
?
21
?
?
?
?
a
n
1
p>
称为
特征多项式
(characteris
tic polynomial)
,
a
1
n
?
p>
a
22
?
?
?
a
2
n
?
?
(9.1.2)
?
?
?
?
?
a
n
2
?
a
nn
?
?
?
?
a
12
<
/p>
det(
A
?
?
I
)
?
0
p>
(9.1.3)
称为
特征方程
(characteristic equation)
.
注
式
p>
(9.1.3)
是以
?
为未知量的一元
n
次代数方程,
p
n
(
?
)
p>
?
det(
A
?<
/p>
?
I
)
是
?
的
n
次
多项式
.
显然,
A
p>
的特征值就是特征方程
(9.1.3)
的根
.
特征方程
(9.1.3)
在复数范围内恒有解,
其个数为方程的次数
(<
/p>
重根按重数计算
)
,
因此
n
阶矩阵
A
< br>在复数范围内有
n
个特征值
.
除特
殊情况
(
如
n
?
2,
3
或
A
为上
(
下
)
三角矩阵
)
外,一般不通过直接求解特征方程
(9.1.3)
来求
A
的特征值
,
原因是这样的算法往往不稳定
.
在计
算上常用的方法是幂法与反幂法和相似变换
方法
.
本章只介绍求矩阵特征值与特征向量的这两种基本方法
.
p>
为此将一些特征值和特征向
量的性质列在此处
.
定理
9.1.2
设
n
阶方阵
A
?
(
a
p>
ij
)
n
?
n
的特征值为
?
1<
/p>
,
?
2
,
?
,
?
n
,那么
?
1
?
?
2
?
?
?
?
n
?
a
11
?
a
22
?
?
?
a
nn
;
(2)
?
1
?
2
?
?<
/p>
n
?
det
A<
/p>
.
定理
9
.1
.3
如果
?
是方阵
A
的特征值,那么
k
k
(1)
?
是
A
p>
的特征值,其中
k
是正整数
;
(1)
?
1
(2)
当
A
是非奇异阵时,
< br>1
是
A
的特征值
.
?
(3)
p
n
(
?
)
是
p
n
(
A
)
的特征值,其中
p
n
(
x
)
是多项式
p
n
(
x
)
< br>?
a
0
?
a
1
x
?
a
2
x
2
?
p>
?
?
a
n
x
n
.
定义
9.1.4
设
A
,
p>
B
都是
n
阶方阵<
/p>
.
若有
n
阶非
奇异阵
P
,使得
P
AP
?
B
,则称矩
阵
A
与
B
相似
(similar)
,
P
p>
AP
称为对
A
进行
相似变换
(similarity transformati
on)
,
P
称
1
?
1
?<
/p>
1
为
相似变换
矩阵
(similarity transformation
matrix)
.
定理
9.1.5
若矩阵
A
与
B
相似,则
A
与
B
的特征值相同
.
定理
9.1.6
如果
A
是<
/p>
n
阶正交矩阵,那么
(1)
A
?
A
,且
det
A
?
1
或
?
1
;
(2)
若
y
?
p>
Ax
,则
y
2
p>
?
1
T
?
x
2
,
即
x
T
?
x
?
y
T
?
y
.
定理
9.1.7
设
A
是任意
n
阶实对称矩阵,则
(1)
A
的特征值都是实数;
(2)
A
有
n
个线性无关的特征向量
.
定理
9.1.8
设
A
是任意
n
阶实对称矩阵,则必存在
n
阶正交矩阵
P
,使得
P
?
1
AP<
/p>
?
P
T
AP
p>
?
?
,
其中
?
?
diag(<
/p>
?
1
,
?
2
,
?
,
?
n
)
是以
A
的
n
个特征值
?
1
,
?
2
,
?
,
?
n
为对角元素的对角矩阵
.
定理
9.1.9
(
圆盘定理
)
矩阵
A
?
(
a
ij
)
n
?
n
的任意一个特征值至少位于复平面上的几个圆
盘
?
D
i
?
?
z
z
?
a
ii
?
?
中的一个圆盘上。
j
?
1,
j
< br>?
i
?
n
?
a
ij
?
,
i
?
1,2,
?
,
n
,
?
9.2
幂法与反幂法
9.2.1
幂法及其加速
9.2.1.1
幂法
幂法是计算矩阵按模最大特征值
(largest
eigenvalue in magnitude)
及相应特征向量的迭代
法
.
该方法稍加修改,也可用来确定其他特征值
.
幂法的一个很有用的特性是:它不仅可以
求特征值,
而
且可以求相应的特征向量
.
实际上,
幂法经常用来求通过其他方法确定的特征
值的特征向量
.
下面探讨幂法的具体过程
.
设矩阵
A
?
R
n
?
n
的
< br>n
个特征值满足
?
1
?
?
2
?
?
3
?
?
?
?
n
?
0
,
(9.2.1)
且有相应的
n
个线性无关的特征向量
x
< br>1
,
x
2
,
?
,
x
n
,则
x
1
,<
/p>
x
2
,
?
,
x
n
构成
n
维向量空间
R
n
p>
?
?
的一组基
,
因此
R
?
?<
/p>
z
z
?
?
?
i
x
i
,
?
i
?
R
,
i
?
1,2,
?
,
n
< br>?
.
i
?
1
?
?
n
n
在
R<
/p>
中选取某个满足
?
1
?
0
的非零向量
2
n
<
/p>
z
0
?
?
?
i
x
i
.
i
?
1
n
用矩阵
A
同时左乘上式两
边,得
Az
0
?
?
?
i
A
x
i
?
?
?<
/p>
i
?
i
x
i
.
i
?
1
i
?
1
n
n
再用矩阵
A
左乘上式两边,得
A
z
0
?
?
?
p>
i
?
i
2
x
i
.
2
i
?
1
n
这样继续下去,一般地有
A
z
0
?
?
?<
/p>
i
?
i
k
x
i
,
k
?
1,2,
?
.
(9.2.2)
k
i
?
1
n
记
z
k
?
Az
k
?
1
?
A
k
z
0
,
p>
n
k
?
1
,2,
?
,则由式
(9
.2.2)
得
k
n
?
?
?
?
?
k
k
k<
/p>
i
z
k
?
A
z
0
?
?
?
i
?
i
x
i
?
?
1
?
?
1
x
1
?
?<
/p>
?
i
?
?
x
i
?
,
k
?
1,2,
?
.
(9.2.3)
?
i
?
1
i
?
2
?
?
1
?
?
?<
/p>
?
由假设
(9.2.1)
,结合式
(9.2.3)
,得
lim
k
??
z
k
k
?
1<
/p>
?
?
1
x
1
,
(9.2.4)
于是对充分大的
k<
/p>
有
z
k
?
?
1
k
?
1
x
1
.
(9.2.5)
k
式
(9.2.4)
表明随着
k
的增大,
序列
z
k
?
1
越来越接近
< br>A
的对应于特征值
?
1
的特征向量
x
1
的
p>
?
?
?
1
倍
,
由此可确定对应于
?
1
的特征向量
x
< br>1
.
当
k
充分大时,可得
x
1
的近似值
p>
.
上述收敛速度取决于比值
?
2
?
1
.
事实上,由式
(9.2.3)
知,
k
?
2
?
?
x
?
?
1
1
2
p>
?
1
k
?
1
z
k
再由式
(9.2.1)
得
?
x
2
?
p>
?
3
3
?
1
k
?
x
3
?
?
?
< br>?
n
n
?
1
k
x
n
.
(9.2.6)
1
?
?
?
2
?
3
?
?
?
?
n
.
(9.2.7)
?
< br>1
?
1
?
1
k
结合式
(9.2.6)
和式
(9.2.7)
知,序列
< br>z
k
?
1
收敛速度取决于比值
?
?
?
2
?
1
.
p>
下面计算
?
1
.
由式
(9.2.3)
知
3
< br>k
?
1
n
?
?
?
?
?
i
k
?
1
p>
k
?
1
z
k
?
1
?
A
z
k
?
< br>A
z
0
?
?
1
?
?
1
x
1
?
?
p>
?
i
?
?
x
i
?
,
?
?
i
< br>?
2
?
?
1
?
?
?
k
?
1
当
k
p>
充分大时
,
z
k
?
1
?
?
p>
1
?
1
x
1
.
结合式
(9.2
.5)
,得
z
k
?
1
?
?
1
z
k
. <
/p>
这表明两个相邻向量大体上只差一个常数倍,这个倍数就是
A
p>
的按模最大特征值
?
1
.
记
(1)
(2)
(
n
)
z
< br>k
?
?
z
k
,
z
k
,
?
,
z
k
p>
?
,
则有
p>
T
(
j
)
z
k
1
lim
(
?
?
?
1
,
j
k
< br>??
z
)
k
j
?
1,2,
?
< br>,
n
,
(9.2.8)
即两个相邻的迭代
向量所有对应分量的比值收敛到
?
1
.
定义
9.2.1
上述由已知非零向量
z
0
及矩阵
A
的乘幂
A
构造向量序列
{
z<
/p>
k
}
来计算
A<
/p>
的
按模最大特征值
?
1
及相应特征向量的方法称为
幂法
(power
method)
,其收敛速度由比值
k
?
?
?
2
?
1
来确定,
?
越小,收敛越快
.
注
由幂法的迭代过程
(9.2.3)
容易看出,如果
,那么迭代向量
z
k
的
?
1
?
1
(或
?
1
?
1
< br>)
各个非零的分量将随着
k
?<
/p>
?
而趋于无穷
(或趋于零)
,
这样在计算机上实现时就可能上溢
(
或下溢
).
为了克服这个缺点,需将每步迭代向量进行规范化
:
记
y
k
p>
?
Az
k
?
1
?
?
y
,
y
,
?
,
y
(1)
k
(2)
k
(
n
)
T
k
?
.
若存在
y
k
< br>的某个分量
y
k
(
j
0
)
,满足
y
k
(
j
< br>0
)
(
j
)
(
j
)
,
则记
max(
y
k
)
?
y
k
0
.
将
y
k
规范化
?
max
y
k
1
?
j
?
n
z
k
p>
?
y
k
max(<
/p>
y
k
)
,
p>
这样就把
z
k
的分
量全部控制在
[
?
1,1]
中
.
例如,设
y
k
?
(
?
2,
3,
0,
?
5,
1)
T
,因为
y
k
的所有分量
中,绝对值最大的的是
?
5
,
所以
max(
y
k
p>
)
?
?
5
,故
z
k
?
y
k
max(
y
k
)
?
(0.4,
p>
?
0.6,
0,
1
,
?
0.2)
T
.
综上所述,得到下列算法:
算法
9
.2.1
(
幂法
)
设
A
是
n
阶实矩
阵,取初始向量
z
0
?
R
n
,通常取
z
0
?
(1
,1
,
?
,1)
T
,
其迭代过程是:对
k
?<
/p>
1,2,
?
,有
?
y
k
?
p>
Az
k
?
1
,
?
?
m
k
?
max(
y
k
),
(9.2.9)
?
z
?
y
m
.
k
k
?
k
定理
9.2.1
< br>对式
(9.2.9)
中的序列
?
z
k
?
和
p>
?
m
k
?
有
lim
z
k
?
k
??
x
1
,
lim
m
k
?
?
1
p>
,
(9.2.10)
k
??
max
?
x
1
?
其收敛速度由
?
?
p>
?
2
?
1
确定
.
4
证明
由
迭代过程
(9.2.9)
知
z
k
?
y
k
m
k
?
< br>Az
k
?
1
m
k
?
A
y
k
?
1
m<
/p>
k
m
k
?
1
?
A
2
z
k
?
1
m
k
m
k
?
1
?
?
?
A
z
0<
/p>
n
k
?
m
?
A
i
i
?
1
k
k
z
0
max(
A
k
z
0
)
< br>,
(9.2.11)
k
n
?
?
?
?
< br>?
k
k
i
其中
z
0
?
?
?
i
x
i<
/p>
.
若
?
1
p>
?
0
,则由
(9.
2.3)
知:
A
z
0
?
?
1
?
?
1
x
1<
/p>
?
?
?
i
?
?
x
i
?
,
代入式
i
?
1
?
i
?
2
?
?
< br>1
?
?
?
?
(9.2.11)
得
?
n
?
?
?
k
1
x
1
?
?
?
i
?
i
z
i<
/p>
?
2
?
?
?
x
i
1
?
k
?
,
p>
n
k
max
?
p>
?
?
?
?
?
?
?
1
x
1
?
?
< br>?
?
i
?
i
?
x
i
?
i
?
2
?
p>
?
1
?
?
?
故
lim
x
1
k
??
z
k
?
max
?
x
.
1
?
而
p>
y
Ay
2
k
?
1
A
z
k
?
2
A
k
z
0
k
?
Az
k
?
1
?
max
?
y
?
?
?
?
k
?
1
p>
k
?
1
?
max
?
A
z
k
?
2
?
max
?
A
z
0
?
?
k
< br>?
n
a
?
?
k
i
?
?
1
?
1
x
p>
1
?
?
a
?
?
p>
i
?
?
x
i
?
?
i
?
2
?
?
< br>1
?
?
?
,
?
k
?
1
?<
/p>
n
?
?
p>
k
?
1
p>
i
?
1
m
a
x
?
?
a
x
?
?
< br>?
a
?
1
1
i
?
?
?
?
x
i
?
p>
i
?
2
1
?
?
?
于是
max
?
n
?
?
a
?
?
k
i
?
x
?
1
x
1
?
?
a
i<
/p>
?
?
i
?
m
?
i
?
2
?
?
1
?
?
?
k
?
max(
y
k
< br>)
?
?
1
,
p>
max
?
n
k
p>
?
1
?
?
a
?
?
i
?
?
1
x
< br>1
?
?
a
i
?
?
x
i
?
?
i
?
p>
2
?
?
1
?
?
?
故
lim
k
??
m
k
?
?
1
.
p>
由式
(9.2.6)
和式
< br>(9.2.7)
知:上述收敛速度由
?
< br>?
?
2
?
1
确定
.
证毕!
例
9.2.1
用幂法求方阵
(9.2.12)
(9.2.13)
(9.2.14)
(9.2.15)
(9.2.16)
5
p>
?
1
2
3
?
?
A
?
?
2
1
< br>3
?
?
?
?
3
3
5
?
?
按模最大特征值及相应的特征向量,要求
m
k
?
m
k
?
1
?
10
?
3
.
解
选取初
始向量
z
0
?
(1
,1
,1)
T
,按式
(9.2.9)
迭代,结果见表
9.2.1.
表
9.2.1
k
1
2
3
4
5
z
k
(0.
545455,
0.545455,
1)
T
(0.560441
,
0.560441
,
1)
T
m
k
11
8.2727
m
k
?
m
k
?
1
2.7273
(0.559787,
0.559787,
1)
T
8.3627
(0.
559818,
0.559818,
1)
T
8.3587
(0.559817,
0.559817,
1)
T
8.3589
9
?
10
?<
/p>
2
4
?
10
?
3
0.2
?
10
?
3
因此,所求按模最大特征值
?
1
?
8.3589
,相应特征向量
x
1
?
(0.559817,
0.559817,
< br>1)
T
.
< br>事
实
上
,
A
的
按
模
最
大
特
征
值
p>
?
1
?
4
?
19
?
8.3588
989
?
,
相
应
特
征
向
量<
/p>
x
1
?
(0.5
5981649
?
,
0.559816
49
?
,
1)
T
,
故所得结果具有较高的精度
.
9.2.1.2
幂法的加速
从上面的讨论可知,由幂
法求按模最大特征值,可归结为求数列
{
m
k
}
的极限值,其收
敛速度由
p>
?
?
?
2
?
1
确定
.
当
?
?
?
2
?
1
接近
1
时
,
收敛速度相当缓慢
.
为了提高收敛速度
,
?
2
?
1
确定,
所以若
{
m
k
}
收敛,
当
k
充分大时,
则有
可以利用外推法进行加速
.
因为序列
{
m
k
}
的收敛速度由
?
?
?
?
m
p>
k
?
?
1
?
O
?
?
2
?
1
?
< br>?
k
?
?
,或改写为
?
m
< br>k
?
?
1
?
?
c
?
2
1
?
?
k
p>
,
其中
c
是与<
/p>
k
无关的常数
.
由此可得
m
k
?
1
?
?<
/p>
1
?
2
?
,
(9.2.17)
m
k
?
?
< br>1
?
1
这表明幂法是线性收敛的
.
由式
(9.2.17)
知
m
k
?
1
?
?
1
m
k
?
2
?
?
1
.
?
m
k
?
p>
?
1
m
k
?
1
?
?
1
6
?
k
?
2
,即
由上式解出
?
1
,并记为
m
?
k
?
2
m
2
m
k
?
2
m
k
?
m
k
(
m
k<
/p>
?
1
?
m
k
)
2
?
1
,
(9.2.18)
?
?
m
k
?
m
< br>k
?
2
?
2
m
k
?
1
?
m
k
m
p>
k
?
2
?
2
m
k
?
1
?
m
k
< br>这就是计算按模最大特征值的加速公式
.
将上面的分析归结为如下算法:
算法
9
.2.2
(
幂法的加速
)
设
A
是
n
阶实矩阵,给定非零初始向量
z
0
?<
/p>
R
n
,通常取
T
.
对
k
?<
/p>
1,
2
,用迭代式
z
0
?
(
1,
1,
?
,
1)
?
y
k
?
Az
k
?
1<
/p>
,
?
?
m
k
?
max(
y
p>
k
),
?
z
?
y
m
k
k
?
k
求出
m
1
,
< br>m
2
及
z
1
,
z
2
.
再对
k
?
3,
4,
?
,迭代过程为
?
y
k
?
Az
k
?
1
,
?
m
?
m
ax(
y
),
k
k
?
?
(
m
k
?
1
?
p>
m
k
?
2
)
2
(9.2.19)
?
?
k
?
< br>m
k
?
2
?
,
?
m
m
k
?
2
m
p>
k
?
1
?
m
k
?
2
?
?
k
.
< br>?
?
z
k
?
y
k
m
?
k
?
m
?
p>
k
?
1
?
?
(
?
?
0
是预先给定的精度
)
时
,迭代结束,并计算
z
k
;否则继续迭
代,
当
m
?
k
?
m
?
k
p>
?
1
?
?
.
直至满足迭代停止条件
m
< br>有关加速收敛技术,读者请参考文献
[11].
8.2.2
反幂法及其加速
反幂法是计算矩阵按模最小特征值及相应特征向量的迭代法
,
其基本思想是:设矩阵
?
1
?
1
A
?
R
n
?
n
非奇异,
用其逆矩阵
A
代替
A
,
矩阵
A
的按模最小特征值就是矩阵
A
的按模
最大特征值
.
这样用
A
代替
A
做幂法,
即可求出
A
的按模最大特征值,
也
就是矩阵
A
的
按模最小特征值;这种方
法称为
反幂法
(inverse power
method)
.
因为矩阵
A
非奇异,所以由
Ax
i
?
?
i
x
i
可知:
A
x
i
?
值满足
?
1
?
1
p>
?
1
1
?
i
x
i
.
这说明:如果
A
的特征
?
1
?
?
2
p>
?
?
?
?
n-
1
?
?
n
?
0
,
(9.2.20)
那么
A
的特征值满足
p>
?
1
1
?
n
?
1
?
n-
1
?
?
?
1
?
2
?
1
?
1
,
(9.2.21)
7
且对
A
的应
于特征值
?
i
的特征向量
x
i
也是
A
的对应于特征值
1
?
i
的特征向量
.
由上述分析知:对
A
应用幂法求按模最大的特征值
1
是求
A
的按模最小的特征值
?
n
及相应的特征向量
x
n
.
算法
9
.2.3
(
反幂法
)
任取初始非零向量
z
0
?
R
n
,通常取
z
0
?
(1
,1
,
?
,1)
T
.
为了避免
求
A
p>
,对
k
?
1,2,
?
,将迭代过程
(9.2.9)
改写为
:
?
p>
1
?
1
?
1
?
n
及相应的特征向
量
x
n
,就
?
Ay
k
?
z<
/p>
k
?
1
,
?
?
m
k
?
max(
y
k
),
(9.2.22)
?
z
?
y
m
.
< br>k
k
?
k
仿定理
9.2.1
的证明,可得:
定理
9.2.2
对式
(9.2.22)
中的序列
?
z
k
?
和
?
m
k
?
有
lim
z
k
?
k
??
x
n
,
max(
x
n
)
lim
m
k
?
k
??
1
?
n
,
(9.2.23)
?
?
其收敛速度由
?
?
n
?
n
?
1
确定
.
注
按
(9.
2.22)
进行计算,每次迭代都需要解一个方程组
Ay
k
?
z
k
?
1
.
若利用三角分解
法
求解方程组,即
A
?
LU
,其中
L
是下三角矩阵,
U
是上三角矩阵,这样每次迭代只需解
两个三角方程组
?
L
?
?
z
k
?
1
,
?
Uy
?
v
.
?
k
?
1
p>
9.2.3
原点平移法
为了提高收敛速度,下面介绍加速收敛的原点平移法
.
设矩阵
B
?
A
?
p
I
p>
,其中
p
是一个待定的常数,
A
与
B
除主对角线上的元素
外,其
他元素相同
.
设
A
的特征值为
?
1
,
?
2
,
?
,
?
n
< br>,则
B
的特征值为
?
1
?
p
,
?
2
?
p
,
?
,
?
n
?
p
,
且<
/p>
A
与
B
的特征向
量相同
.
9.2.3.1
原点平移法的幂法
设
?
1
是
A
的按模最大的特征值,选择
p
,
使
?
1
?<
/p>
p
?
?
2
?
p
?
?
i
?
p
,
i
?
3,4,
?
,
n
,
及
?
2
p>
?
p
?
2
?
.
(9.2.24)
?
1
?
p
< br>?
1
对
B
应用幂法,得
算法
9
.2.4
对
k
?
p>
1,2,
?
,有
8
?<
/p>
y
k
?
(
A
?
p
I
)
z
k
?
1
,
?
?
m
k
?
max(
< br>y
k
),
(9.2.25)
?
z
?
y
m
,
< br>k
k
?
k
且
lim
m
k
?
?
1
?
p
,
lim
z
k
?
k
??<
/p>
k
??
x
1
p>
,
(9.2.26)
max
?
x
1
?
其收敛速度由<
/p>
(
?
2
?
p
)
(
?
1
?
p
)
确定
.
由式
(9.2.24)
知:在计算
B
的按模最大特征值
?
1
?
p
的过程
(9.2.25)
中,收敛速度得到
加速;算法
9.2.4
< br>又称为
原点平移下的幂法
(power method
with shift)
.
9.2.3.2
原点平移下的反幂法
设
?
n
是
A
< br>的按模最小的特征值,选择
p
,使
?
n
?
p
?
?
n
?
p>
1
?
p
?
?
i
?
p
,
i
?
1
< br>,2,
?
,
n
< br>?
2
.
(9.2.27)
及
?
n
?
p
< br>?
?
n
.
(9.2.28)
?
n
?
1
?
p
< br>?
n
?
1
若矩阵
B
?
A
?
p
I
可逆,则
B
的特征值为
?
1
1
1
,
,
?
,
?
1
?
p
?
2<
/p>
?
p
且有
1
,
?
p>
n
?
p
1
1
1
?
?
,
i
?
1,2,
?
,
n
?
2
.
(9.2.29)
?
n
?
p
?
< br>n
?
1
?
p
?
i
?
p
对
B
应用反幂法,得:
算法
9
.2.5
对
k
?
p>
1,2,
?
,有
?
(
A
?
p>
p
I
)
y
k
?
z
k
?
1
,
?
< br>?
m
k
?
max(
y
k
),
< br>
(9.2.30)
?
z
?
y
m
,
< br>k
k
?
k
且
lim
m
k
?
k
??
x
1
1
,
li
m
z
k
?
,<
/p>
(9.2.31)
?
n
?
p
k
??
max
?
x
1
?
其收敛速度由
(
?
n
?
p
)
(
?
n
?
< br>1
?
p
)
确定
.
9
由式
(9.2.28)
知:在计算
B
的按模最大特征值
?
1
1
的过程
(9.2.30)
中,收敛速度得
?
n
?
p
到加速
.
算法
9.2.5
又称为
原点平移下的反幂法
(inverse
power method with shift)
.
定义
9
.2.8
原点平移下的幂法与原点平移下的反幂法统称为
原点平移法
.
注
有的资料上的原点平移法专指原点
平移下的反幂法;而有的资料上的反幂法指的
就是原点平移下的反幂法
< br>.
原点平移下的反幂法
(
算法
9.2.7)
的主要应用是:已知矩阵的近似特征值后,求矩阵
的
?
,则
A
?
?
?
I
的特征
特征向量
.
其主要思想是:若已知<
/p>
A
的特征值
?
m
的近似特征值为
?
m
< br>m
?
(
i
?
1
,2,
?
,
m
)
,其中
?
(
i
?
1
,2,
?
,
m
)
是
A
的特征值
.
而按模最小的特征
值就是
?
p>
i
?
?
i
m
?
,相应的特征向量与
A
的特征向量相同
.
值是
?
m
?
?
m
利用公式
(9.2.3
0)
可求出
?
m
k
?
,
?
z
k
?
,并且有
lim
m
k
?
k
??
1
,<
/p>
?
?
m
?
?
m
lim
z
k
?
k
??
x
m
,
max(
x
m
)
?
?
m
?
?
m
?
?
min
?
?
?
?
)
确定
.
于是当
k
充分大时,可取
(
?
i
?
?
其收敛速度由
m
i
p>
m
?
1
?
i
?
m
?
1
?
i
?
< br>?
m
0
0
?
?
1
,
x
m
?
z
k<
/p>
,
?
m
?
?
m
m
k
?
适当好,收敛速度就很快,往往迭代几次就能
得到满意的结果
.
只要近似值
?
p>
m
?
?
?
0.3589
,用原点平移下的反幂法求方阵
例
9.2.2
已知特征值
?
的近似值
?
?
1
2
3
?
?
A
?
?
2
p>
1
3
?
?
?
?
3
3
5
?
?
得对应特征值
?
的特征向量
.
?
?
?
0.3589
,对矩阵
解
取
p
?
p>
?
2
3
?
?
1.3589
?
.
A
?
pI
?<
/p>
A
?
0.3589
I
?
?
2
1
.3589
3
?
?
?
3
5.3589
?
?
3
?
迭代公式
(9.2.30)
中的
y
k
是通过解方程组
(
A
?
p
I
)
y
k
?
z
k
?
1
求得的
.
为了节省工作量,可先将<
/p>
A
?
p
I
进行
LU
分解
. <
/p>
在
LU
分解中尽量避免较小的
u
rr
当除数,
通常可以
先对矩阵
A
?
p
I
的行进行调换后再
10
?
0
0
p>
1
?
?
?
分解
.
为此,我们可用
P
?
0
1
0
乘
A
?
p
p>
I
后再进行
LU
分
解,即
?
?
?
?
1
0
0<
/p>
?
?
2
3
?
?
3
3
5.3589
?
?
0
p>
0
1
?
?
1.3589
?
?
2<
/p>
?
?
?
2
?
?
LU
P
(
A
?
p
I
)
?
?
< br>0
1
0
1.3589
3
1.3589
3
?
p>
?
?
?
?
?
?
3
5.3589<
/p>
?
2
3
?
?
1
0
0
?
?
?
?
3
?
?
?
1.3589
?
0
0
?
?
3
3
< br>5.3589
?
?
1
?
?
0
?
0.6411
?
,
?
?
0.6667
1
0<
/p>
?
0.5726
?
?
?
?
?
6
?
0
?
3.0
7
?
10
?
?
0.4530
?
1
1
?
?
?
?
0
?
P
(<
/p>
A
?
1.2679
I
)
y
k
?
Pz
k
?
1<
/p>
,
即
LUy
k
?
Pz
k
?
1<
/p>
.
令
Lv
k<
/p>
?
Pz
k
?
p>
1
,
Uy
k
?
v
k
?
1
.
选取
z
0
,使
Uy
1
?
L
?
< br>1
Pz
0
?
(1
,1
,1)
T
,得
y
1
< br>?
(290929.45,
290927.56,
?
325732.90)
T
< br>,
m
1
?
max(
y
1
)
?
?
325732.90
,
z
1
?
y
1
m
1
?
(
?
0.8932,
?
0.8931
,
1)
T
.
由
Uy
2
?
L<
/p>
?
1
Pz
1
p>
得
y
2
?
(
?
845418.
49,
?
845418.49,
946
558.42)
T
,
m
2
?
max(
y
2
)
?
< br>946558.42
,
z
p>
2
?
y
2
m
2
?
(
?
0.8932,
?
0.8
932,
1
)
T
.
由于
z
1
与
z
2
的对应分量几乎相等,故
p>
A
的特征值为
?
?
?
?
?
p>
1
1
?
?
0.3589
?
?
?<
/p>
0.35889894354117
,
m
2
946558.42
相应的特征向
量为
x
?
z
2
?
(
?
0.8
932,
?
0.8932,
1)
T
.
而矩阵
p>
A
的一个特征值为
?
?
4
?
19
?
0.358898943540674
?
,相应的特征向量为
(
?
0.89
315,
?
0.89315,
1)
p>
T
,由此可见得到的结果具有较高的精度
.
9.3
QR
算法
上一节我们介绍了求矩阵特征值的幂法和反幂法
.
幂法主要用来求矩阵的按模最大特
征值,而反幂法主要用于求特征向量
p>
.
本节将介绍幂法的推广和变形
——
p>
QR
算法,它是求
一般中小型矩阵全部特征
值的最有效的方法之一,其基本思想就是利用矩阵的
QR
分解<
/p>
.
矩
阵
A
p>
?
R
n
?
n
的
QR
分解就是:用
Householder
变换将矩阵
A
分解成正交矩阵
Q
与上三角矩
阵
R
的乘积,即
A
p>
?
QR
.
下面首
先介绍
Householder
变换
.
9.3.1
豪斯霍尔德变换
11
定义
9.3.1
设
B
?
p>
(
b
ij
)
n
?
n
是
n
阶方阵,若当
i
?
p>
j
?
1
时,
b
ij
?
0
,则称矩阵
B
为上
He
ssenberg
矩阵
(
Hessen
berg matrix
)
,
又称为<
/p>
准上三角矩阵
,
它的一般形式为
?
b
11
b
12
?
b
?
21
b
22
b
32
B
?
?
?
?
?
?
矩阵
——
Householde
r
矩阵
.
?
?
?
?
?
?<
/p>
?
b
n
,
n
?
1
b
1
n
?
b
2
n
?
?
b
3
n
?
.
(9.3.1)
?
?
?
b
nn
< br>?
?
下面讨论如何将矩阵
A
p>
用正交相似变换化成
(9.3.1)
的形式
.
为此先介绍一个对称正交
定义
p>
9.3.2
设
向量
u
?
R
的
欧氏长度
u
2
?
1
,
I
为
n
阶单位矩阵,则称
n
阶方阵
n
H
?
H
(
u
)
?
I
?
2
uu
T
(9.3.2)
为
Househol
der
矩阵
(Householder
matrix)
.
对任何
x
?
R
,称由
H
?
H
(
u
)
确定的变换
n
y
?
Hx
(9.3.3)
为
镜
面
反
射
变
换
(specular
reflection
transformation)
p>
,
或
Householder
(Householder
transformation)
.
注
Hou
seholder
变换,最初由
A.C Aitken
在
1932
年提出
.
Alston Scott Householder
在
195
8
年指出了这一变换在数值线性代数上的意义
.
这一变换将一个向量变换为由一个超平
面反射的镜像,是一种线性变换
.
运用线性代数的知识,很容易证明:
定理
9.3.3
(9.3.2)
式定义的矩阵
H
是对称正交矩阵;对任何
x
?
R
,由线性变换
n
y
?
H
x
得到
y
的欧氏长度满足
y<
/p>
2
?
x
2
.
反之,有下列结论:
定理
9.3.4
设
x
,
p>
y
?
R
n
,
x
?
y
.
若
x
反射矩阵
H
(
u
)
,使得
H
x
?
y
.
证
令
u
?
p>
2
?
y
2
,
则一定存在由单位向量确定的镜面
x
?
y
,显然
u
2
?
1
.
构造单位向量
u
确定的镜面反射矩阵
x
?
y
< br>2
H
?
H
(
u
)
?
I
?
2
uu
T<
/p>
,
?
(
x
?
y
)(
x
?
y
)
T
?
2(
x
?
y
)(
x
< br>T
x
?
y
T
x
)
.
?
x
?
x
?<
/p>
Hx
?
?
I
p>
?
2
uu
?
x
?
?
I
?
2
2
2
x
?
y
2
x
?
y
2
?
?
?
?
T<
/p>
又因为
x
2
?<
/p>
y
2
,即
x
p>
T
x
?
y
T
y
,所以
————————————————————
海森伯格
(
Karl
Adolf
Hessenberg
1904
年
9
月
8
日
–<
/p>
1959
年
2
月
22
日
)
是德国数学家和工程师
.
豪斯霍尔德
(
Alston
Scott
Householder
1904
年
5
月
5
日
–
1993
年
7
月
4
日
)
是美国数学家,
在数学生物学和
数值分析等领域卓有建树
.
12
2
x
?
y
2
?
p>
(
x
?
y
)
T
(
x
?
y
)
?
< br>(
x
T
?
y
T
)(
x
?
y
)
?
x<
/p>
T
x
?
x
T
y
?
y
T
x
?
y
T
y
?
x
T
x
?
x
T
y
?
x
T<
/p>
y
?
x
T
x
?
2(
x
T
x
?
y
T
x
),
于是
Hx
?
x
?
证毕!
由定理
9.3.4
得:
p>
2(
x
?
y
)(
x
T
x
?
y
T
x
)
x
?
y
< br>2
2
?
x
?
(
x
?
y
)
?
y
.
算法
9.3.1
若
x
?
p>
(
x
1
,
x
2
,
?
,
x
n
)
< br>T
,其中
x
2
< br>,
?
,
x
n
不全为零,则由
?
?
?
sgn(
x
1
)
x
2
,
?
T
n
?
u
?
x
?
?
e
1
,<
/p>
其中
e
1
?
p>
(1,0,
?
,0)
?
R
,
?
1
2
(9.3.4)
?
?
=
2
u
2
?
?
(
?
?
x
1
),
?
T
?
H
?
p>
H
(
u
)
?
I
?
2
uu
?
I
?
?
?
1
uu
< br>T
2
?
u
?
2
?
1,
a
?
0,
?
a
?
0,
确定
的镜面反射矩阵
H
,使得
Hx
?
?
e
1
,其中
sgn(
a
)
?
?
0,
?
?
1,
a
?
0.
?
例
9.3.1
p>
设
x
?
(
?
1,
2,
?
2)
T
,
按
式
(9.3.4)<
/p>
的
方
法
构
造
镜
面
反
射
矩
阵
H
,
使
得
Hx
< br>?
(
*
,
0,
0)
T
(
*
表示某非零元素
).
解
?
?
p>
sgn(
x
1
)<
/p>
x
2
?
(
?
1)
(
?
1)
2
?
2
2
?
(
?
2)
2
?
?
< br>3
;
u
?
x
?
?
e
1
?
(
?
p>
1
,
2,
?
2)
T
?
(3,
p>
0,
0)
T
?
p>
(
?
4,
2,
p>
?
2)
T
,
p>
其中
e
1
?
(1
,
0,
0)
p>
T
;
?
?
1
u
2
?
?
(
?
< br>?
x
1
)
?
?
3[
?
3
?
(
?
1)
]
?
12
,
2
则所求镜面反射矩阵为
2
?
1
0
0
?
?
?
4
?
?
?
1
3
2
3
?<
/p>
2
3
?
?
?
1
?
2
?
(
?
4,
2,
?
2)
?
?
2
3
2
< br>3
1
3
?
,
H
?
I
?
?
?
1
p>
uu
T
?
?
0
1
0
?
?
12
?
?
?
?
?
?
< br>?
?
0
0
1
?
?
?
?
2
?
?
?
p>
?
2
3
1
3
2
3
?
?
且
?
< br>?
1
3
2
3
?
2
3
?
?
?
1
?
p>
?
3
?
?
?
2
?
?
?
0
?
.
Hx
?
?
2
< br>3
2
3
1
3
?
?
?
?
?
?
?
?
p>
?
2
3
1
3
2
3
?
?
?
?
?
< br>2
?
?
?
?
0
?
?
可
以证明:
13
定理
9.3.6
对任意
n
阶
方阵
A
?
(
a
ij
)
n
?<
/p>
n
,存在正交矩阵
Q
,使得
B
?
Q
T
AQ
为形如式
(9.3.1)
的上
Hes
senberg
阵
.
证
记
?
a
p>
11
?
a
A
?
?
21
?
?
?
?
a
n
1
(1)
a
12
?
a
1
n
?
?
a
11
?
(1)
a
< br>22
?
a
2
n
?
?
?
?
a
21
?
?
?
?
?
?
p>
?
(1)
a
n
p>
2
?
a
nn
?
?
?
a
n
1
(1)
(1)
?
?
?
a
12
?
a
1
(1)
a
n
21
?
(1)
?
(1)
(1)
?
a
22
p>
?
a
2
a
n
?
?
A
1
,
x
1
< br>?
?
31
?
.
?
?
?
?
?
?
?
(
1)
?
(1)
(1)
< br>?
a
n
2
?
a
nn
?
?
?
?
a
n<
/p>
1
?
?
由式
p>
(9.3.4)
可构造
n
< br>?
1
阶对称正交矩阵
H
1
:
1
2
?
?
n
< br>2
?
?
?
1
?
sgn(
a
21
)
x
1
2
?
?
sgn(
a
21
)
?
?
a
i
1
?
,
?
i
?
p>
2
?
?
?
u
?
x
?
?
e
,
其中
e
?
(1,0,
?
,0)
T
?
R
n
?
1
,
(9.3.5)
1
1
1
1
?
1
?
1
u
2
?
?
(<
/p>
?
?
a
),
p>
?
=
1
1
21
?
1
2
1
2
?
H
?
I
?
?
?
1
u
u
T
,
?
1
1<
/p>
1
1
使得
H
p>
1
x
1
?
?
1
e
1
.
记
Q
1
?
?
?
I
1
?
?
n
?
n
,
且
,<
/p>
I
1
是
1
阶单位矩阵
.
显然
Q
1
是对称正交矩阵
.
用
Q
1
Q
< br>?
R
1
?
H
1
?
(1)
?
a
11
?
?
?
1
Q
1<
/p>
AQ
1
?
1
p>
?
Q
1
A
1
Q
1
?
?
0
?
?
< br>?
?
0
?
(2)
?
a
12
?
a
1
(2)
n
(2)
(2)
?
a
22
?
a
< br>2
n
?
(2)
< br>(2)
?
a
32
?
a
3
n
?
?
?
?
(2)
(2)
?
a
n
?
a
nn
2
?
对
A
作
相似变换,得
记
A
2
.
(9.3.6)
(2)
(2)
(2)
T
n
?
2
记
x
2
?
(
a
32
,同理可构造
n
?
2
阶对称正交矩阵
H
2
,使得
,
a
42
,
?
,
a
n
2
)
?
p>
R
H
2
x
2
?
?
2
e
1
(
< br>e
1
?
(1
,
0,
?
,
0)
T
?
R
n
?
2
).
记
Q
2
?
p>
?
?
I
2
?
?
,其中
I
2
为
2
阶单位矩阵,则
Q
2
仍是对称正交矩阵,用
Q
2
对
A
2
H
2
?
?
作相似变换,得
14
(1)
?
a
11
?
?
?
1
?
0<
/p>
?
1
Q
2
A
2
Q
2
?
Q
2
A
2
Q
2
?
?
?
0
?
?
?
?
0
?<
/p>
(2)
a
12
(
2)
a
22
?
2
0
?
0
(3
)
?
a
13
?
a
1
(3)
n
(3)
(3)
?
a
23
?
a
2
n
?
(3)
(3)
?
a
33
?
a
3
n
(
3)
(3)
?
a
43
?
a
4
n
?
?
?
?<
/p>
?
(3)
(3)
a
n
3
?
a<
/p>
nn
?
?
记
p>
A
3
.
(9.3.7)
依此类推,经过
k
步对称正交相似变换,得
1
Q
k
?
1
A
k
?
1
Q
k
?
?
1
?
Q<
/p>
k
?
1
A
k
?
1
Q
k
?
1
(1)
?
a
11
?
?
?
1
< br>?
0
?
?
0
?
?
0
?
?
0
?
0
p>
?
?
?
?
?
0
(2)
a
12
(2)
a
22
p>
(3)
a
13
?<
/p>
(3)
a
23
?
(3)
a
33
?
(
k
?
1)
a
1,
k
?<
/p>
1
(
k
?
1)
a
2,
k
?
1
(
k
?
1)
a
3,
k
?
1
a
1
(
k
k
)
(
k
)
a
2
k
(
k<
/p>
)
a
3
k
(
k
)
a
1,
k
?
1
(
k
)
a
< br>2,
k
?
1
(
k
)
a
3,
k
?
1
?
?
?
?
2
p>
0
0
0
0
?
0
?
3
0
0
0
?
< br>0
?
?
?
?
?
?
(
k
)
a
k
?
p>
1,
k
(
k
)
a
kk
(
k
)
a
k
?
1,
k
?
(
k
)
a
k
?
?
1,
k
?
1
(
k
)
a
k
,
p>
k
?
1
(
k
?
1)
?
a
k
?
1,
k
?
1
?
< br>k
?
1
0
?
0
?
(
k
)
a
k
?
p>
?
1,
k
?
1
?
(
k
)
a
nk
?
(
k
)
a
< br>n
,
k
?
1
?
k
)
?
a
1
(
n
p>
(
k
)
?
a
2
n
?
(
k
)
?
< br>a
3
n
?
?
?
(
k
)
?
a
k
?
p>
1,
n
?
(
k
)
a
kn
?
(
k
)
?
a
k
?
< br>1,
n
?
?
?
(
k
)
?
a
nn
?
记
A
k
.
(9.3.8)
重复上述过程,则有
(1)
?
a
1
1
?
?
?
1<
/p>
?
?
?
?
?
?
?
?
(2)
a
12
(2)
p>
a
22
(3)
a<
/p>
13
?
(3)
a
23
?
(3)
a
33
?
(
n
?
1)
a
1,
n
?
1
(
p>
n
?
1)
a
2,
n
?
1
(
n
?
1)
a
3,
n
?
1
n
)
?
< br>a
1
(
n
(
n
)
?
a
2
n
?
(
p>
n
)
?
a
3
n
?
?
?
(
n
)
< br>?
a
n
?
1,
n
?
(
n
)
a
nn
?
?
Q
n
?
p>
2
A
n
?
2
Q
?
1
n
?
2
?
< br>Q
n
?
2
A
n
?
2
Q
n
?
2
?
p>
2
记
?
3
?
?
(
n
?
1)
?
a
n
?
1,
n
< br>?
1
A
n
?
1
.
(9.3.9)
?
n
?
1
由式
(9.3.6)-(
9.3.9)
,可得
A
n
?
1
?
< br>Q
n
?
2
A
n
?
2
Q
n
?
2
?
p>
Q
n
?
2
Q
n
?
3
A
n
?
3
< br>Q
n
?
3
Q
n
?
2
?
Q
n
?
2
p>
Q
n
?
3
?
Q
1
AQ
1
?
Q
n
?
3
Q
n
?
2
.
若记
B
?
A
n
?
1
,
Q
?<
/p>
Q
1
Q
2
?
Q
n
?
2
,则
Q
为正交矩阵,且
有
B
?
Q
AQ
.
证毕!
注
1
p>
由定理
9.3.6
知:
因为任意
n
阶方阵
A
与
n
阶上
Hessenbe
rg
矩阵
B
相似,
所以求
T
A
的阵特征值的问题,就
可转化为求上
Hessenberg
矩阵
B
的特征值的问题
.
注
2
p>
若
A
是对称矩阵,则
B
也是对称矩阵
.
再由
B
的形式
(9.3.1)
知,此时
B
一定是
对称三对角阵
.
于是,求对称矩阵
A
的阵特征值的问题,便可转化为求对称三对角阵
B
的
特征值问题
.
例
9.3.2
设矩阵
15
?
1
2
1
?
2
p>
2
?
1
A
?
?
?
1
?
1
1
?
< br>?
2
1
1
解
第
1
步
记
A
1
p>
镜面反射阵
H
1
:
2
?
1
p>
?
?
1
?
?
1
?
试用镜面反射变换求正交矩阵
Q
,使
Q
T
AQ
为上
Hessenberg
矩阵
.
(1)
(1)
(1)
T
?
A
,
x
< br>1
?
(
a
21
,
a
31
,
a
41
)
?
(2,
1
,
2)
T
,利用式
(9.3.4)
构造三阶
?
1
?
p>
sgn(2)
x
1
2
?
2
2
?<
/p>
1
2
?
2
2
?
3
;
u
1
?
x
1
?
?
1
e
1
?
(2,
1
,
2)
T
?
(3,
0,
0)
T
?
(5,
< br>1
,
2)
T
,
其中
e
1
?
(1
,
0,
0)
T
;
(1)
?
1
?
1
u
?
?
(
?
?
a
)
p>
?
3(3
?
2)<
/p>
?
15
;
p>
1
1
1
21
2
2
2
则所求镜面反
射矩阵为
?
1
0
0
?
?
5
?
?
?
0.6
667
?
0.3333
?
0.6667
?
?
?
1
?
1
?
(5,
1,
2)
?
?
?
0.3333
0
.9333
?
0.1333
?
,
H
1
?
I
?
?
< br>1
?
1
u
1
u
1
T
?
?
0
1
0
p>
?
?
15
?
?
?
?
?
?
?
?
0
0
1
?
?
?
2
?
?
?
?
0.6667
?
< br>0.1333
0.7333
?
?
?
I
Q
1
p>
?
?
1
?
0
0
0
?
?
1
?
?
< br>?
?
0
?
0.6667
?
0.3333
?
p>
0.6667
?
?
,
?
H<
/p>
1
?
?
0
?
0.3333
0.9333
?
0.1333
?
?
?
0
?
0.6667<
/p>
?
0.1333
0.7333
?
?
?
3
0
0
?
?
1
?
?
3
2.3333
0.4667
?
0.06
67
?
?
.
A
2
?
Q
1<
/p>
T
A
1
Q
1
?
Q
1
AQ
1
?
?
?
0
0.4667
1.57
33
1.3467
?
?
?
0
?
0.0667
1.3467
0.0933
?
< br>?
(2)
(2)
T
第
2
步
记
x
2
p>
?
(
a
32
,
a
42
)
?
(0.4667,
?
0.0667)
T
,
利用式
(9.3.4)
构造
2
阶镜面反
射阵
H
2
:
?
2
?
sgn(0.4667)
x
2
2
?
0.4667
2
?
(
?
p>
0.0667)
2
?
0.4714
;
u
2
?
x
2
?
?
2
e
1
?
(0.4667,
?
0.0667)
T
?
(0.
4714,
0)
T
?
< br>(0.9381
,
?
0.066
7)
T
,
其
中
e
1
?
(1
,
0)
;
<
/p>
T
?
2
?
1
u
2
2
2
2
(2)
?
?
2
(
?
2
?
a
32
< br>)
?
0.4714(0.4714
?
0.4667)
?
0.4422<
/p>
;
则所求镜面反射矩阵为
16
?
1
0
?
?
?
p>
0.9901
0.1415
?
1
?
0.9381
?
?
1
T
H
2
?
I
?
?
2
u
2
u
2
?
?
?
(0.9381,
?
0.0667)
?
,
?
?
p>
?
?
?
?
0
1
?
0.4422<
/p>
?
?
0.0667
?
?
0.1415
0.9899
p>
?
?
I
Q
2
?
?
2
?
?
1
?
< br>?
?
0
?
H
2
?
?
?
0
?
?
0
p>
0
0
0
?
1
0
0
?
?
,
0
< br>?
0.9901
0.1415
?
?
0
0.1415
0.9899
?
?
3
0
0
?
?
< br>1
?
?
3
2.3333
?
0.4714
?
p>
0
T
?
.
A
3
?
Q
2
A
2
Q
2
?
Q
2
A
2
Q
2
?
?
?
0
?<
/p>
0.4714
1.5733
?
1.5000
?
?
?
p>
0
0
?
1.500
0
0.5000
?
?
< br>由此得正交矩阵
0
0
0
?
?
1
?
0
?
0.6667
0.2357
?
0.7071
?
?
,
Q
?
Q
1
Q
2
?
?
?<
/p>
0
?
0.3333
?
0.9429
0.0001
?
p>
?
?
0
?
0.6667
0.2357
0.7070
?
?
使得
?
3
0
0
?
?
1
?
< br>?
3
2.3333
?
0.4714
?
0
T
p>
?
Q
AQ
?
A
3
?
?
?
0
?
0.4714
1.1667
?
1.5000
?
?
< br>?
0
0
?
1.5000
0.5000
?
?
p>
为上
Hessenberg
矩阵。
9.3.2
QR
算法
Q
R
算法的基本思想是:利用
QR
分解得
到一系列与
A
相似的矩阵
?
A
k
?
,在一定的条
p>
件下,当
k
?
?<
/p>
时,
?
A
k
p>
?
收敛到一个以
A
的特征值
?
i
(
i
?
1,2,
?
,
n
)
为主对角线元素的上
三
角
矩
阵
.
首
先
介
绍
用
QR
分
解
(QR
decomposition
或
QR
factorization)
;
即
用
Householder
变换将矩阵
A
分解成正交矩阵
Q
< br>与上三角矩阵
R
的乘积,即
A<
/p>
?
QR
.
9.3.2.1
QR
分解
算法
9.3.2
(
QR
分解
)
(1)
(1)
(1)
T
第
1
步
记
A
的第
1
列为
x
1
?
(
a
p>
11
,
a
21
,
?
,
a
n
1
)
,
A
?
A
< br>1
?
(
a
ij
)
n
?
n
.
利用式
(9.3.4)
:
1
2
n
?
2
?
?
< br>(1)
(1)
?
?
1
?
sgn(
a
11
)
?
?
?
a
i
1
< br>?
?
,
?
i
?
1
?
?
?
T
n
?
p>
u
1
?
x
1
?
?
1
e
1
,
其中
e
1
?
(1,0,
?
,0)
?
R
,
?
(1)
?
=
?
(
?
?
a
?
1
1
1
11
),
?
H
?
I
?
?
?
1<
/p>
u
u
T
,
?
1
1
1
1
(1)
17
构造出的
H
1
是
n
阶对称正交矩阵,使得
H
1
x
1
?
?
1
e
< br>1
,从而有
?
?
1
?
0
A
2
?
H
1
A
1
?
?<
/p>
?
?
?
?
0
?
(2)
?
a
12
?
a
1
(2)
n
(2)
p>
(2)
?
a
22<
/p>
?
a
2
n
?
.
?
?
?
(2)
(2)
?
a
n
?
a
?
2
nn
?
(2)
(2)
(2)
p>
T
?
,使得
p>
n
?
1
阶对称正交
矩阵
H
第
2
步
记
x
p>
2
?
(
a
22
,
a
32
,
?
,
a
n
2
)
,同理可构造出
p>
2
n
?
(2)
p>
?
?
H
2
x
2
?
?
2
e
2
,其中
?
2
?
sgn(
a
22
)
?
?
a
i
2
< br>2
?
,
e
2
?
(1
,
0,
?
,
0)
T
?
R
n
?<
/p>
1
.
?
i
p>
?
2
?
1
2
?
1
H
?
若记
2
?
?
?
?
,
它仍是对称正交矩阵,于是有
?
H
2
?
(2
)
(2)
?
?
1
a
12
?
a
13
?
a
1<
/p>
(2)
n
?
(3
)
(3)
?
0
?
a
?
a
2<
/p>
23
2
n
?
p>
?
(3)
(3)
?
A
3
?
H
p>
2
A
2
?
?
0
.
0
a
33
?
a
3
n
?
?
?
?
?
?
?
?
(3)
(3)
?
?
0
0
a
n
?
a<
/p>
nn
3
?
?
p>
如此继续下去,直到完成第
n
?
1
步后,得到上三角矩阵
?
?
1
?
?
A
n
?
H
p>
n
?
1
A
n
?
1
?
?
?
?
?
< br>?
于是有
(2)
a
12
(2)
a
13
?
(3)
a
23
?
2
?
?
a
1
(2)
n
(3)
?
?
a
2
n
?
?
?
?
.
(
n
)
?
?
a
n
?<
/p>
1,
n
?
?
p>
n
?
?
A
n
?
H
n
?
1
A
n
< br>?
1
?
H
n
?
1
H
n
?
2
A
n
p>
?
2
?
?
?
H
n
?
1
H
n
?
< br>2
?
H
1
A
1
?
H
n
?
1
H
n
p>
?
2
?
H
1
A
.
令
R
?
A
< br>n
,
Q
?
H
1
H
2
?
H
n
?
1
p>
,
其中
Q
是对称正
交矩阵,
则
R
?
QA
.
因为
Q
是对称正交
矩阵,所以得
A
?
QR
.
注
1
若
A
非奇异
,则上三角矩阵
R
也非奇异,从而
R<
/p>
的主对角线元素不为零
.
注
2
若要求
R
的
主对角线元素取正数,则
A
的
QR
p>
分解是唯一的
.
例
9.3.3
求矩阵
?
1
0
?
1
p>
?
?
A
?
?
2
1
4
?
?
?
< br>?
?
2
3
0
?
?
18
的
Q
R
分解
A
?
Q
R
,并使矩阵
R
的主对角线上的元素都
是正数。
解
对
A
运用算
法
9.3.2.
第
1
步
记
A
1
p>
?
A
,
x
1
?
(1
,
2,
?
2)
T
,则
?
1
?
sign(1)
1
2<
/p>
?
2
2
?
(
?
2)
2
?
3
,
?
1
=
?
1
(
?
1
?
a
11
)
?
3(3
?
1)
?
12
,
u
1
?
x
1
?
?
1
e
1<
/p>
?
(1
,
2,<
/p>
?
2)
T
?
p>
(3,
0,
0)
T
?
(4,
2,
?
2)
T
,
e
1
?
(1
,<
/p>
0,
0)
T
,
?
1
0
p>
0
?
?
4
?
?
?
1
?
2
2
?
< br>?
?
1
?
2
?
(4,
2,
?
2)
?
1
?
?
2
2
1
?
,
H
1<
/p>
?
I
?
?
1
?
1
u
1
u
1
T
?
?
0
1
0
?
?
12
?
?
?
3
?
?
?
?
?
p>
0
0
1
?
?
?
?
2
?
?
?
2
< br>1
2
?
?
?
?
3
4
3
?
7
3
?
p>
?
.
A
2
?
H
1
A
1
?
?
0
5
3
10
3
< br>?
?
?
?
0
7
3
2
3
?
?
2
2
p>
第
2
步
记
x
2
p>
?
(5
3,
7
p>
3)
T
,
?
2
?
sign(5
3
)
(5
3)
?
(
7
3)
?
2
.86744
,
(2)
?
2
=
?
2
(
?
2
?
a
22
)
?
2.86744(2.86744
?
5
3)
?
13.0013
,
u
2
?
x
2
< br>?
?
2
e
1
?
(5
3,
7
3)
T
?
(2.86744,
0)
T
?
(4.53411
,
2.33333)
T
,
?
4.53411
?
1
?<
/p>
?
I
?
?
?
1
u
u
T
?
?
1
0
?
?
H
2
2
2
2
?
0
1
?
13
.0013
?
2.33333
?
(4.53411,
2.33333)
?
?
?
p>
?
?
?
0.581
24
?
0.81373
?
?
?
?
,
?
0.81373
0.58124
?
?
记
0
0
?
1
?<
/p>
1
0
?
?
?
?
,
H
2
?
?
?
0
?
0.58124
?
0.81373
?
?
?
?
?
0<
/p>
H
2
?
?
0
?
0.81373
0.58124
?
?
?
则
?
?
3
1.33333
?
2.3333
3
?
?
.
A
3
?
H
2
p>
A
2
?
?
0
?
2.86744
?
2.47995
?
?
< br>?
0
?
2.32494
?
?
0
?
?
?
1
0
< br>0
?
?
?
为了使
R
的主对角线上的元素都是正数,取
< br>H
3
?
0
?
1
0
,显然
H
3
是正交阵,
?
< br>?
?
?
0
0
?
1
?
?
且
19
-
-
-
-
-
-
-
-
-
上一篇:公司中文名称
下一篇:电磁流量计ModBus通讯协议