程可欣,彭振赟,杜丹丹,肖憲偉

(桂林電子科學大學數(shù)學與計算科學學院廣西高校數(shù)據(jù)分析與計算重點實驗室，桂林 541004)

1 引言

非線性矩陣方程

其中A∈Rn×n，Q為n×n階正定陣，在控制理論、動態(tài)規(guī)劃、插值理論和隨機濾波等領域中具有廣泛的應用[1-3]．例如，很多具有實際應用的偏微分方程問題進行離散化處理后，需要求解線性方程組[4]M x=f，其中

其中Q是對稱正定矩陣．若將系數(shù)矩陣M分解為

則對稱正定矩陣X必須是非線性矩陣方程(1)的解．一旦將系數(shù)矩陣M分解為上述形式，那么，利用Woodbury公式[5]就可以有效地求得線性方程組M x=f的解．

關于這類非線性矩陣方程問題的研究文獻有很多，如Zhan和Xie[1]給出了方程X+ATX?1A=I有正定解的條件及計算極解的無求逆迭代算法．Engwerda[6,7]給出了矩陣方程X+A?X?1A=Q正定解存在的條件，并提出了計算其極解的遞歸算法．Ferrante,Levy[2]和Mein[8]分別研究了矩陣方程X=Q+NX?1N?和X ?A?X?1A=Q矩陣方程，給出了計算其極解的不動點迭代算法．劉新國[9]給出了求解矩陣方程I=X+AHX?1A的簡單迭代并進行誤差估計分析．Ramadan和EIShazly[10]給出了方程有正定對稱解的充要條件，以及計算其解的不動點迭代算法．Peng等[11]給出了計算方程X+A?X?αA=Q(0< α≤1)的極解的無求逆迭代算法，并討論了算法的收斂性問題．段雪峰和廖安平[12]給出了矩陣方程X+A?X?qA=Q(q≥1)有解的條件同時給出了解的誤差估計式．Liu和Zhang[13]研究了利用牛頓迭代法求解矩陣方程X=Q+AH(I?X?C)?1A，并給出了解的誤差估計．關于矩陣方程X±A?X?1A=Q的解的擾動問題，在文獻[14–16]中進行了比較系統(tǒng)地討論．

本文研究矩陣方程X?ATX?1A=Q的牛頓迭代解法．在給定初值的條件下，證明了迭代方法產生的矩陣序列包含在確定的閉球內并收斂到閉球內唯一的矩陣方程的解，同時給出了近似解與真解的誤差估計式．此外，還給出了說明算法有效性的數(shù)值例子．

記號：Rn×n表示階矩陣的集合，I表示n階單位陣，AT表示矩陣A的轉置，∥A∥F表示矩陣A的Frobenius范數(shù)，∥A∥表示矩陣A的譜范數(shù)，A?B表示矩陣A和B的K ronecker積．

2 牛頓迭代法及其收斂性

將正定矩陣Q做Cholesky分解，即Q=LTL，其中L∈Rn×n為n×n階上三角陣，則矩陣方程(1)轉化為

令Y=(LT)?1X L?1，則有X=LTYL．因而矩陣方程(1)等價于

其中P=(LT)?1AL?1∈ Rn×n．

記F(Y)=Y?PTY?1P?I，則求矩陣方程(2)的解等價于求方程F(Y)=0的解．注意到矩陣函數(shù)F(Y)在Y處方向為的Fréchet-導數(shù)為

牛頓法求解非線性矩陣方程(2)的迭代格式可以描述為如下算法1．

算法1(牛頓法求矩陣方程(2)的迭代格式)

第1步給出初始矩陣Y0=I，誤差容許值ε>0．令k:=0；

第2步如果∥F(Yk)∥F≤ε停止；否則，求Ek∈Rn×n，使得；即求Ek∈Rn×n，使得E+BTEB=C，其中

第3步計算Yk+1=Yk+Ek；

第4步置k=k+1，轉第2步．

下面將討論算法1的有關收斂性結果．首先給出如下引理．

引理1[13,17]設S,T ∈ Cn×n．若S可逆，且∥S?1∥≤β,∥S?T∥≤α,αβ<1，則T也可逆，且．

