矩陣方程X?ATX?1A=Q的牛頓迭代解法?

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

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

1 引言

非線性矩陣方程

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

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

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

關(guān)于這類非線性矩陣方程問題的研究文獻(xiàn)有很多，如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的簡單迭代并進(jìn)行誤差估計分析．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，并給出了解的誤差估計．關(guān)于矩陣方程X±A?X?1A=Q的解的擾動問題，在文獻(xiàn)[14–16]中進(jìn)行了比較系統(tǒng)地討論．

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

記號：Rn×n表示階矩陣的集合，I表示n階單位陣，AT表示矩陣A的轉(zhuǎn)置，∥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)轉(zhuǎn)化為

令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-導(dǎo)數(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，轉(zhuǎn)第2步．

下面將討論算法1的有關(guān)收斂性結(jié)果．首先給出如下引理．

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

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

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

因此

引理3假設(shè)Y0可逆，且有，則有

證明 由，可知

上式表明，當(dāng)且僅當(dāng)E=0．因此，矩陣

是可逆的，并且有

定義矩陣函數(shù)

則

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

引理4假設(shè)

則當(dāng)Y0=I時，下列結(jié)論成立：

證明 (a)由于

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

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

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

所以由引理1可證得

(d)由(b)式可知

因此，有

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

定理1設(shè)∥P∥2=σ．若

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

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

1)當(dāng)k=1時

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

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

此時(h)成立．

2)假設(shè)

均成立，則有

從而有

進(jìn)而有

所以有

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

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

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

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

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

可得

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

顯然．

假設(shè),? k=1,2,…,n成立，則

即．因此

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

綜上所述，定理1得證．

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

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

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

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

其中

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

第4步置k=k+1，轉(zhuǎn)第2步．

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

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

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

3 數(shù)值例子

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

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)的解；否則，轉(zhuǎn)步驟2.6；

步驟2.6令k=k+1，轉(zhuǎn)步驟2.1．

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

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

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

并且

參考文獻(xiàn)：

[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]劉新國.關(guān)于求解矩陣方程X=X+AHX?1A的簡單迭代[J].高等學(xué)校計算數(shù)學(xué)學(xué)報，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]劉新國.關(guān)于求解矩陣方程X=X+AHX?1A的簡單迭代[J].高等學(xué)校計算數(shù)學(xué)學(xué)報，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].高等學(xué)校計算數(shù)學(xué)學(xué)報，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].北京：科學(xué)出版社，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