高藝萌， 王衛(wèi)國

(中國海洋大學數(shù)學科學學院， 山東 青島 266100)

本文主要考慮Sylvester矩陣方程

AX+XB=C,

(1)

其中A∈Rn×n，B∈Rm×m，C∈Rn×m，且系數(shù)矩陣滿足下列條件之一：

(Ⅰ)A,B是非奇異M-矩陣，C≥0；

(Ⅱ)A,B是不可約M-矩陣且其中至少有一個是非奇異矩陣，C≥0。

方程(1)稱為M-矩陣Sylvester方程(MSE)[1]。當系數(shù)矩陣沒有結(jié)構(gòu)要求時，稱為一般的Sylvester矩陣方程；當A=BT時，稱為Lyapunov矩陣方程。此類廣泛應(yīng)用控制理論領(lǐng)域[2]，海洋平臺減振系統(tǒng)中反饋增益模型需要求解Sylvester方程[3]。

求解方程(1)的常用數(shù)值算法有不動點迭代法[1，4]、牛頓法[2]、交替方向加倍算法[5]、子空間迭代法[6]和保結(jié)構(gòu)加倍算法[2]等。當系數(shù)矩陣A,B階數(shù)n,m很大時，MSE的求解速度會非常慢，且占用大量內(nèi)存。P. Benner等針對大型矩陣代數(shù)Riccati方程，結(jié)合牛頓法和交替方向隱式迭代的思想提出了低秩牛頓交替方向法[7]；D. Kressner等針對線性矩陣方程在低秩情況下提出了簡單而快速的數(shù)值方法，基于這種思想，對系數(shù)矩陣具有特定的低秩分級結(jié)構(gòu)提出了分而治之方法[8]。

本文將考慮系數(shù)矩陣具有分級非對角低秩(Hierarchically off-diagonal low-rank，HODLR)結(jié)構(gòu)M-矩陣Sylvester方程，提出了分而治之方法，把方程化為低階MSE和低秩MSE分別求解。具體思路如下：

設(shè)矩陣A,B,C可以分裂為

A=A0+δA，B=B0+δB，C=C0+δC，

(2)

其中δA,δB,δC是低秩矩陣。假設(shè)X0是矩陣方程

A0X0+X0B0=C0

(3)

的解，求δX使得X=X0+δX是方程(4)的解。

A0+δAX0+δX+
X0+δXB0+δB=(C0+δC)

(4)

由式(3)和(4)得

(5)

故式(5)是右端矩陣為低秩矩陣的方程，當式(5)是M-矩陣Sylvester方程時，我們簡稱為更新MSE問題。當式(3)中已知矩陣都是塊對角陣時，式(3)將轉(zhuǎn)化為多個低階MSE問題的求解。

注1求解方程(1)轉(zhuǎn)化為求解方程(3)和(5)。

1 M-矩陣及相關(guān)性質(zhì)

本節(jié)主要介紹M-矩陣的定義與性質(zhì)，以及對應(yīng)的矩陣方程的可解性。

定義1[1]設(shè)A∈Rn×n，如果所有非對角元素都小于等于0即對?i≠j,aij≤0，稱A是Z-矩陣。如果存在數(shù)s和非負矩陣N，使得A=sI-N且s≥ρ(N)，就稱矩陣A是M-矩陣。當s>ρ(N)時，稱A是非奇異M-矩陣；當s=ρ(N)時，稱A是奇異M-矩陣。

關(guān)于M-矩陣有如下性質(zhì)：

引理1[1](Ⅰ)設(shè)A是Z-矩陣，且存在v≥0，使得Av≥0，則A是M-矩陣；

(Ⅱ)M-矩陣的主子矩陣仍然是M-矩陣；

(Ⅲ)如果A是非奇異M-矩陣，且Z-矩陣B滿足B≥A，則B是非奇異M-矩陣；

(Ⅳ)如果A是不可約奇異M-矩陣，且Z-矩陣B滿足B≥A，當B≠A時，則B是非奇異M-矩陣；

(Ⅴ)如果A是非奇異的M-矩陣或不可約奇異M-矩陣，分塊如下：

其中A11和A22是方陣，則A11和A22都是非奇異的M-矩陣。

引理2[1]設(shè)A∈Rn×n是Z-矩陣，則下列幾個命題等價：

(Ⅰ)A是非奇異M-矩陣；

(Ⅱ)A-1≥0；

(Ⅲ)存在v∈Rn滿足v>0,使得Av>0；

(Ⅳ)A的全部特征值都有正的實部。

關(guān)于矩陣方程(1)有如下定理。

定理1(Ⅰ)設(shè)A和B都是非奇異M-矩陣，且C≥0，則方程(1)存在唯一非負解。

(Ⅱ)設(shè)A和B是不可約M-矩陣且至少一個是非奇異矩陣，C≥0，則方程(1)存在唯一非負解。