引理2對任意可逆矩陣X,Y有下列不等式成立

證明 對任意可逆矩陣X,Y和任意矩陣E有

因此

引理3假設Y0可逆，且有，則有

證明 由，可知

上式表明，當且僅當E=0．因此，矩陣

是可逆的，并且有

定義矩陣函數(shù)

則

是恒等算子，且關于算子H(Y)的牛頓法生成的矩陣序列{Yk}與關于算法1生成的矩陣序列{Yk}是相同的，并且有如下引理．

引理4假設

則當Y0=I時，下列結論成立：

證明 (a)由于

則由引理3及(4)式，有

(b)因為Y∈B(Y0,σ)，則有∥Y?Y0∥<σ．注意到∥Y0∥=1，則由引理1可得．同理有．因此，由引理2，有

(c)因為，且對任意的Y∈B(Y0,δ)，有

所以由引理1可證得

(d)由(b)式可知

因此，有

對于算法1，有如下收斂性定理．

定理1設∥P∥2=σ．若

則由算法1生成的矩陣序列收斂于F(Y)在閉球B(Y0,δ)內的唯一零點Y?，且有．

證明 首先利用歸納法證明如下結論(e)–(h)成立：

1)當k=1時

(e)和(f)此時成立．由引理4中結論(c)，可得

此時(g)成立．由引理4中的結論(d)，可得

此時(h)成立．

2)假設

均成立，則有

從而有

進而有

所以有

綜上可知，結論(e),(f),(g),(h)均成立．

其次證明{Yk}收斂于，且,?k>1．由(e)式，對任意的k,l≥1，有

因此可知，矩陣序列{Yk}是柯西列．又由于，故存在，使得．因此，矩陣序列{Yk}收斂于Y?．注意到

且H(Y)是關于Y的連續(xù)函數(shù)，故有

因此Y?是F(Y)在內的零點．此外，由

可得

最后證明Y?是F(Y)在閉球內的唯一零點．設Z?∈B(Y0,δ)是H(Y)的零點，即H(Z?)=0．下證,?k≥1．

顯然．

假設,? k=1,2,…,n成立，則

即．因此

即有Y?=Z?，亦即Y?是F(Y)在閉球B(Y0,δ)內的唯一零點．

綜上所述，定理1得證．

注意到矩陣方程(1)和(2)的等價關系，牛頓法求解非線性矩陣方程(1)的迭代格式可以描述為如下算法．

算法2(牛頓法求矩陣方程(1)的迭代格式)

第1步給出初始矩陣X0=Q，誤差容許值ε>0．令k:=0；

第2步如果ε，停止；否則，求Ek∈Rn×n，使得

其中

第3步計算Xk+1=Xk+Ek；

第4步置k=k+1，轉第2步．

對于算法2，有如下收斂性定理．其證明與定理1類似，故省略．

定理2設∥(LT)?1AL?1∥2= σ．若

則由算法2生成的矩陣序列收斂于矩陣方程(1)在閉球B(X0,δ)內的唯一零點X?，且有．

3 數(shù)值例子

本節(jié)將給出說明算法有效性的數(shù)值例子．對于算法2中的子問題(5)，算例中采用如下LSQRM算法[18]進行計算．

LSQRM算法(LSQR算法計算子問題(5)的近似解Ek)：

步驟1給定初始矩陣計算E0=0：

步驟2對于k=1,2,…，執(zhí)行步驟2.1至步驟2.6：

步驟2.1計算

步驟2.2計算

步驟2.3令

步驟2.4計算

步驟2.5若?kαk|ck|< ε，則運行終止，Ek即為問題(6)的解；否則，轉步驟2.6；

步驟2.6令k=k+1，轉步驟2.1．

例1給定矩陣A=randn(9,9),Q=randn(9,9)(Matlab格式)如下：

則有，因此，矩陣方程(1)在閉球

內有唯一解．利用算法2迭代5步，得到矩陣方程(1)的近似解為

并且

參考文獻：

[1]Zhan X Z,Xie J J.On the matrix equation X+ATX?1A=I[J].Linear Algebra and its Applications,1996,247(1):337-345

