李勝利，王川龍，溫瑞萍

（1.太原理工大學(xué) 數(shù)學(xué)學(xué)院，山西 太原 030024；2.工程科學(xué)計(jì)算山西省重點(diǎn)實(shí)驗(yàn)室，山西 太原 030012）

1 預(yù)備知識(shí)

線性方程組

是非奇異矩陣，其解法已有了廣泛的研究，主要包括直接法［1-3］和迭代法［4-5］.其中直接法包括高斯消去法，Gram-Schmidt正交化QR方法，雙對角化方法等.單從浮點(diǎn)數(shù)（Flop Numbers）角度看，高斯消去法是解線性方程組（1）最經(jīng)濟(jì)的方法［1-2，6-7］.但是高斯消去法數(shù)值穩(wěn)定性較差［7］，我們不得不考慮“增長因子”的影響.Gram-Schmidt正交化QR［3，8-11］與雙對角化方法［6］對于病態(tài)問題更可靠，但是計(jì)算量又過大.本文提出的正交降階方法既能保證較好的數(shù)值穩(wěn)定性，又比Gram-Schmidt正交化QR以及雙對角化方法所需計(jì)算量要少.

為方便，用Rn×n表示n階實(shí)矩陣全體，Rn表示n維實(shí)向量全體，AT和xT分別表示矩陣A和向量x的轉(zhuǎn)置，（x，y）表示向量x和向量y的內(nèi)積.A＝［v1，v2，…，vn］是A∈Rn×n的一個(gè)列劃分，即vi＝（a1，i，a2，i，…，an，i）T，i＝1，2，…，n.

下面是本文將要用到的一個(gè)重要結(jié)果.

定理1 設(shè)A∈Rn×n為非奇異矩陣，將A的第1列中第1個(gè)不為零的元素記為ai，1，將A的后n-1 列分別與第1列正交后去掉第1列和第i行剩下的（n-1）×（n-1）矩陣仍然非奇異.

證明 設(shè)A的各列向量分別記為v1，v2，…，vn-1，vn，現(xiàn)在將后n-1列分別和第1列正交化，用Gram-Schmidt正交化得

用矩陣表示為

將矩陣分塊，令p＝（a1，1，a2，1，…，ai-1，1）T，q＝（ai+1，1，ai+2，1，…，an，1）T，r＝（ai，2，ai，3，…，ai，n）T，l＝（l1，2，l1，3，…，l1，n）T，I表示（n-1）×（n-1）單位矩陣，則式（2）變成

經(jīng)第1步變換后得到第1列，與后n-1 列分別正交，于是有

這里0表示n-1 維零向量，所以可得到

所以向量-ai，1lT+rT可由-plT+B和-qlT+C線性表示.又因?yàn)閷的后n-1列分別與第1列正交后的矩陣的秩仍然為n，去掉第1列后剩余的矩陣的秩為n-1.所以剩余矩陣的行秩也為n-1.由式（3）知第i行可由其余行線性表示，則去掉第i行后剩余的n-1行線性無關(guān)，所以剩余的n-1 階矩陣非奇異.

2 正交降階分解算法

2.1 算法的構(gòu)造

為求解線性方程組（1），記A的第1列中第1個(gè)不為零的元素為ai，1，先將矩陣A的第1列v1與后n-1列v2，…，vn分別正交，這等價(jià)于將矩陣A右乘矩陣P，于是線性方程組變成下列的等價(jià)形式：

對矩陣A＝［v1，v2，…，vn-1，vn］∈Rn×n，現(xiàn)在將后n-1列v2，…，vn分別和第1列v1正交化，將Gram-Schmidt正交化后A的各列記為則原方程變?yōu)?/p>

將矩陣A第1次正交化后的n階矩陣去掉第1列和第i行，剩下的n-1階非奇異方陣記作A2.由定理1知，線性方程組（2）的系數(shù)矩陣中第i行與其它行線性相關(guān)，去掉第i行的n-1階非奇異方陣即為A2.將b-v1y1去掉第i個(gè)元素后的向量記作，于是n階線性方程組（1）轉(zhuǎn)化為n-1階線性方程組，以上過程稱為正交降階的過程.對得到的新的低階線性方程組重復(fù)以上步驟，直到降為一階線性方程.在重復(fù)正交降階的過程中分別設(shè)yk＝xk+lk，k+1xk+1+…+lk，nxn，k＝1，…，n.最終可求解得到y(tǒng)＝（y1，y2，…，yn）T.因而有

最終將一個(gè)n階線性方程組轉(zhuǎn)化為一個(gè)n階上三角線性方程組，利用回代法可求得其解.

算法1 正交降階分解算法

第1步：令k＝1，＝x.

第2步：記A的第1列中第1個(gè)不為零的元素為ai1，將矩陣A的第1列與其它列分別正交（也就是對A作列變換），正交后的各列記作v1，v2，…，vn-k+1，得到即v1y（k）+v2xk+1+…+vn-k+1xn＝b，其中y（k）是xk，xk+1，…，xn的線性組合，且各項(xiàng)系數(shù)構(gòu)成R的第k行元素.

