浅谈矩阵的LU 分解和QR 分解及其应用
基于理论研究和计算的需要,往往有必要把矩阵分解为具有某种特性的矩阵之积,这就是我们所说的矩阵分解.
本文将介绍两种常用的矩阵分解方法,以及其在解线性方程组及求矩阵特征值中的应用.
1. 矩阵的LU 分解及其在解线性方程组中的应用 1.1 高斯消元法
通过学习,我们了解到利用Gauss 消去法及其一些变形是解决低阶稠密矩阵方程组的有效方法. 并且近些年来利用此类方法求具有较大型稀疏矩阵也取得了较大进展. 下面我们就通过介绍Gauss 消去法,从而引出矩阵的LU 分解及讨论其对解线性方程组的优越性. 首先通过一个例子引入:
(1.1)
例1, 解方程组
(1.2) (1.3)
解. Step 1(1.1)⨯(-2) +(1.3) 消去(1.3)中未知数, 得到
-4x 2-x 3=-11(1.4)
Shep 2. (1.2)+(1.4) 消去(1.4)中的未知数x 2
⎧x 1+x 2+x 3=6⎛1⎫⎪ ⎪*
有⎨4x 2-x 3=5 显然方程组的解为x = 2⎪ 上述过程相当于 ⎪-2x =-6 3⎪
3⎩⎝⎭
6⎫⎛1116⎫⎛1116⎫⎛111
⎪ ⎪ ⎪
04-1504-1504-15~~⎪ ⎪ ⎪ 2-211⎪ 0-4-1-1⎪ 00-2-6⎪
⎭⎝⎝⎭⎝⎭
(-2)+ (r i 表示矩阵的i 行)
由此看出,消去法的基本思想是:用逐次消去未知数的方法把原方程化为与其等价的三角方程组.
下面介绍解一般n 阶线性方程组的Gauss 消去法.
⎛a 11 a 1n ⎫⎛x 1⎫⎛b 1⎫
⎪ ⎪ ⎪
X = ⎪b = ⎪ 则n 阶线性方程组 设A = ⎪
a ⎪ x ⎪ b ⎪⎝n1 a nn ⎭⎝n ⎭⎝n ⎭
AX =b (1.5)
并且A 为非奇异矩阵.
通过归纳法可以将AX =b 化为与其等价的三角形方程,事实上: 及方程(1.5)为A ()X =b (), 其中
1
1
A (1)=A b ()=b
1
(1)
(1) 设a 11
a i (1)
≠0,首先对行计算乘数m i1=1. 用-m i 1乘(1.5)的第一个方程加到第
m 11
1
i (i =2,3, ⋯, n )个方程上. 消去方程(1.5)的第2个方程直到第n 个方程的未知数x 1.
(1)(1)⎫(1)⎛a 11⋯a 1n ⎛x 1⎫⎛b 1⎫ ⎪ ⎪ ⎪
得到与(1.5)等价的方程组 ⎪ ⎪= ⎪简记作
(2)⎪ x ⎪ (2)⎪0⋯a nn ⎭⎝n ⎭⎝⎝b n ⎭
A (2)=b (2)(1.6)
(2)(1)(1)
其中a ij =a ij -m i 1a ij b i (2)=b i (1)-m i 1b 1(1)
(2) 一般第k (1≤k ≤n -1)次消去,设第k -1步计算完成. 即等价于
A (k )X =b (k ) (1.7)
(1)⎛a 11
= ⎝
()a 12
1
1
a 1(n )⎫
⎪(2)
a 2n ⎪
⎪
⎪
(k )(k )⎪a kk a kn
⎪
⎪(k )(k )⎪a nk a nna ⎭
(2)
a 22
且消去未知数x 1, x 2, ⋯, x k -1. 其中A (k )
(k )
设a kk ≠0计算m ik =a
(k )
ik
()a ik
/a (i k =+1, ⋯n , , ) 用m ik =k (i =k +1, ⋯, n ) 消去第
a kk
(k )kk
n
k +1个方程直到第n 个方程的未知数x k . 得到与(1. 7等) 价的方程组
最后可以把原方程化成一个与原方程等价的三A (k +1)X =b (k +)1故由数学归纳法知,
(i )角方程组. 但是以上分析明显存在一个问题,即使A 非奇异也无法保证a ii ≠0,
需要把非奇异的条件加强.
(i )
引理1 约化主元素a ii ≠(0i =1, ⋯, k ) 的充要条件是矩阵A 的顺序主子式
D i ≠0. 即
a 11 a ik
D 1=a 11≠0, D k = ≠0
a k 1⋯a kk
证明 利用数学归纳法证明引理的充分性. 显然,当k =1 时引理的充分性是成立的,现在假设引理对k -1是成立的,求证引理对k 亦成立. 有归纳法,设
(i )a ii ≠0(i =1,2⋯k -1)于是可用Gauss 消去法将中,即
(1)⎛a 11
= ⎝
()a 12
1
(2)a 22
A (1)→A (k )
(1)⎫ a 1n
⎪(2)
a 2n ⎪
⎪
⎪
(k )(k )⎪a kk a kn
⎪
⎪(k )(k )⎪a nk a nn ⎭
即
D 2=
(1)
a 11(1)a 12(2)a 22()a 12
1
1
(1)(2)(1)(2)(3)=a 11a 22 D 3=a 11a 22a 33
()
a 11
a 1(k )
1
D k =
()a 22
2(2)a 2k
()a kk
k
(1)(2)(k )=a 11a 22⋯a kk (1.8)
()
由设D i ≠0(i =1, ⋯, k ) 及式(1.8)有a kk ≠0
k
()
显然,由假设a ii ≠0(i =1,2⋯k ),利用(1.8)亦可以推出D i ≠0(i =1, ⋯, k )
i
从而由此前的分析易得;
定理1 如果n 阶矩阵A 的所有顺序主子式均不为零,则可通过Gauss 消去法(不进行交换两行的初等变换),将方程组(1.5)约化成上三角方程组,即
(1)⎛a 11
⎝
()a 12
1
()a 22
2
1(1)⎫a 1(n )⎫⎛x 1⎫⎛b 1⎪ (2)⎪(2) x ⎪a 2n ⎪ 2⎪ b 2⎪⎪ ⎪= ⎪ (1.9)
⎪ ⎪ ⎪
(n )⎪x b (n )⎪n ⎭⎝a nn ⎝n ⎭⎭
1.2 矩阵LU 分解
从而由以上讨论即能引出矩阵的LU 分解,通过高等代数我们得知对A 施行
(1)(2)
行初等变换相当于用初等矩阵左乘A ,即L 1A (1)=A (2) Lb 其中 =b 1
⎛1⎫
⎪-m 121⎪ L 1= ⎪ ⎪-m 0 1⎝n1⎭
一般第k 步消元,,相当于
L k A ()=A (
k
k +1)
L k b ()=b (
k
k +1)
重复这一过程,最后得到
(1)(n )⎧⎪L n -1⋯L 2L 1A =A
(1.10) ⎨(1)(n )
⎪⎩L n -1⋯L 2L 1b =b
其中
⎛ L k =
⎝
n
1
1-m k +1, k
-m n k
⎫⎪⎪⎪⎪
1⎪ ⎪
⎪1⎪⎭
-1-11
将上三角形矩阵A ()记作U , 由式(1.9)得到A =L 1L 2⋯L -n -1U =LU , 其中
⎛1⎫ ⎪m 1-1-11 21⎪ L =L 1L 2⋯L -n -1= ⎪ ⎪m m 1n 2⎝n 1⎭
由以上分析得;
定理2 (LU 分解) 设A 为n 阶矩阵,如果A 的顺序主子式
D i ≠0(i =1,2,, n -1) . 则A 可分解为一个单位下三角矩阵L 和一个上三角矩阵U 的乘积,且这种分解是唯一的.
证明 由先前的分析得出存在性是显然的,即A =LU . 下证唯一性, 设
A =LU =CD 其中L , C 为单位下三角矩阵,U , D 为上三角矩阵. 由于D -1L -1C =U D -1上式右端为上三角矩阵,左端为单位下三角矩阵,从而上式两端都
必须等于单位矩阵,故U =D , L =C . 证毕.
⎛111⎫ ⎪
例2 对于例子1 系数矩阵矩阵A = 04-1⎪由Gauss 消去法,得
2-21⎪⎝⎭
结合例1,故
⎛100⎫⎛111⎫
⎪⎪
A =LU = 010⎪04-1⎪
2-11⎪00-2⎪⎝⎭⎝⎭
对于一般的非奇异矩阵,我们可以利用初等排列矩阵I ki k (由交换单位矩阵I 的第k 行与第i k 行得到) ,即
(1)(2)(1)(2)⎧L I A =A , L I b =b 11i 1⎪11i 1
(1.11) ⎨(k )(k +1)(k )(k +1)
⎪⎩L k I ki k A =A , L k I k i k b =b
n
利用(1.11)得L n -1I n -1, i n -1L 1I 1i 1A =A ()=U . 简记做. 其中
下面就n 情况来考察一下矩阵.
A =A (5)=L 4I 4i 4L 3I 3i 3L 2I 2i 2L 1I 1i 1A =L 4I 4i 4L 3I 4i 4(I 4i 4I 3i 3L 2I 4i 4I 3i 3) ⨯(I 4i 4I 3i 3I 2i 2L 1I 4i 4I 3i 3I 2i 2) ( I 4i 4I 3i 3I 2i 2I 1i 1) A
()
从而记
从而
容易的为单位下三角矩阵, 总结以上讨论可得如下定理.
定理3 如果A 非奇异矩阵,则存在排列矩阵P 使PA =LU 其中L 为单位下三角矩阵,U 为上三角矩阵.
1.3 矩阵LU 分解的应用
以上对非奇异矩阵A 的LU 分解进行了全面的讨论,一下我们就简单介绍一下应用.
对于矩阵A 一旦实现了LU 分解,则解线性方程的问题, 便可以等价于:
(1)Ly =b 求y (2)Ux =y , 求x (1.12)
即,设A 为非奇异矩阵,且有分解式A =LU ,其中L 为单位下三角矩阵,U 为上三角矩阵。即
⎛1⎫⎛u 11u 12 ⎪ l 1u 2221⎪ A = ⎪ ⎪ l l 1⎝n 1n 2⎭⎝
u 1n ⎫
⎪
u 2n ⎪
(1.13)
⎪
⎪
u nn ⎭
下面说明L , 的元素可以由n 步直接计算出: 由(1.13)有
a 1i =u 1i (1.14)
再由矩阵乘法得a 1i =l i 1u 11, 故有
l i 1=a 1i /u 11
于是得U 的第一列元素.
设已经得U 的第一行到r -1行元素与L 的第一列到r -1列元素,由(1.13)有
a ri =∑l rk u ki =∑l rk u ki +u ri
k =1
k =1
n
r -1
故有
u ri =a ri -∑l rk u ki (i ) (1.15)
k =1r -1
再由(1.13)得
a ir =∑l ik u kr =∑l ik u kr +l ir u rr .
k =1
k =1
n
r -1
得
l ir =(a ir -∑l ik u k r ) /u r r (1.16)
k =1r -1
故有以上分析结合(1.12)得;
y 1=b 1⎧⎪i -1 ⎨
⎪y i =b i -∑l ik y k
k =1⎩
x n =y n /u nn ⎧
⎪n
(1.18) ⎛⎫⎨
x =b -u x /u ⎪i i ∑ik k ⎪i i
k =i +1⎝⎭⎩
⎛123⎫⎛x 1⎫⎛14⎫
⎪⎪ ⎪
例3. 求解 252⎪x 2⎪= 18⎪
315⎪x ⎪ 20⎪⎝⎭⎝3⎭⎝⎭
解. 由 (1.15),(1.16)计算,得
⎛100⎫⎛123⎫
⎪⎪
A = 210⎪01-4⎪=LU
3-51⎪00-24⎪⎝⎭⎝⎭
⎛14⎫⎛14⎫⎛14⎫⎛1⎫
⎪ ⎪ ⎪ ⎪
求解Ly = 18⎪ 得y = -10⎪ 求解Ux = -10⎪ 得x = 2⎪
3⎪ 20⎪ -72⎪ -72⎪
⎝⎭⎝⎭⎝⎭⎝⎭
以上便是通过介绍Gauss 消去法,引出矩阵的LU 分解,这种分解在数值分析中,在设计算法求解高维线性方程组能提高效率,不像Grammar 法则只是从理论上解决了非奇异线性方程组的解法,实际操作起来是十分不方便的,而利用
LU 分解的基本原理能大大减少计算过程的繁琐.
2.矩阵的QR 分解及其在计算矩阵特征值的应用 2.1 转化非奇异矩阵为上Hessebberg
定义1 一方阵,如果i >j +1时有b ij =0. 则称矩阵B 为上Hessebberg 阵,即
⎛b 11
b 21 0B =
0⎝
b 12 b 22b 320 0
b 1n ⎫
⎪
b 2n ⎪ ⎪
⎪
⎪ ⎪
⎪
b n , n -1b nn ⎪⎭
定义2 设向量w 满足w 2, 矩阵H =I -2ww T 称为初等反射矩阵,记作
H (w ), 即
⎛1-2w 12-2w 1w 2⋯-2w 1w n ⎫
⎪2
⋯-2w w -2w w 1-2w 2n ⎪212
H (w )= 其中
⎪ ⎪2⎪ -2w w 1-2w n ⎭n 1⎝
定理2.1 初等反射阵H 是对称矩阵,正交矩阵和合同矩阵.
定理2.2 设x , y 为两个不相等的n 维向量,但x 2=y 2, 则存在一个初等反射矩阵H ,使Hx =y .
证明 令w =(x -y ) /x -y 2, 则得到一个初等反射矩阵
H =I -2w w T =I -2(x T -y T )(x -y )/x -y 2
2
而且
Hx =x -2(x -y )
x x -y
22
(x
T
-y T )=x -2(x -y )(x T x -y T x )/x -y 2
2
因为
x -y
22
=(x -y ) (x -y ) T =2(x T x -y T x )
所以Hx =x -(x -y ) =y . 并且由代数学知识易知,w =(x -y ) /x -y 2. 成立的唯一长度等于1的向量.
推论 设向量x ∈R n (x ≠0), σ=±x 2, 且x ≠-σe 1则存在一个初等反射矩阵
H =I -2uu T /u
22
=I -ρ-1uu T , 使Hx =-σe 1, 其中u =x +σe 1 ,ρ=u 2/2
2
设x =(α1, α2, α3, ⋯, αn )u =(u 1, u 2, ⋯, u n ), 则u =(α1+σ, α2, ⋯, αn ) ,
ρ=u 2/2=σ(σ+α1). 如果σ,那么计算σ+α1时有效数字可能损失,取有相同符号,即取
2
σ=sgn (α1)x 2
有了以上关于初等反射矩阵的性质接下来讨论正交相似约化一般矩阵和对角矩阵. 设
a 1n ⎫⎛a 11a 12
⎪⎛
⋯a a a a 2n ⎪2122
A = = 11
(1) ⎪ a
⎪⎝21⎝a n 1a n 2 a nn ⎭
()⎫A 12
⎪, (1)⎪A 22⎭
1
(1)(1)
不妨设a 21,其中 ≠0,否则这一步不需约化,选择初等矩阵R 1使R 1a 21
1n
⎧22⎪σ1=sgn (a 21)(∑a i 1) ,
i =2⎪
(1)⎪u 1=a 21+σ1e 1,
(2.1) ⎨
12⎪ρ=u ⎪1212=σ1(σ1+a 21), ⎪T
R 1=I -ρ1-1u 1u 1. ⎩
⎛10⎫
令U 1= ⎪则
0R ⎝1⎭
⎛a 11
A 2=U 1AU 11= R a (1)
⎝121
2
(1)(2)A 21R 1⎫⎛A 11
⎪= (1)⎪ 0R 1A 22R 1⎭⎝
(2)
a 12(2)a 22
()⎫A 13
⎪, (2)⎪A 23⎭
2
()
其中A 11∈R 2⨯1, 一直重复这一过程得;
定理
2.1 如果A =R n ⨯n ,则存在初等反射矩阵, 使
U n -2⋯U 2U 1AU 1U 2⋯U n -2=C (上Hessenberg ).
推论A 为对称矩阵,则存在初等反射矩阵 其中A n -1为
⎛c 1b 1⎫ ⎪b c b 122 ⎪ ⎪ A n -1= ⎪
⎪
b n -2c n -1b n -1⎪ ⎪ ⎪b c n -1n ⎭⎝
2.2 矩阵的QR 分解
有了以上的讨论我们知道任何一个方阵A 都可以正交相似化为Hessenberg 或更为特殊的对称三角矩阵,接下来便要使用矩阵的QR (Q 为正交矩阵,R 为三角矩阵)分解来解决Hessenberg 阵和对称角矩阵的全部特征值问题.
引理1 设, 其中αi ,αj 不全为零,则可选取一平面矩阵P ij 使
P ij x =y =(α1, α2, ⋯, αi (1), ⋯, α(j 1), ⋯, αn ) T , 其中
αi (
1)=(2.2)
(2.3)
⎧
c =⎪⎪
(2.4) ⎨
α⎪s =⎪⎩
证明 事实上P ij x 只改变x 的第i 个及第j 个元素,且有
αi (1)=c αi +s αj ,,
于是可选P 即按式(2.4)求c , s ,且有式(2.2)及(2.3). ij
定理2.4 如果A 为非奇异矩阵,则即存在平面旋转矩阵P 1, P 2, ⋯, P n -1(一系列平面旋转矩阵)使
且r ii >0(i =1,2, ⋯, n -1).
证明 由于A 第一列一定存在a j 1≠0,于是,如果a j 1≠0(j =2,3, ⋯, n ), 应用引理1,即存在平面旋转矩阵P 12, P 13, ⋯, P 1n , 使
2(2)⎛r 11a 12a 1(n )⎫
⎪(2)(2)
0a 22⋯a 2n ⎪(2)
P ⋯P P A =1n 1312 ⎪=A
⎪
0a (2) a (2)⎪
n 2nn ⎭⎝
且记P 1n ⋯P 13P 12=P 1,同理重复上述过程,最后得到:存在正交矩阵P 1, P 2, ⋯, P n -1使式(2.5)成立.
定理2.5(矩阵的QR 分解) 如果n 阶方阵A 为非奇异矩阵,则A 可分解为一正交矩阵Q 与一上三角矩阵R , 且当R 对角元素都为正数时分解唯一. 证明 由定理2.4知,存在正交矩阵P 1, P 2, ⋯, P n -1, 使
(2.6)
为上三角阵,记
Q T =P n -1P n -2⋯P 1
T T T
于是式(2.6)为Q T A , 即A =QR , 其中Q =P 为正交矩阵. P ⋯P 12n -1
现证明唯一性. 设有A =Q 1R 1=Q 2R 2,其中R 1,R 2为上三角矩阵(显然为非奇
异矩阵)且对角元素都为正. Q 1, Q 2为正交矩阵. 于是
T
Q 2Q 1=R 2R 1-1(2.7)
由式(2.7)知上三角阵R 2R 1-1为正交矩阵,故R 2R 1-1为对角矩阵,即
R 2R 1-1=D =diag (d 1, d 2, ⋯d n )
因为R 2R 1-1是正交矩阵,所以D 2=I , 又因为R 1,R 2对角元素都为正数,故
d i >0(i =1,2, ⋯, n ),即D =I 于是R 1=R 2,由式(2.6)得到Q 1=Q 2.
2.3 矩阵特征值QR 算法理论依据
设A =A 1=(a ij )∈R n ⨯n 非奇异, 且对A 直接QR 分解,即A =QR , 其中R 为上三角矩阵, Q 为正交矩阵,于是得到一个新矩阵 B =RQ =Q T AQ . 显然, B 是由A 经过正交相似变换得到, 因此A 和B 有相同的的特征值, 在对B 进行QR 分解,又可得到一新矩阵, 重复可得矩阵序列.
设A =A 1, 进行QR 分解, 得A 1=Q 1R 1, 作矩阵A 2=R 1Q 1=Q 1T A 1Q 1 ,,求得A K 后
T 将其进行QR 分解, 得A k =Q k R k , 作矩阵A k +1=R k Q k =Q k A k Q k 而利用QR 分解求
矩阵的特征值, 就是利用矩阵QR 分解, 按上述递推法则构造矩阵序列{A k }的过程, 只要A 非奇异, 利用上述方法就能构造出{A k }.
定理2.6(基本的QR 方法) 设A =A 1=(a ij )∈R n ⨯n 非奇异,QR 基本方法为:
T
⎧R k 为上三角矩阵), ⎪A k =Q k R k (Q k Q k =I
⎨
A k +1=R k Q (⎪k k =1,2, ⋯), ⎩
且记=Q 1⋯Q k ,则有;
T
a:A k +1相似与A k 即 A k +1=Q k A k Q k T b: A k +1=(Q 1⋯Q k )A 1(Q 1⋯Q k )
c: A k 的QR 分解式为 A k
T 证明 a:A k +1=R k Q k =Q k A k Q k
b:
T A k +1=Q k A k Q k
c: 显然,k =1 时有A 1,设A k -1有分解式,于是
=A .
结合定理2.5,可以得出A k +1可由A k 按下述方法求得.
T T T
(1)左变换P n -1P n -2⋯P 1A k =R k (上三角阵) (2)右变换R k P 1P 2⋯P n -1=A k +1
引理2 设M k , 其中Q k ,R k 为具有正对角元素的上三角矩阵,如果
M k →I (k →∞), 则Q k →I , 及R k →I .
(k )T T T
证明 设R k R k =M k M k →I (k →∞) , 记R k =(r ij ) , 矩阵R k R k
()()()
, r 11∙(r 11, r 12, ⋯, r 1(n ))
k
k
k
k
(k )
因此有 r 11,,, (2.8)
()()()()()
r 12∙(r 11, r 12, ⋯, r 1(n ))+r 22∙(0, r 22, ⋯, r 2(n ))
k
k
k
k
k
k
k
利用式(2.8)的结果,则有
(k )
,,, 0(k →∞) r 22
-1T
对于R k R k 其他行同理可得, 故R k , 且易知有R k →I (k →∞), 因此
Q k =M k R k -1→I (k →∞).
定理2.7 设A =(a ij )∈R n ⨯n 如果A 的特征值满足:;
A 有标准型A =XDX -1其中(), 且设X -1有三角分解(L 为单位下三角矩阵,U
为单位上三角矩阵),则有QR 算法产生{A k }本质上收敛于上三角阵,即
(k →∞)
这个定理证明过程比较繁琐,具体证明内容参考文献[2],这里需要说明在证明过
程中A k 的收敛速度依赖r n =λn /λn -1, 当r n 越小时, 收敛越快.
推论2.1 如果A 是对称矩阵满足定理2.7条件, 则由QR 方法产生{A k }本质上收敛于对角矩阵.
不管是上三角矩阵,还是对角矩阵,其特征值得计算是很方便的,注意到对角矩阵或上三角矩阵是极限形式,故在给定的误差范围内,计算满足定理2.7的矩阵特征值是比较容易的. 以上对于理论知识的解释完毕,接下来看一下如何具体在应用这一原理计算矩阵特征值.
2.3 矩阵特征值QR 算法实际操作
(k )
在定理2.7的证明中,我们发现a nn (k →∞)速度很大程度依赖于比值
当r n 很小时,收敛很快. 如果s 为r n 一个估计,且对A -sI 应用QR r n =λn /λn -1, ,
方法,收敛速度恰当,收敛速度会更快.
由于QR 计算矩阵特征值是一种利用极限近似的方法,故如果{A k }收敛的越快,计算将会更方便. 为此,为了加速收敛,选择数列{s k }, 按照下述方法构造矩阵序列{A k }.
Shep 1 设A =A 1=(a ij )∈R n ⨯n
Shep 2 将A k -s k I 进行QR 分解,即A k -s k I , k =1,2, ;
T
Shep 3 构造新矩阵A k +1=R k Q k +s k I =Q k A k Q k ;
Shep 4A k +1 ,其中 ;
Shep 5 矩阵(A -s 1I )(A -s 2I ) ⋯(A -s n I ) =φ(A )有QR 分解式
φ(A );
Shep 6 结合定理2.5的结论:首先利用正交变换(左变换)将A k -s k I 化为上
三角阵, 即
P n -1P n -2⋯P 1(A k -s k I )
T
其中Q k 为一系列平面旋转矩阵的乘积. 于是
A k +1(2.9)
有前面的讨论我们知道任何一个方阵都可以正交相似化为一个上Hessenberg 矩阵,而有正交过程中矩阵特征值不变性质,故下面我们只考虑Hessenberg 矩阵的特质值得计算.
设A =A 1为上Hessenberg g 结合(2.1)的具体P i 的构造方法shep 1~6, 构造出上Hessenberg 矩阵序列{A k },在构造的过程中对{s k }的选择也是十分重要的,
(k )一般选取s k =a nn . 具体方法我们看下面这个例题的解决.
例2.1 用QR 方法计算对称三角矩阵的全部特征值
⎛210⎫
⎪
A =A 1= 131⎪
014⎪⎝⎭
解 首先选取s 1 结合2.1的方法构造旋转平面
⎛2.2361-1.3420.4472⎫ ⎪
P 23P (A -s I )=R =01.0954-0.36511211 ⎪
0-0.36510.81650⎪⎝⎭0⎫⎛1.40000.4899 ⎪T T
A 2=RP 3.26670.7454⎪ 12P 23+s 1I = 0.4899
00.74544.3333⎪⎝⎭
重复以上的操作
0⎫⎛1.29150.2017
⎪A 3= 0.20173.02020.2724⎪
00.27244.6884⎪⎝⎭0⎫⎛1.27370.0993 ⎪A 4= 0.09932.99430.0072⎪
00.00724.7320⎪⎝⎭
0⎫⎛1.26940.0498
⎪A 5= 0.04982.99860⎪
004.7321⎪⎝⎭
记
现在收缩,继续对A 5的子矩阵进行变换,得到
-s 5I
故求得A 近似特征值为
1.2680
且A 的实际特征值是
=3.0 =
结束语
以上便是从矩阵的两个常见分解LU , QR 分解出发,进而把它们应用到解线性方程组,以及求矩阵的特征值,这种方法在解决阶数比较高的矩阵时利用这种思想设计出相应的的算法,往往能十分有效的减缓繁重的计算过程.
参考文献
[1]王萼芳,石生明. 高等代数[M].北京:高等教育出版社,1978.162-203 [2]李庆阳,王能超,易大义. 数值分析[M].武汉:华中科技大学出版社,1982.173-245
[3]周金土. 高等代数解题思想和方法[M].杭州:浙江大学出版社,2008.96-130
[4]徐仲,陆全,张凯院. 高等代数考研教案[M].西安:西北工业大学出版社,2005.155-198
[5]郑崇友,樊磊,崔宏斌.Frame 与连续格[M].北京:首都师范大学出版社,1994
[6]Gierz G ,Hoffmann K H,et al.A Compendium of Continuous Lattice[M].Berlin:Springer-Verlag,1980.
[7]吴耀强. 关于代数L-domain 的一个刻画定理. 江西师范大学学报(自然科学版)[J].2006,30(6):573-576.
浅谈矩阵的LU 分解和QR 分解及其应用
基于理论研究和计算的需要,往往有必要把矩阵分解为具有某种特性的矩阵之积,这就是我们所说的矩阵分解.
本文将介绍两种常用的矩阵分解方法,以及其在解线性方程组及求矩阵特征值中的应用.
1. 矩阵的LU 分解及其在解线性方程组中的应用 1.1 高斯消元法
通过学习,我们了解到利用Gauss 消去法及其一些变形是解决低阶稠密矩阵方程组的有效方法. 并且近些年来利用此类方法求具有较大型稀疏矩阵也取得了较大进展. 下面我们就通过介绍Gauss 消去法,从而引出矩阵的LU 分解及讨论其对解线性方程组的优越性. 首先通过一个例子引入:
(1.1)
例1, 解方程组
(1.2) (1.3)
解. Step 1(1.1)⨯(-2) +(1.3) 消去(1.3)中未知数, 得到
-4x 2-x 3=-11(1.4)
Shep 2. (1.2)+(1.4) 消去(1.4)中的未知数x 2
⎧x 1+x 2+x 3=6⎛1⎫⎪ ⎪*
有⎨4x 2-x 3=5 显然方程组的解为x = 2⎪ 上述过程相当于 ⎪-2x =-6 3⎪
3⎩⎝⎭
6⎫⎛1116⎫⎛1116⎫⎛111
⎪ ⎪ ⎪
04-1504-1504-15~~⎪ ⎪ ⎪ 2-211⎪ 0-4-1-1⎪ 00-2-6⎪
⎭⎝⎝⎭⎝⎭
(-2)+ (r i 表示矩阵的i 行)
由此看出,消去法的基本思想是:用逐次消去未知数的方法把原方程化为与其等价的三角方程组.
下面介绍解一般n 阶线性方程组的Gauss 消去法.
⎛a 11 a 1n ⎫⎛x 1⎫⎛b 1⎫
⎪ ⎪ ⎪
X = ⎪b = ⎪ 则n 阶线性方程组 设A = ⎪
a ⎪ x ⎪ b ⎪⎝n1 a nn ⎭⎝n ⎭⎝n ⎭
AX =b (1.5)
并且A 为非奇异矩阵.
通过归纳法可以将AX =b 化为与其等价的三角形方程,事实上: 及方程(1.5)为A ()X =b (), 其中
1
1
A (1)=A b ()=b
1
(1)
(1) 设a 11
a i (1)
≠0,首先对行计算乘数m i1=1. 用-m i 1乘(1.5)的第一个方程加到第
m 11
1
i (i =2,3, ⋯, n )个方程上. 消去方程(1.5)的第2个方程直到第n 个方程的未知数x 1.
(1)(1)⎫(1)⎛a 11⋯a 1n ⎛x 1⎫⎛b 1⎫ ⎪ ⎪ ⎪
得到与(1.5)等价的方程组 ⎪ ⎪= ⎪简记作
(2)⎪ x ⎪ (2)⎪0⋯a nn ⎭⎝n ⎭⎝⎝b n ⎭
A (2)=b (2)(1.6)
(2)(1)(1)
其中a ij =a ij -m i 1a ij b i (2)=b i (1)-m i 1b 1(1)
(2) 一般第k (1≤k ≤n -1)次消去,设第k -1步计算完成. 即等价于
A (k )X =b (k ) (1.7)
(1)⎛a 11
= ⎝
()a 12
1
1
a 1(n )⎫
⎪(2)
a 2n ⎪
⎪
⎪
(k )(k )⎪a kk a kn
⎪
⎪(k )(k )⎪a nk a nna ⎭
(2)
a 22
且消去未知数x 1, x 2, ⋯, x k -1. 其中A (k )
(k )
设a kk ≠0计算m ik =a
(k )
ik
()a ik
/a (i k =+1, ⋯n , , ) 用m ik =k (i =k +1, ⋯, n ) 消去第
a kk
(k )kk
n
k +1个方程直到第n 个方程的未知数x k . 得到与(1. 7等) 价的方程组
最后可以把原方程化成一个与原方程等价的三A (k +1)X =b (k +)1故由数学归纳法知,
(i )角方程组. 但是以上分析明显存在一个问题,即使A 非奇异也无法保证a ii ≠0,
需要把非奇异的条件加强.
(i )
引理1 约化主元素a ii ≠(0i =1, ⋯, k ) 的充要条件是矩阵A 的顺序主子式
D i ≠0. 即
a 11 a ik
D 1=a 11≠0, D k = ≠0
a k 1⋯a kk
证明 利用数学归纳法证明引理的充分性. 显然,当k =1 时引理的充分性是成立的,现在假设引理对k -1是成立的,求证引理对k 亦成立. 有归纳法,设
(i )a ii ≠0(i =1,2⋯k -1)于是可用Gauss 消去法将中,即
(1)⎛a 11
= ⎝
()a 12
1
(2)a 22
A (1)→A (k )
(1)⎫ a 1n
⎪(2)
a 2n ⎪
⎪
⎪
(k )(k )⎪a kk a kn
⎪
⎪(k )(k )⎪a nk a nn ⎭
即
D 2=
(1)
a 11(1)a 12(2)a 22()a 12
1
1
(1)(2)(1)(2)(3)=a 11a 22 D 3=a 11a 22a 33
()
a 11
a 1(k )
1
D k =
()a 22
2(2)a 2k
()a kk
k
(1)(2)(k )=a 11a 22⋯a kk (1.8)
()
由设D i ≠0(i =1, ⋯, k ) 及式(1.8)有a kk ≠0
k
()
显然,由假设a ii ≠0(i =1,2⋯k ),利用(1.8)亦可以推出D i ≠0(i =1, ⋯, k )
i
从而由此前的分析易得;
定理1 如果n 阶矩阵A 的所有顺序主子式均不为零,则可通过Gauss 消去法(不进行交换两行的初等变换),将方程组(1.5)约化成上三角方程组,即
(1)⎛a 11
⎝
()a 12
1
()a 22
2
1(1)⎫a 1(n )⎫⎛x 1⎫⎛b 1⎪ (2)⎪(2) x ⎪a 2n ⎪ 2⎪ b 2⎪⎪ ⎪= ⎪ (1.9)
⎪ ⎪ ⎪
(n )⎪x b (n )⎪n ⎭⎝a nn ⎝n ⎭⎭
1.2 矩阵LU 分解
从而由以上讨论即能引出矩阵的LU 分解,通过高等代数我们得知对A 施行
(1)(2)
行初等变换相当于用初等矩阵左乘A ,即L 1A (1)=A (2) Lb 其中 =b 1
⎛1⎫
⎪-m 121⎪ L 1= ⎪ ⎪-m 0 1⎝n1⎭
一般第k 步消元,,相当于
L k A ()=A (
k
k +1)
L k b ()=b (
k
k +1)
重复这一过程,最后得到
(1)(n )⎧⎪L n -1⋯L 2L 1A =A
(1.10) ⎨(1)(n )
⎪⎩L n -1⋯L 2L 1b =b
其中
⎛ L k =
⎝
n
1
1-m k +1, k
-m n k
⎫⎪⎪⎪⎪
1⎪ ⎪
⎪1⎪⎭
-1-11
将上三角形矩阵A ()记作U , 由式(1.9)得到A =L 1L 2⋯L -n -1U =LU , 其中
⎛1⎫ ⎪m 1-1-11 21⎪ L =L 1L 2⋯L -n -1= ⎪ ⎪m m 1n 2⎝n 1⎭
由以上分析得;
定理2 (LU 分解) 设A 为n 阶矩阵,如果A 的顺序主子式
D i ≠0(i =1,2,, n -1) . 则A 可分解为一个单位下三角矩阵L 和一个上三角矩阵U 的乘积,且这种分解是唯一的.
证明 由先前的分析得出存在性是显然的,即A =LU . 下证唯一性, 设
A =LU =CD 其中L , C 为单位下三角矩阵,U , D 为上三角矩阵. 由于D -1L -1C =U D -1上式右端为上三角矩阵,左端为单位下三角矩阵,从而上式两端都
必须等于单位矩阵,故U =D , L =C . 证毕.
⎛111⎫ ⎪
例2 对于例子1 系数矩阵矩阵A = 04-1⎪由Gauss 消去法,得
2-21⎪⎝⎭
结合例1,故
⎛100⎫⎛111⎫
⎪⎪
A =LU = 010⎪04-1⎪
2-11⎪00-2⎪⎝⎭⎝⎭
对于一般的非奇异矩阵,我们可以利用初等排列矩阵I ki k (由交换单位矩阵I 的第k 行与第i k 行得到) ,即
(1)(2)(1)(2)⎧L I A =A , L I b =b 11i 1⎪11i 1
(1.11) ⎨(k )(k +1)(k )(k +1)
⎪⎩L k I ki k A =A , L k I k i k b =b
n
利用(1.11)得L n -1I n -1, i n -1L 1I 1i 1A =A ()=U . 简记做. 其中
下面就n 情况来考察一下矩阵.
A =A (5)=L 4I 4i 4L 3I 3i 3L 2I 2i 2L 1I 1i 1A =L 4I 4i 4L 3I 4i 4(I 4i 4I 3i 3L 2I 4i 4I 3i 3) ⨯(I 4i 4I 3i 3I 2i 2L 1I 4i 4I 3i 3I 2i 2) ( I 4i 4I 3i 3I 2i 2I 1i 1) A
()
从而记
从而
容易的为单位下三角矩阵, 总结以上讨论可得如下定理.
定理3 如果A 非奇异矩阵,则存在排列矩阵P 使PA =LU 其中L 为单位下三角矩阵,U 为上三角矩阵.
1.3 矩阵LU 分解的应用
以上对非奇异矩阵A 的LU 分解进行了全面的讨论,一下我们就简单介绍一下应用.
对于矩阵A 一旦实现了LU 分解,则解线性方程的问题, 便可以等价于:
(1)Ly =b 求y (2)Ux =y , 求x (1.12)
即,设A 为非奇异矩阵,且有分解式A =LU ,其中L 为单位下三角矩阵,U 为上三角矩阵。即
⎛1⎫⎛u 11u 12 ⎪ l 1u 2221⎪ A = ⎪ ⎪ l l 1⎝n 1n 2⎭⎝
u 1n ⎫
⎪
u 2n ⎪
(1.13)
⎪
⎪
u nn ⎭
下面说明L , 的元素可以由n 步直接计算出: 由(1.13)有
a 1i =u 1i (1.14)
再由矩阵乘法得a 1i =l i 1u 11, 故有
l i 1=a 1i /u 11
于是得U 的第一列元素.
设已经得U 的第一行到r -1行元素与L 的第一列到r -1列元素,由(1.13)有
a ri =∑l rk u ki =∑l rk u ki +u ri
k =1
k =1
n
r -1
故有
u ri =a ri -∑l rk u ki (i ) (1.15)
k =1r -1
再由(1.13)得
a ir =∑l ik u kr =∑l ik u kr +l ir u rr .
k =1
k =1
n
r -1
得
l ir =(a ir -∑l ik u k r ) /u r r (1.16)
k =1r -1
故有以上分析结合(1.12)得;
y 1=b 1⎧⎪i -1 ⎨
⎪y i =b i -∑l ik y k
k =1⎩
x n =y n /u nn ⎧
⎪n
(1.18) ⎛⎫⎨
x =b -u x /u ⎪i i ∑ik k ⎪i i
k =i +1⎝⎭⎩
⎛123⎫⎛x 1⎫⎛14⎫
⎪⎪ ⎪
例3. 求解 252⎪x 2⎪= 18⎪
315⎪x ⎪ 20⎪⎝⎭⎝3⎭⎝⎭
解. 由 (1.15),(1.16)计算,得
⎛100⎫⎛123⎫
⎪⎪
A = 210⎪01-4⎪=LU
3-51⎪00-24⎪⎝⎭⎝⎭
⎛14⎫⎛14⎫⎛14⎫⎛1⎫
⎪ ⎪ ⎪ ⎪
求解Ly = 18⎪ 得y = -10⎪ 求解Ux = -10⎪ 得x = 2⎪
3⎪ 20⎪ -72⎪ -72⎪
⎝⎭⎝⎭⎝⎭⎝⎭
以上便是通过介绍Gauss 消去法,引出矩阵的LU 分解,这种分解在数值分析中,在设计算法求解高维线性方程组能提高效率,不像Grammar 法则只是从理论上解决了非奇异线性方程组的解法,实际操作起来是十分不方便的,而利用
LU 分解的基本原理能大大减少计算过程的繁琐.
2.矩阵的QR 分解及其在计算矩阵特征值的应用 2.1 转化非奇异矩阵为上Hessebberg
定义1 一方阵,如果i >j +1时有b ij =0. 则称矩阵B 为上Hessebberg 阵,即
⎛b 11
b 21 0B =
0⎝
b 12 b 22b 320 0
b 1n ⎫
⎪
b 2n ⎪ ⎪
⎪
⎪ ⎪
⎪
b n , n -1b nn ⎪⎭
定义2 设向量w 满足w 2, 矩阵H =I -2ww T 称为初等反射矩阵,记作
H (w ), 即
⎛1-2w 12-2w 1w 2⋯-2w 1w n ⎫
⎪2
⋯-2w w -2w w 1-2w 2n ⎪212
H (w )= 其中
⎪ ⎪2⎪ -2w w 1-2w n ⎭n 1⎝
定理2.1 初等反射阵H 是对称矩阵,正交矩阵和合同矩阵.
定理2.2 设x , y 为两个不相等的n 维向量,但x 2=y 2, 则存在一个初等反射矩阵H ,使Hx =y .
证明 令w =(x -y ) /x -y 2, 则得到一个初等反射矩阵
H =I -2w w T =I -2(x T -y T )(x -y )/x -y 2
2
而且
Hx =x -2(x -y )
x x -y
22
(x
T
-y T )=x -2(x -y )(x T x -y T x )/x -y 2
2
因为
x -y
22
=(x -y ) (x -y ) T =2(x T x -y T x )
所以Hx =x -(x -y ) =y . 并且由代数学知识易知,w =(x -y ) /x -y 2. 成立的唯一长度等于1的向量.
推论 设向量x ∈R n (x ≠0), σ=±x 2, 且x ≠-σe 1则存在一个初等反射矩阵
H =I -2uu T /u
22
=I -ρ-1uu T , 使Hx =-σe 1, 其中u =x +σe 1 ,ρ=u 2/2
2
设x =(α1, α2, α3, ⋯, αn )u =(u 1, u 2, ⋯, u n ), 则u =(α1+σ, α2, ⋯, αn ) ,
ρ=u 2/2=σ(σ+α1). 如果σ,那么计算σ+α1时有效数字可能损失,取有相同符号,即取
2
σ=sgn (α1)x 2
有了以上关于初等反射矩阵的性质接下来讨论正交相似约化一般矩阵和对角矩阵. 设
a 1n ⎫⎛a 11a 12
⎪⎛
⋯a a a a 2n ⎪2122
A = = 11
(1) ⎪ a
⎪⎝21⎝a n 1a n 2 a nn ⎭
()⎫A 12
⎪, (1)⎪A 22⎭
1
(1)(1)
不妨设a 21,其中 ≠0,否则这一步不需约化,选择初等矩阵R 1使R 1a 21
1n
⎧22⎪σ1=sgn (a 21)(∑a i 1) ,
i =2⎪
(1)⎪u 1=a 21+σ1e 1,
(2.1) ⎨
12⎪ρ=u ⎪1212=σ1(σ1+a 21), ⎪T
R 1=I -ρ1-1u 1u 1. ⎩
⎛10⎫
令U 1= ⎪则
0R ⎝1⎭
⎛a 11
A 2=U 1AU 11= R a (1)
⎝121
2
(1)(2)A 21R 1⎫⎛A 11
⎪= (1)⎪ 0R 1A 22R 1⎭⎝
(2)
a 12(2)a 22
()⎫A 13
⎪, (2)⎪A 23⎭
2
()
其中A 11∈R 2⨯1, 一直重复这一过程得;
定理
2.1 如果A =R n ⨯n ,则存在初等反射矩阵, 使
U n -2⋯U 2U 1AU 1U 2⋯U n -2=C (上Hessenberg ).
推论A 为对称矩阵,则存在初等反射矩阵 其中A n -1为
⎛c 1b 1⎫ ⎪b c b 122 ⎪ ⎪ A n -1= ⎪
⎪
b n -2c n -1b n -1⎪ ⎪ ⎪b c n -1n ⎭⎝
2.2 矩阵的QR 分解
有了以上的讨论我们知道任何一个方阵A 都可以正交相似化为Hessenberg 或更为特殊的对称三角矩阵,接下来便要使用矩阵的QR (Q 为正交矩阵,R 为三角矩阵)分解来解决Hessenberg 阵和对称角矩阵的全部特征值问题.
引理1 设, 其中αi ,αj 不全为零,则可选取一平面矩阵P ij 使
P ij x =y =(α1, α2, ⋯, αi (1), ⋯, α(j 1), ⋯, αn ) T , 其中
αi (
1)=(2.2)
(2.3)
⎧
c =⎪⎪
(2.4) ⎨
α⎪s =⎪⎩
证明 事实上P ij x 只改变x 的第i 个及第j 个元素,且有
αi (1)=c αi +s αj ,,
于是可选P 即按式(2.4)求c , s ,且有式(2.2)及(2.3). ij
定理2.4 如果A 为非奇异矩阵,则即存在平面旋转矩阵P 1, P 2, ⋯, P n -1(一系列平面旋转矩阵)使
且r ii >0(i =1,2, ⋯, n -1).
证明 由于A 第一列一定存在a j 1≠0,于是,如果a j 1≠0(j =2,3, ⋯, n ), 应用引理1,即存在平面旋转矩阵P 12, P 13, ⋯, P 1n , 使
2(2)⎛r 11a 12a 1(n )⎫
⎪(2)(2)
0a 22⋯a 2n ⎪(2)
P ⋯P P A =1n 1312 ⎪=A
⎪
0a (2) a (2)⎪
n 2nn ⎭⎝
且记P 1n ⋯P 13P 12=P 1,同理重复上述过程,最后得到:存在正交矩阵P 1, P 2, ⋯, P n -1使式(2.5)成立.
定理2.5(矩阵的QR 分解) 如果n 阶方阵A 为非奇异矩阵,则A 可分解为一正交矩阵Q 与一上三角矩阵R , 且当R 对角元素都为正数时分解唯一. 证明 由定理2.4知,存在正交矩阵P 1, P 2, ⋯, P n -1, 使
(2.6)
为上三角阵,记
Q T =P n -1P n -2⋯P 1
T T T
于是式(2.6)为Q T A , 即A =QR , 其中Q =P 为正交矩阵. P ⋯P 12n -1
现证明唯一性. 设有A =Q 1R 1=Q 2R 2,其中R 1,R 2为上三角矩阵(显然为非奇
异矩阵)且对角元素都为正. Q 1, Q 2为正交矩阵. 于是
T
Q 2Q 1=R 2R 1-1(2.7)
由式(2.7)知上三角阵R 2R 1-1为正交矩阵,故R 2R 1-1为对角矩阵,即
R 2R 1-1=D =diag (d 1, d 2, ⋯d n )
因为R 2R 1-1是正交矩阵,所以D 2=I , 又因为R 1,R 2对角元素都为正数,故
d i >0(i =1,2, ⋯, n ),即D =I 于是R 1=R 2,由式(2.6)得到Q 1=Q 2.
2.3 矩阵特征值QR 算法理论依据
设A =A 1=(a ij )∈R n ⨯n 非奇异, 且对A 直接QR 分解,即A =QR , 其中R 为上三角矩阵, Q 为正交矩阵,于是得到一个新矩阵 B =RQ =Q T AQ . 显然, B 是由A 经过正交相似变换得到, 因此A 和B 有相同的的特征值, 在对B 进行QR 分解,又可得到一新矩阵, 重复可得矩阵序列.
设A =A 1, 进行QR 分解, 得A 1=Q 1R 1, 作矩阵A 2=R 1Q 1=Q 1T A 1Q 1 ,,求得A K 后
T 将其进行QR 分解, 得A k =Q k R k , 作矩阵A k +1=R k Q k =Q k A k Q k 而利用QR 分解求
矩阵的特征值, 就是利用矩阵QR 分解, 按上述递推法则构造矩阵序列{A k }的过程, 只要A 非奇异, 利用上述方法就能构造出{A k }.
定理2.6(基本的QR 方法) 设A =A 1=(a ij )∈R n ⨯n 非奇异,QR 基本方法为:
T
⎧R k 为上三角矩阵), ⎪A k =Q k R k (Q k Q k =I
⎨
A k +1=R k Q (⎪k k =1,2, ⋯), ⎩
且记=Q 1⋯Q k ,则有;
T
a:A k +1相似与A k 即 A k +1=Q k A k Q k T b: A k +1=(Q 1⋯Q k )A 1(Q 1⋯Q k )
c: A k 的QR 分解式为 A k
T 证明 a:A k +1=R k Q k =Q k A k Q k
b:
T A k +1=Q k A k Q k
c: 显然,k =1 时有A 1,设A k -1有分解式,于是
=A .
结合定理2.5,可以得出A k +1可由A k 按下述方法求得.
T T T
(1)左变换P n -1P n -2⋯P 1A k =R k (上三角阵) (2)右变换R k P 1P 2⋯P n -1=A k +1
引理2 设M k , 其中Q k ,R k 为具有正对角元素的上三角矩阵,如果
M k →I (k →∞), 则Q k →I , 及R k →I .
(k )T T T
证明 设R k R k =M k M k →I (k →∞) , 记R k =(r ij ) , 矩阵R k R k
()()()
, r 11∙(r 11, r 12, ⋯, r 1(n ))
k
k
k
k
(k )
因此有 r 11,,, (2.8)
()()()()()
r 12∙(r 11, r 12, ⋯, r 1(n ))+r 22∙(0, r 22, ⋯, r 2(n ))
k
k
k
k
k
k
k
利用式(2.8)的结果,则有
(k )
,,, 0(k →∞) r 22
-1T
对于R k R k 其他行同理可得, 故R k , 且易知有R k →I (k →∞), 因此
Q k =M k R k -1→I (k →∞).
定理2.7 设A =(a ij )∈R n ⨯n 如果A 的特征值满足:;
A 有标准型A =XDX -1其中(), 且设X -1有三角分解(L 为单位下三角矩阵,U
为单位上三角矩阵),则有QR 算法产生{A k }本质上收敛于上三角阵,即
(k →∞)
这个定理证明过程比较繁琐,具体证明内容参考文献[2],这里需要说明在证明过
程中A k 的收敛速度依赖r n =λn /λn -1, 当r n 越小时, 收敛越快.
推论2.1 如果A 是对称矩阵满足定理2.7条件, 则由QR 方法产生{A k }本质上收敛于对角矩阵.
不管是上三角矩阵,还是对角矩阵,其特征值得计算是很方便的,注意到对角矩阵或上三角矩阵是极限形式,故在给定的误差范围内,计算满足定理2.7的矩阵特征值是比较容易的. 以上对于理论知识的解释完毕,接下来看一下如何具体在应用这一原理计算矩阵特征值.
2.3 矩阵特征值QR 算法实际操作
(k )
在定理2.7的证明中,我们发现a nn (k →∞)速度很大程度依赖于比值
当r n 很小时,收敛很快. 如果s 为r n 一个估计,且对A -sI 应用QR r n =λn /λn -1, ,
方法,收敛速度恰当,收敛速度会更快.
由于QR 计算矩阵特征值是一种利用极限近似的方法,故如果{A k }收敛的越快,计算将会更方便. 为此,为了加速收敛,选择数列{s k }, 按照下述方法构造矩阵序列{A k }.
Shep 1 设A =A 1=(a ij )∈R n ⨯n
Shep 2 将A k -s k I 进行QR 分解,即A k -s k I , k =1,2, ;
T
Shep 3 构造新矩阵A k +1=R k Q k +s k I =Q k A k Q k ;
Shep 4A k +1 ,其中 ;
Shep 5 矩阵(A -s 1I )(A -s 2I ) ⋯(A -s n I ) =φ(A )有QR 分解式
φ(A );
Shep 6 结合定理2.5的结论:首先利用正交变换(左变换)将A k -s k I 化为上
三角阵, 即
P n -1P n -2⋯P 1(A k -s k I )
T
其中Q k 为一系列平面旋转矩阵的乘积. 于是
A k +1(2.9)
有前面的讨论我们知道任何一个方阵都可以正交相似化为一个上Hessenberg 矩阵,而有正交过程中矩阵特征值不变性质,故下面我们只考虑Hessenberg 矩阵的特质值得计算.
设A =A 1为上Hessenberg g 结合(2.1)的具体P i 的构造方法shep 1~6, 构造出上Hessenberg 矩阵序列{A k },在构造的过程中对{s k }的选择也是十分重要的,
(k )一般选取s k =a nn . 具体方法我们看下面这个例题的解决.
例2.1 用QR 方法计算对称三角矩阵的全部特征值
⎛210⎫
⎪
A =A 1= 131⎪
014⎪⎝⎭
解 首先选取s 1 结合2.1的方法构造旋转平面
⎛2.2361-1.3420.4472⎫ ⎪
P 23P (A -s I )=R =01.0954-0.36511211 ⎪
0-0.36510.81650⎪⎝⎭0⎫⎛1.40000.4899 ⎪T T
A 2=RP 3.26670.7454⎪ 12P 23+s 1I = 0.4899
00.74544.3333⎪⎝⎭
重复以上的操作
0⎫⎛1.29150.2017
⎪A 3= 0.20173.02020.2724⎪
00.27244.6884⎪⎝⎭0⎫⎛1.27370.0993 ⎪A 4= 0.09932.99430.0072⎪
00.00724.7320⎪⎝⎭
0⎫⎛1.26940.0498
⎪A 5= 0.04982.99860⎪
004.7321⎪⎝⎭
记
现在收缩,继续对A 5的子矩阵进行变换,得到
-s 5I
故求得A 近似特征值为
1.2680
且A 的实际特征值是
=3.0 =
结束语
以上便是从矩阵的两个常见分解LU , QR 分解出发,进而把它们应用到解线性方程组,以及求矩阵的特征值,这种方法在解决阶数比较高的矩阵时利用这种思想设计出相应的的算法,往往能十分有效的减缓繁重的计算过程.
参考文献
[1]王萼芳,石生明. 高等代数[M].北京:高等教育出版社,1978.162-203 [2]李庆阳,王能超,易大义. 数值分析[M].武汉:华中科技大学出版社,1982.173-245
[3]周金土. 高等代数解题思想和方法[M].杭州:浙江大学出版社,2008.96-130
[4]徐仲,陆全,张凯院. 高等代数考研教案[M].西安:西北工业大学出版社,2005.155-198
[5]郑崇友,樊磊,崔宏斌.Frame 与连续格[M].北京:首都师范大学出版社,1994
[6]Gierz G ,Hoffmann K H,et al.A Compendium of Continuous Lattice[M].Berlin:Springer-Verlag,1980.
[7]吴耀强. 关于代数L-domain 的一个刻画定理. 江西师范大学学报(自然科学版)[J].2006,30(6):573-576.