證明 由A和B是非奇異M-矩陣或是不可約M-矩陣且至少一個是非奇異矩陣，可得

P=Im?A+BT?In

是非奇異M-矩陣，因此方程(1)有唯一解。

由引理2中(Ⅱ)得P-1≥0，注意到C≥0，所以解是非負的。

定義2[8]若A∈Rn×n，分裂p層記為Γp,則：

(Ⅰ)對于k∈N,稱A為Γp,k)-HODLR矩陣當且僅當每一個非對角塊A(Ili,Ilj)和A(Ilj,Ili)的秩不超過k，其中i,j={2,1,4,3,…,(2l,2l-1)}，l=1,2,…,p。

(Ⅱ)對于Γp，矩陣A的HODLR秩指：最小正整數(shù)k使得A是Γp,k)-HODLR矩陣。

2 低秩更新

本節(jié)主要考慮矩陣C是非負低秩時，MSE的求解方法。對于系數(shù)矩陣是小型矩陣或稠密矩陣情形，已有多種求解算法，例如矩陣分解法、Hoskins迭代法、交替方向法等[2]。文獻[4]中提出了Smith方法，但Smith方法中的參數(shù)不易選擇，文獻[9]提出了求解MSE方程的交替方向Smith方法(Alternating directional Smith method，ADSM)，迭代序列是非負單調(diào)上升的，且平方收斂。ADSM是Smith方法的改進。具體算法如下：

算法1[9] :ADSM(A,B,C) % 求解AX+XB=C1.令α=max1≤i≤naii, β=max1≤i≤mbii。2.計算 E0=(B-βI)(B+αI)-1, F0=(A+βI)-1(A-αI), X0=(α+β)(A+βI)-1C(B+αI)-1。3.對k=1,2,…,重復(fù)下列步驟直到矩陣序列Xk收斂 Xk+1=Xk+FkXkEk, Ek+1=E2k,Fk+1=F2k。

對于方程(5)，給定如下形式的低秩分解

則

式中U,V均為低秩且令

s=rank(δA)+rank(δB)+rank(δC)。

則方程(5)可以表示為:

AδX+δXB=UVT。

(6)

式中：A=A0+δA,B=B0+δB,U∈Rn×s,V∈Rm×s,s?min{n,m}。低秩近似解Xk=Zk(Yk)T可以直接利用低秩ADI方法求得[5，10]，但迭代中的參數(shù)需要事先給定，且最優(yōu)參數(shù)不易獲得。本文將使用文獻[8]中的擴展Krylov子空間方法，此類算法不需要事先選定參數(shù)?；舅悸肥抢脡KArnoldi過程得到擴展Krylov子空間：

Ut:=εKt(A,U)=

{U,A-1U,AU,A-2U,…,At-1U,A-tU},

Vt:=εKt(BT,V)=

式中2ts

具體算法如下：

算法2 [8]:LRMSE(A,B,U,V) %求解AδX+δXB=UVT1.U1,- =qr(U,A-1 ),V1,- =qr(V,B-1 )。2.對k=1,2,…, 計算 A=UTkAUk,B=VTkBVk,C=UTkUVTVk, Xk←利用稠密方法計算AXk+XkB=C的解, 判斷Xk是否收斂,如果收斂,令 X=UTkXkVk,停止。 否則,令 Uk=U0,U+,U- , Vk=[V0,V+,V-], 計算 U=[AU+,A-1U-], V=[BV+,B-1V-], U←U-UkUTkU, V←U-VkVTkV, U,- =qrU , [V,-]=qr(V), Uk+1=[Uk,U], Vk+1=[Vk,V]。

3 分而治之法

設(shè)方程(1)的系數(shù)矩陣具有分級非對角低秩矩陣結(jié)構(gòu)，首先進行如下分裂：

A=A0+δA，B=B0+δB，C=C0+δC。

(7)

式中:A0,B0,C0是塊對角矩陣；δA,δB,δC是低秩矩陣；此時方程(1)可化為方程(3)和(5)。當對角塊矩陣仍具HODLR結(jié)構(gòu)時，將遞歸地進行分裂。具體遞歸分裂的層數(shù)將根據(jù)矩陣的HODLR結(jié)構(gòu)、計算機的存儲能力和計算速度確定，相關(guān)內(nèi)容見文獻[11]。圖1 展示了HODLR結(jié)構(gòu),對角塊黑色部分為低階矩陣,其余為低秩矩陣。具體如下：

圖1 HODLR形式[8]Fig.1 HODLR form [8]

式中:A12,A21為低秩矩陣，A11,A22仍具有HODLR結(jié)構(gòu)特征，對主對角塊矩陣繼續(xù)分裂，直至最小塊階數(shù)滿足要求后停止。

設(shè)方程(1)的系數(shù)矩陣A,B,C都是HODLR矩陣，則類似文獻[8]，得出MSE的分而治之方法。對系數(shù)矩陣進行如下分裂：

