浅谈矩阵的LU分解和QR分解及其应用

浅谈矩阵的LU 分解和QR 分解及其应用基于理论研究和计算的需要，往往有必要把矩阵分解为具有某种特性的矩阵之积，这就是我们所说的矩阵分解.本文将介绍两种常用的矩阵分解方法，以及其在解线性方程组及求矩阵特征值中的应用.1.矩阵的LU 分解及其在解线性方程组中的应用 1.1 高斯消元法通过学习，我们了解到利用Gauss 消去法及其一些变形是解决低阶稠密矩阵方程组的有效方法.并且近些年来利用此类方法求具有较大型稀疏矩阵也取得了较大进展.下面我们就通过介绍Gauss 消去法，从而引出矩阵的LU 分解及讨论其对解线性方程组的优越性. 首先通过一个例子引入：例1,解方程组(1.1)(1. 2)(1.3)解. 1Step (1.1)(2)(1.3)⨯-+ 消去(1.3)中未知数,得到23411x x --=- (1.4)2Shep . (1.2)(1.4)+ 消去(1.4)中的未知数2x有12323364526x x x x x x ++=-=-=-⎧⎪⎨⎪⎩ 显然方程组的解为*x =123⎛⎫ ⎪ ⎪ ⎪⎝⎭上述过程相当于 111604152211⎛⎫ ⎪- ⎪ ⎪-⎝⎭~111604150411⎛⎫ ⎪- ⎪ ⎪---⎝⎭~111604150026⎛⎫⎪- ⎪ ⎪--⎝⎭2-（）+ ()i i r 表示矩阵的行由此看出，消去法的基本思想是：用逐次消去未知数的方法把原方程化为与其等价的三角方程组.下面介绍解一般n 阶线性方程组的Gauss 消去法.设111n n1nn a a a a A ⎛⎫ ⎪=⎪ ⎪⎝⎭ 1n x X x ⎛⎫ ⎪= ⎪ ⎪⎝⎭ 1n b b b ⎛⎫ ⎪= ⎪ ⎪⎝⎭则n 阶线性方程组AX b = (1.5)并且A 为非奇异矩阵.通过归纳法可以将AX b =化为与其等价的三角形方程，事实上： 及方程(1.5)为()()11A X b =,其中()1A A = ()1b b =(1) 设(1)110a≠，首先对行计算乘数()()11i1111i a m m =.用1i m -乘(1.5)的第一个方程加到第()2,3,,i i n =⋯个方程上.消去方程(1.5)的第2个方程直到第n 个方程的未知数1x .得到与(1.5)等价的方程组()()()11n 12n 111nn 0a a x x a ⎛⎫⋯⎛⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⋯⎝⎭⎝⎭=()()112n b b ⎛⎫ ⎪⎪ ⎪⎝⎭简记作()()22Ab = (1.6)其中()()()()()()211211111 ijij i ij i i i a m b b m a a b =-=- (2) 一般第()11k k n ≤≤-次消去，设第1k -步计算完成.即等价于()()k k AX b = (1.7)且消去未知数121,,,k x x x -⋯.其中()()()()()()()()()()1111112122222k k k k kk knk nknna n n a a a a a A a a a a ⎛⎫⎪ ⎪⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎝⎭设()0k kk a ≠计算()() (i=/1,,)k k ik ikkkaa k m n =+⋯,用()()(1,,)a n ikik k kki k n a m ==+⋯消去第1k +个方程直到第n 个方程的未知数k x .得到与(1.7)等价的方程组()()1k 1k A X b ++= 故由数学归纳法知，最后可以把原方程化成一个与原方程等价的三角方程组.但是以上分析明显存在一个问题，即使A 非奇异也无法保证()0i ii a ≠,需要把非奇异的条件加强.引理1 约化主元素()01,,)i ii a k ≠=⋯（i 的充要条件是矩阵A 的顺序主子式0i D ≠.即1111110,0ikk kkk a a D a a D a =≠=≠⋯证明 利用数学归纳法证明引理的充分性.显然，当1k = 时引理的充分性是成立的，现在假设引理对1k -是成立的，求证引理对k 亦成立.有归纳法，设()()01,21iii i a k ≠=⋯-于是可用Gauss 消去法将中，即()()()()()()()()()()()11111121n22222n 1k k k k k kk knnknn a a a a a A a a a a A ⎛⎫ ⎪ ⎪⎪ ⎪→= ⎪ ⎪ ⎪ ⎪⎝⎭即()()()()()()()()11121231112211223112233222a a D a D a a a a a ===()()()()()()()()()11111122212222k 11122k k k kkk kka a a a a a D a a a ==⋯ (1.8) 由设0(1,,)i i D k ≠=⋯及式(1.8)有()0k kk a ≠显然，由假设()()01,2iiii k a ≠=⋯，利用(1.8)亦可以推出0(1,,)i i D k ≠=⋯ 从而由此前的分析易得；定理1 如果n 阶矩阵A 的所有顺序主子式均不为零，则可通过Gauss 消去法（不进行交换两行的初等变换），将方程组(1.5)约化成上三角方程组，即()()()()()()()()()1111111121122222222b b b n n n n nn n n a a a x x a a x a ⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪⎝⎭⎝⎭⎝⎭ (1.9) 1.2 矩阵LU 分解从而由以上讨论即能引出矩阵的LU 分解，通过高等代数我们得知对A 施行行初等变换相当于用初等矩阵左乘A ，即()()()()121211L A Lb A b == 其中 211n11101L m m ⎛⎫⎪- ⎪= ⎪⎪-⎝⎭一般第k 步消元，，相当于()()()()11k kk k k kL A A L b b ++==重复这一过程，最后得到()()()()11211121n n n n Ab L L L A L L L b --⎧⋯=⎪⎨⋯=⎪⎩ (1.10) 其中k 1,111m 1n k k km L +⎛⎫ ⎪ ⎪ ⎪=⎪- ⎪ ⎪ ⎪ ⎪-⎝⎭将上三角形矩阵()n A U 记作,由式(1.9)得到111121=U n A L L L LU ----⋯=,其中211111211211m 1n n n m L L m L L ----⎛⎫⎪⎪=⋯= ⎪⎪⎝⎭由以上分析得；定理2 （LU 分解） 设A 为n 阶矩阵，如果A 的顺序主子式i 0(1,2,,1)D i n ≠=-.则A 可分解为一个单位下三角矩阵L 和一个上三角矩阵U的乘积，且这种分解是唯一的.证明 由先前的分析得出存在性是显然的，即A LU =.下证唯一性,设A LU CD == 其中L , C 为单位下三角矩阵，U ,D 为上三角矩阵.由于1D -11D C L U --=上式右端为上三角矩阵，左端为单位下三角矩阵，从而上式两端都必须等于单位矩阵，故U D =,L C =.证毕.例2 对于例子1 系数矩阵矩阵111041221A ⎛⎫ ⎪=- ⎪ ⎪-⎝⎭由Gauss 消去法，得结合例1，故100111010041211002A LU ⎛⎫⎛⎫⎪⎪==- ⎪⎪ ⎪⎪--⎝⎭⎝⎭对于一般的非奇异矩阵，我们可以利用初等排列矩阵kki I (由交换单位矩阵I的第k 行与第k i 行得到)，即()()()()()()()()111212111111,,kk k k k ki k i k k i i k A L b L I A I b L I A I b A L b++⎧==⎪⎨==⎪⎩ (1.11) 利用(1.11)得()1111,11n nn n i i L I L U I A A ---==.简记做.其中下面就n 情况来考察一下矩阵.()()4321444343544332211443443243)(i i i i i i i i i i I L I L I L I A I L I I L I I A A L L I ===⨯4324324321432()i i i i i i I I I L I I I 43214321 )(i i i i I I I I 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 为上三角矩阵。

即A =111122212212111n n n n nn u u u u l u l l u ⎛⎫⎛⎫⎪⎪ ⎪⎪ ⎪ ⎪ ⎪⎪⎝⎭⎝⎭(1.13)下面说明L ,的元素可以由n 步直接计算出： 由(1.13)有11i i a u = (1.14)再由矩阵乘法得1111i i a u l =,故有11i i l a =/11u于是得U 的第一列元素.设已经得U 的第一行到1r -行元素与L 的第一列到1r -列元素，由(1.13)有111nr ri rk ki rk ki ri k k l l u u u a -====+∑∑故有11r ri ri rk ki k a u u l -==-∑ (i ) (1.15)再由(1.13)得111ir ik kr ik kr nr k k ir rr a l u l l u u -====+∑∑ .得11()/ir ir ik k r r r k r a l u l u -==-∑ (1.16)故有以上分析结合(1.12)得；111i i 1ik i k ky y l y b b -==⎧⎪⎨=-⎪⎩∑ 1//n n nn i i ik k n i i ik x x b u x y u u =+=⎧⎪⎛⎫⎨=- ⎪⎪⎝⎭⎩∑ (1.18) 例3． 求解123123142521831520x x x ⎛⎫⎛⎫⎛⎫⎪⎪ ⎪= ⎪⎪ ⎪ ⎪⎪ ⎪⎝⎭⎝⎭⎝⎭解. 由 (1.15)，(1.16)计算，得1001232100143510024A LU ⎛⎫⎛⎫⎪⎪=-= ⎪⎪ ⎪⎪--⎝⎭⎝⎭求解141820Ly ⎛⎫ ⎪= ⎪ ⎪⎝⎭ 得141072y ⎛⎫ ⎪=- ⎪ ⎪-⎝⎭ 求解141072Ux ⎛⎫ ⎪=- ⎪ ⎪-⎝⎭得123x ⎛⎫⎪= ⎪ ⎪⎝⎭以上便是通过介绍Gauss 消去法，引出矩阵的LU 分解，这种分解在数值分析中，在设计算法求解高维线性方程组能提高效率，不像Grammar 法则只是从理论上解决了非奇异线性方程组的解法，实际操作起来是十分不方便的，而利用LU 分解的基本原理能大大减少计算过程的繁琐.2．矩阵的QR 分解及其在计算矩阵特征值的应用 2.1 转化非奇异矩阵为上Hessebberg定义1 一方阵，如果1i j >+时有0ij b =.则称矩阵B 为上Hessebberg 阵，即111212122232,10000b b n nn n n n b b b b b b b B -⎛⎫⎪ ⎪⎪=⎪ ⎪ ⎪ ⎪ ⎪⎝⎭定义2 设向量w 满足2w ,矩阵2T H I ww =-称为初等反射矩阵，记作()H w ,即()21112222122n 121222212122n n n w w w w w w H w w w w w w w w ⋯-⎛⎫-- ⎪⋯---⎪= ⎪ ⎪ ⎪--⎝⎭其中 定理2.1 初等反射阵H 是对称矩阵，正交矩阵和合同矩阵.定理2.2 设,x y 为两个不相等的n 维向量，但22x y =,则存在一个初等反射矩阵H ，使Hx y =.证明 令2()/w x y x y =--,则得到一个初等反射矩阵()()2222/T T T H x x y I I y w y x w --=-=--而且()()()()2T T 22222/TT Hx x x y y x x xx x x x y x y x yy =---=-----因为22(()2)()T T T x y x y y x y x x x -=---=所以()Hx x x y y =--=.并且由代数学知识易知，2()/y w x y x --=.成立的唯一长度等于1的向量.推论 设向量()0n x R x ∈≠,2σx =±,且1σx e ≠-则存在一个初等反射矩阵212,/2I T T H I uu uu u ρ-=-=-使1σe Hx =-,其中1σu e x =+ ,22/2u ρ=设()123,,,,n x αααα=⋯ ()12,,,n u u u u =⋯,则12n (σ,,,)u ααα=+⋯，()212/2σσu ρα==+.如果σ，那么计算1σα+时有效数字可能损失，取有相同符号，即取()12σsgn x α=有了以上关于初等反射矩阵的性质接下来讨论正交相似约化一般矩阵和对角矩阵.设()()()11112122122111211212212nn n n n n A a a a a a a a A a A a a a ⎛⎫ ⎪⎛⎫⋯ ⎪==⎪ ⎪ ⎪⎝⎭⎪⎝⎭, 不妨设()1210a ≠，否则这一步不需约化，选择初等矩阵1R 使1R ()121a ，其中 ()()()122121121121112111121211111σsgn (),σ,1,.σσ2ni i Ta a I a a u e u R u u ρρ=-⎧=⎪⎪⎪=+⎨⎪==+⎪⎪=-⎩∑ (2.1) 令11100U R ⎛⎫= ⎪⎝⎭则()()()()()()()()21221311211121121112112231211221220A a A R a A A AU A R a R A R a U ⎛⎫⎛⎫=== ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭，其中()22111A R ⨯∈,一直重复这一过程得;定理 2.1 如果n n A R ⨯= ,则存在初等反射矩阵,使221122n n U U AU U U U C --⋯⋯=(上Hessenberg ).推论 A 为对称矩阵，则存在初等反射矩阵 其中1n A -为1112212111n n n n n n c b b c b A b c b b c -----⎛⎫⎪ ⎪ ⎪⎪ ⎪ ⎪⎪ ⎪⎝⎭=2.2 矩阵的QR 分解有了以上的讨论我们知道任何一个方阵A 都可以正交相似化为Hessenberg 或更为特殊的对称三角矩阵，接下来便要使用矩阵的QR （Q 为正交矩阵，R 为三角矩阵）分解来解决Hessenberg 阵和对称角矩阵的全部特征值问题.引理1 设,其中i j αα，不全为零，则可选取一平面矩阵ij P 使()()1112(,,,,,,,)Tij i j nP y x ααααα==⋯⋯⋯,其中 ()1i α= (2.2)(2.3)c s α⎧=⎪⎪⎨⎪=⎪⎩(2.4)证明 事实上ij P x 只改变x 的第i 个及第j 个元素，且有()1i i j c s ααα=+,,于是可选ij P 即按式(2.4)求,c s ，且有式(2.2)及(2.3).定理2.4 如果A 为非奇异矩阵，则即存在平面旋转矩阵121,,,n P P P -⋯（一系列平面旋转矩阵）使且()01,2,,1ii r i n >=⋯-.证明 由于A 第一列一定存在10j a ≠，于是，如果()102,3,,j j n a ≠=⋯,应用引理1，即存在平面旋转矩阵12131,,,n P P P ⋯,使()()()()()()()22111122222221131222200nn nn n n a r a a a P P A a a P A ⎛⎫ ⎪⋯ ⎪⋯== ⎪⎪⎪⎝⎭ 且记113121n P P P P ⋯=，同理重复上述过程，最后得到：存在正交矩阵211,,,n P P P -⋯使式(2.5)成立.定理2.5（矩阵的QR 分解） 如果n 阶方阵A 为非奇异矩阵，则A 可分解为一正交矩阵Q 与一上三角矩阵R ,且当R 对角元素都为正数时分解唯一. 证明 由定理2.4知，存在正交矩阵211,,,n P P P -⋯,使(2.6)为上三角阵，记T 112n n P Q P P --=⋯于是式(2.6)为T Q A ,即A QR =,其中121T T n TQ P P P -=⋯为正交矩阵.现证明唯一性. 设有1122A Q R R Q ==，其中12R R ，为上三角矩阵（显然为非奇异矩阵） 且对角元素都为正. 1Q , 2Q 为正交矩阵.于是12121T R Q Q R -= (2.7)由式(2.7)知上三角阵121R R -为正交矩阵，故121R R -为对角矩阵，即()12112,,n R D diag d R d d -==⋯因为121R R -是正交矩阵，所以2D I =,又因为12R R ，对角元素都为正数，故()01,2,,i i n d >=⋯，即D I = 于是12R R =，由式(2.6)得到12Q Q =.2.3 矩阵特征值QR 算法理论依据设()1ij n n A A a R ⨯==∈非奇异,且对A 直接QR 分解，即A QR =,其中R 为上三角矩阵,Q 为正交矩阵，于是得到一个新矩阵 T AQ B RQ Q ==.显然,B 是由A 经过正交相似变换得到,因此A 和B 有相同的的特征值,在对B 进行QR 分解，又可得到一新矩阵,重复可得矩阵序列.设1A A =,进行QR 分解,得111A Q R =,作矩阵211111=T A Q Q A R Q = ,,求得K A 后将其进行QR 分解,得k k k A Q R =,作矩阵1Tk k k k k k R Q A Q A Q +==而利用QR 分解求矩阵的特征值,就是利用矩阵QR 分解,按上述递推法则构造矩阵序列{}k A 的过程,只要A 非奇异,利用上述方法就能构造出{}k A .定理2.6(基本的QR 方法) 设()1ij n n A A a R ⨯==∈非奇异，QR 基本方法为：()k 1,1,2,Tk k k k k k k k A R Q Q A Q IR R Q k +⎧==⎪⎨==⋯⎪⎩为上三角矩阵（）, 且记=1k Q Q ⋯ ,则有;a: 1k A +相似与k A 即 1Tk k k k Q A A Q +=b: 1k A +111()T k k Q A Q Q Q =⋯⋯（）c: k A 的QR 分解式为 k A证明 a ：1T k k k k k k Q Q A R A Q +==b:1Tk k k k Q A A Q +=c: 显然，1k = 时有1A ，设1k A -有分解式，于是=A .结合定理2.5，可以得出1k A +可由k A 按下述方法求得.(1)左变换211n k k n P A R P P --⋯=(上三角阵) (2)右变换1121n T T Tk k P P R P A -+=⋯引理2 设k M ,其中k Q ，k R 为具有正对角元素的上三角矩阵，如果()k M I k ∞→→,则k ,k Q I R I →→及.证明 设T T k k k k M R I R M =→ ()k ∞→,记()()kk ij R r =,矩阵T k k R R()()()()1111121,,,kkkn kr r r r ⋯•（）, 因此有 ()11kr ,,, (2.8)()()()()()()()k121112122222,,,0,,,n n kkkkkkr r r r r r r ⋯+•⋯•（）（）利用式(2.8)的结果，则有()22k r ,,,0()k ∞→ 对于T k k R R 其他行同理可得,故k R ,且易知有()1k R I k ∞-→→,因此k Q =()1k k I R k M ∞-→→.定理2.7 设()ij n n a A R ⨯=∈ 如果A 的特征值满足：；A 有标准型1A XDX -=其中（）,且设1X -有三角分解（L 为单位下三角矩阵，U 为单位上三角矩阵），则有QR 算法产生{}k A 本质上收敛于上三角阵，即()k ∞→这个定理证明过程比较繁琐，具体证明内容参考文献[2],这里需要说明在证明过程中k A 的收敛速度依赖1/n n n r λλ-=,当n r 越小时,收敛越快.推论2.1 如果A 是对称矩阵满足定理2.7条件,则由QR 方法产生{}k A 本质上收敛于对角矩阵.不管是上三角矩阵，还是对角矩阵，其特征值得计算是很方便的，注意到对角矩阵或上三角矩阵是极限形式，故在给定的误差范围内，计算满足定理2.7的矩阵特征值是比较容易的.以上对于理论知识的解释完毕，接下来看一下如何具体在应用这一原理计算矩阵特征值.2.3 矩阵特征值QR 算法实际操作在定理2.7的证明中，我们发现()k nn a ()k ∞→速度很大程度依赖于比值1/n n n r λλ-=,，当n r 很小时，收敛很快.如果s 为n r 一个估计，且对A sI -应用QR 方法，收敛速度恰当，收敛速度会更快.由于QR 计算矩阵特征值是一种利用极限近似的方法，故如果{}k A 收敛的越快，计算将会更方便.为此，为了加速收敛，选择数列{}k s ,按照下述方法构造矩阵序列{}k A .1Shep 设()1ij n n A A a R ⨯==∈2Shep 将k k A s I -进行QR 分解，即k k A s I -,1,2,k =;3Shep 构造新矩阵1Tk k k k k k k A R s Q A Q I Q +=+=;4Shep 1k A + ,其中 ；5Shep 矩阵()12()()()φn A I A I s s A A s I --⋯-=有QR 分解式()φA ；6Shep 结合定理2.5的结论：首先利用正交变换（左变换）将k k A s I -化为上三角阵, 即112n n P P P --⋯()k k A s I -其中Tk Q 为一系列平面旋转矩阵的乘积.于是1k A + (2.9)有前面的讨论我们知道任何一个方阵都可以正交相似化为一个上Hessenberg 矩阵，而有正交过程中矩阵特征值不变性质，故下面我们只考虑Hessenberg 矩阵的特质值得计算.设1A A =为上Hessenberg g 结合(2.1)的具体i P 的构造方法16shep ~,构造出上Hessenberg 矩阵序列{}k A ，在构造的过程中对{}k s 的选择也是十分重要的，一般选取()k s =knn a .具体方法我们看下面这个例题的解决.例2.1 用QR 方法计算对称三角矩阵的全部特征值1210131014A A ⎛⎫⎪== ⎪ ⎪⎝⎭解 首先选取1s 结合2.1的方法构造旋转平面231211 2.2361 1.3420.44720 1.09540.365100.36510.81650P P I A s R -⎛⎫ ⎪-==- ⎪ ⎪-⎝⎭（）212231 1.40000.489900.4899 3.26670.745400.7454 4.3333T TA P I RP s ⎛⎫ ⎪=+= ⎪ ⎪⎝⎭重复以上的操作3 1.29150.201700.2017 3.02020.272400.2724 4.6884A ⎛⎫ ⎪= ⎪ ⎪⎝⎭4 1.27370.099300.0993 2.99430.007200.0072 4.7320A ⎛⎫ ⎪= ⎪ ⎪⎝⎭51.26940.04980 0.0498 2.9986000 4.7321A⎛⎫ ⎪= ⎪ ⎪⎝⎭记现在收缩，继续对5A的子矩阵进行变换，得到-5s I故求得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.。