第3步：利用各列與第1 列的正交性，求得y（k），得到v2xk+1+…+vn-k+1xn＝b-v1y（k），去掉其系數(shù)矩陣的第i行，去掉右端b-v1y（k）的第i個(gè)元素，得到新的n-k階方程.系數(shù)矩陣、未知向量和右端項(xiàng)仍記作和b.k＝k+1，若k＜n返回第2步，若k＝n則進(jìn)入第4步.

第4步：得到等價(jià)上三角線性方程組Rx＝y(tǒng)，求解上三角線性方程組.

2.2 算法的復(fù)雜性分析

下面考慮正交降階分解算法的運(yùn)算量.運(yùn)用正交降階法，運(yùn)算量為個(gè)flops；而利用Gram-Schmidt正交化，將n階系數(shù)矩陣QR 分解需要個(gè)flops.

由表1可知，從運(yùn)算量來看，高斯消去法是解線性方程組最經(jīng)濟(jì)的方法.盡管如此，有下面理由說明為什么仍然要考慮正交化方法.

對于病態(tài)問題，正交化方法更加可靠，也就是說正交化方法能更好地保持?jǐn)?shù)值穩(wěn)定性，不必像高斯消去法那樣擔(dān)心“增長因子”.

表1 各種直接法運(yùn)算量比較Tab.1 Computational comparison of various direct methods

3 數(shù)值實(shí)驗(yàn)

例1 首先在MATLAB 上隨機(jī)生成n階非奇異的矩陣，然后將運(yùn)用Gram-Schmidt正交化進(jìn)行QR 分解，將所占用的CPU 的時(shí)間與正交降階分解算法所占用的CPU 的時(shí)間進(jìn)行對比，結(jié)果如表2 所示.

表2 Gram-Schmidt QR 分解法與算法1CPU 占用時(shí)間對比Tab.2 Comparison of CPU time between Gram-Schmidt QR decomposition method and algorithm 1 （s）

從CPU 的占用時(shí)間來看，當(dāng)n＝500，1 000，2 000，5 000時(shí)，Gram-Schmidt QR 算法大約是算法1所占用時(shí)間的1.23，2.31，2.82，3.61倍，且階數(shù)越大，優(yōu)勢越明顯.

例2 設(shè)矩陣

a＝ones（1 000，1），A＊a＝b，Ax＝b的精確解為x＝ones（1 000，1）.

LU 分解方法和算法1求解求線性方程組（1）的過程對比如圖1 所示.

在圖1中，C表示解的分量，S表示解.

例2 表明與運(yùn)算量最少的高斯消去法相比，運(yùn)算量次少的算法1具有較強(qiáng)的數(shù)值穩(wěn)定性.

例3A取病態(tài)矩陣的典型Hilbert矩陣.A＝hilb（8），x＝ones（8，1），b＝A＊x，此方程的精確解為x＝ones（8，1）.用不同方法求解該方程組，得到的結(jié)果如表3 所示.

表3 算法1與QR 方法的求解數(shù)值結(jié)果Tab.3 Comparison result for algorithm 1and QR decomposition method

例3 表明對病態(tài)Hilbert 矩陣，算法1 較Gram-Schmidt正交化QR 方法數(shù)值更穩(wěn)定.

從上述數(shù)值實(shí)驗(yàn)的結(jié)果可以看出，本文提出的正交降階法不僅運(yùn)算速度得到了很好的改善，運(yùn)算精度也得到了提高，而且具有良好的數(shù)值穩(wěn)定性.

圖1 算法1和LU 分解法的對比Fig.1 Comparison between LU decomposition method and algorithm 1

［1］Golub G H，Van Loan C F.Matrix Computations［M］.Baltimore：The Johns Hopkins University Press，1996.

［2］徐樹方.數(shù)值線性代數(shù)［M］.北京：北京大學(xué)出版社，2011.

［3］徐樹方.矩陣計(jì)算的理論與方法［M］.北京：北京大學(xué)出版社，1995.

［4］Strang G.Linear Algebra and Its Applications［M］.2nd ed.New York：Academic Press，1980.

［5］Ruhe A.Numerical aspects of gram-schmidt orthogonal-ization of vectors［J］.Linear Algebra and its Applica-tions，1983，52-53：591-601.

［6］Giraud L，Langou J.Rounding error analysis of the classical gram-schmidt orthogonalization process［J］.Numer.Math.，2005，101（1）：87-100.

［7］Rice J R.Experiments on gram-schmidt orthogonalization［J］.Math.Comput.，1966，20（94）：325-328.

［8］Daniel J W，Gragg W B，Kaufman L，et al.Reorthogonalizationand stable algorithms for updating the gramschmidt QR factorization［J］.Math.Comput.，1976，136（30）：772-795.

［9］Giraud L，Langou J.When modified gram-schmidt generates a well-conditioned set of vectors［J］.IMA Journal of Numerical Analysis，2002，22（4）：521-528.

［10］Wang C L，Yan X H.On convergence of splitting iteration methods for non-Hermitian positive definite linear systems［J］.International Journal of Computer Mathematics，2013，90（2）：292-305.

［11］Meng G Y，Wang C L，Yan X H.Self-Adaptive Non-Stationary Parallel Multisplitting Two-Stage Iterative Methods for Linear Systems［M］.Berlin Heidelberg：Springer-Verlag，2012：38-47.