[2]Ferrante A,Levy B C.Hermitian solutions of the equation X=Q+NX?1N?[J].Linear Algebra and its Applications,1996,247(1):359-373

[3]Helton J W,Sakhnovich L A.Extremal problems of interpolation theory[J].Rocky Mountain Journal of Mathematics,2005,35(3):819-841

[4]Buzbee B L,Golub G H,Nielson C W.On direct methods for solving Poisson’s equations[J].SIAM Journal on Numerical Analysis,1970,7(4):627-656

[5]Housholder A S.The Theory of Matrices in Numerical Analysis[M].New York:Blaisdell,1964

[6]Engwerda J C.On the existence of a positive definite solution of the matrix equation X+ATX?1A=I[J].Linear Algebra and its Applications,1993,194(15):91-108

[7]Engwerda J C,Ran A M,Rijkeboer A L.Necessary and sufficient conditions for the existence of a positive definite solution of the matrix equation X+A?X?1A=Q[J].Linear Algebra and its Applications,1993,186:255-275

[8]Meini B.Efficient computation of the extreme solutions of X+A?X?1A=Q and X ?A?X?1A=Q[J].Mathematics of Computation,2002,71(239):1189-1204

[9]劉新國.關于求解矩陣方程X=X+AHX?1A的簡單迭代[J].高等學校計算數(shù)學學報，1998,(2):163-167 Liu X G.On the simple iteration for solving matrix equation X=X+AHX?1A[J].Numerical Mathematics:A Journal of Chinese Universities,1998,(2):163-167

[10]Ramadan M A,EI-Shazly N M.On the matrix equation X+AT2m√X?1A=I[J].Applied Mathematics and Computation,2006,173(2):992-1013

[9]劉新國.關于求解矩陣方程X=X+AHX?1A的簡單迭代[J].高等學校計算數(shù)學學報，1998,(2):163-167 Liu X G.On the sim p le iteration for solvingm atrix equation X=X+AHX?1A[J].Num ericalM athem atics:A Jou rnal of Chinese Universities,1998,(2):163-167

[10]Ram adan M A,EI-Shazly N M.On the m atrix equation.App lied M athem atics and Com pu tation,2006,173(2):992-1013

[11]Peng Z Y,EI-sayed S M,Zhang X L.Iterative methods for the extremal positive definite solution of the matrix equation X+AX?αA=Q[J].Journal of Computational and Applied Mathematics,2007,200(2):520-527

[12]段雪峰，廖安平.矩陣方程X+A?X?qA=Q(q≥1)的Hermitian正定解及其擾動分析[J].高等學校計算數(shù)學學報，2008,30(3):280-288 Duan X F,Liao A P.The Hermitian positive definite solution of matrix equation X+A?X?qA=Q(q≥1)and its perturbation analysis[J].Numerical Mathematics:A Journal of Chinese Universities,2008,30(3):280-288

[13]Liu P P,Zhang S G.Newton’s method for solving a class of nonlinear matrix equations[J].Journal of Computational and Applied Mathematics,2014,256:254-267

[14]Hasanov V I,Ivanov I G.On two perturbation estimates of the extreme solutions to the equations X±A?X?1A=Q[J].Linear Algebra and its Applications,2006,413(1):81-92

[15]Hasanov V I.Notes on two perturbation estimates of the extreme solutions to the equations X±A?X?1A=Q[J].Applied Mathematics and Computation,2010,216(5):1355-1362

[16]Hasanov V I,Ivanov I G,Uhlig F.Improved perturbation estimates for the matrix equations X±A?X?1A=Q[J].Linear Algebra and its Applications,2004,379(1):113-135

[17]袁亞湘，孫文瑜.最優(yōu)化理論與方法[M].北京：科學出版社，2004:7-9 Yuan Y X,Sun W Y.Optimization Theory and Methods[M].Beijing:Science Press,2004:7-9

[18]Peng Z Y.A matrix LSQR iterative method to solve matrix equation AXB=C[J].International Journal of Computer Mathematics,2010,87(8):1820-1830