(8)

式中：

類似地分裂B和C：

(9)

(10)

注3分塊后的子矩陣需要保證對應(yīng)的矩陣方程有意義。若對角塊Aii,Bii,Cii,i=1,2仍是HODLR矩陣，則繼續(xù)進行分裂。低秩部分對應(yīng)的更新MSE的近似解將用低秩算法求解，如Krylov子空間方法、 低秩Smith法、低秩ADSM法等。若每一遞歸層均有唯一的非負解，則方程(1)有唯一的非負解。

下列定理可以保證解存在唯一性。

定理2設(shè)A和B是非奇異M-矩陣或是不可約M-矩陣且至少一個是非奇異矩陣，C≥0，矩陣分裂由方程(8)、(9)和(10)得到，則下列結(jié)論成立：

(Ⅰ)矩陣方程

AiiXii+XiiBii=Cii,i=1,2

(11)

是M-矩陣Sylvester方程，且有唯一非負解X11,X22，即X0=diag(X11,X22)≥0是式(3)的非負解；

(Ⅱ)此時得到的矩陣方程(5)是M-矩陣Sylvester方程，且有唯一非負解δX；

(Ⅲ)方程(1)的非負解為X=X0+δX；

(Ⅳ)分裂k次時解為

其中，

證明 (Ⅰ)由引理1中(Ⅴ)可知，A11,B11,A22,B22都是非奇異M-矩陣。由分裂方程(10)可知C11,C22都是非負矩陣，所以方程(11)是MSE，故有唯一非負解。

故公式(5)是M-矩陣Sylvester方程，且有唯一非負解δX。

(Ⅲ)由(Ⅰ)、(Ⅱ)證明可知(Ⅲ)成立。

(Ⅳ)由(Ⅰ)～(Ⅲ)知，當k=1時成立。

k=2時，原方程的解為

設(shè)第k-1次分裂時，原方程的解也滿足

定理3在定理2條件下，且A,B,C具有HOLDR結(jié)構(gòu)，則遞歸層對應(yīng)的矩陣方程都是M-矩陣Sylvester方程。

證明 由引理1中(Ⅱ)可知，A和B的對角塊子矩陣都是非奇異M-矩陣。從而A和B的非對角塊子矩陣都是非正矩陣，C≥0，由定理2可得結(jié)論成立。

下面給出求解MSE的分而治之方法。

算法3 DACMSE(A,B,C) % 求解AX+XB=C1.給定A,B,C。2.當A,B,C具有HOLDR結(jié)構(gòu)時,分塊。3.A=A1100A22+UAVTA,B=B1100B22+UBVTB, C=C1100C22+UCVTC。4.Xii←DACMSE(Aii,Bii,Cii),i=1,2。5.令X0=X11X22。6.令U=[UC,-UA,-X0UB],V=[VC,XT0VA,VB]。7.利用低秩算法求解δX←LRMSE(A,B,U,V)。8.返回:X0+δX。

4 數(shù)值算例

在本節(jié)中，我們通過數(shù)值算例來比較兩種不同分而治之法求解MSE的效果?；诠ぞ甙黨m-toolbox[11]在MATLAB2019b運行所有算法。相對誤差如下:

例1 系數(shù)矩陣分別為：

我們設(shè)置最小塊階數(shù)不超過100，分別使用拓展Krylov分而治之方法(dac-ek)和拓展Krylov ADSM分而治之方法(dac-ek ADSM)求解，從表1 可知dac-ek ADSM在時間、精度方面均得到改進，隨著階數(shù)的增大，計算時間節(jié)省。圖2是求解不同階數(shù)MSE的計算時間。

表1 求解不同階數(shù)下MSE的計算時間與相對誤差Table 1 The times (in seconds) and relative residuals for the algorithms applied to the different orders matrices equations

圖2 求解不同階數(shù)MSE的計算時間Fig.2 The times (in seconds) for the algorithms applied to the different orders matrices equations

例2 系數(shù)矩陣為隨機帶狀M-矩陣，

x=randn-1,1,

A=2eye(n)-diag(x,-1)-diag(x,1),

B=2AT，

C0=diag(rand(n,1))+diag(x,-1)+diag(x,1),

UC=rand(n,2),VC=eye(n,2),

圖3表明拓展Krylov ADSM分而治之法比拓展Krylov分而治之法計算速度更快。

圖3 求解不同階數(shù)MSE的計算時間Fig. 3 The times (in seconds) for the algorithms applied to the different orders matrices equations

5 結(jié)語

本文給出求解大型M-矩陣Sylvester方程化為低階方程和低秩方程分別求解的分而治之算法，證明了遞歸層對應(yīng)的低階方程和低秩矩陣方程都是M-矩陣Sylvester方程。數(shù)值實驗表明本文提出的算法對大型M-矩陣Sylvester方程有效。