張 潔,唐 嘉,馬昌鳳

(福建師范大學(xué) 數(shù)學(xué)與信息學(xué)院,福建 福州,350117)

矩陣方程經(jīng)常出現(xiàn)在許多應(yīng)用領(lǐng)域,例如控制理論、系統(tǒng)理論、 穩(wěn)定性理論,以及模型還原和圖像處理等。由于矩陣方程廣泛的應(yīng)用背景,使得它吸引了許多研究者的廣泛關(guān)注。求解矩陣方程的方法有很多,例如參數(shù)迭代算法、共軛梯度迭代算法、矩陣分解法等,其中參數(shù)迭代法倍受關(guān)注[1-3]。 在過去的幾十年,線性矩陣方程得到了廣泛的研究[4-14]。

考慮控制理論中提出的矩陣方程

AXBT+BXAT=F

(1)

其中A,B,F∈Rn×n是給定的矩陣,X∈Rn×n是待求的未知矩陣。方程(1)有唯一解X當(dāng)且僅當(dāng)(B?A+A?B)的特征值不為零。若A和B均為可逆矩陣,則矩陣方程(1)等價于

(A-1B)X+(A-1B)T=A-1FA-T

(2)

方程(2)有唯一解X當(dāng)且僅當(dāng)A-1B沒有相反的特征值。特別地,當(dāng)A,B是對稱正定矩陣時,(B?A+A?B)也是對稱正定矩陣,所以矩陣方程(1)有唯一解。

1 參數(shù)迭代法與收斂性分析

當(dāng)矩陣方程(1)有唯一解X*時的參數(shù)迭代法以及最優(yōu),近似最優(yōu)參數(shù)的選取問題[3]。選取實數(shù)α≠0,使αA-1B-I為可逆矩陣,令

C=(αA-1B-I)-1(αA-1B+I)=I+2(αA-1B-I)-1

(3)

那么,C-I=2(αA-1B-I)-1也是可逆矩陣,且有

αA-1B=(C-I)-1(C+I)

(4)

方程(2)兩端乘以α,再將(4)代入可得

可得

X=CXCT+Y0

(5)

因為α≠0且C-I可逆,所以方程(5)與(1)等價。根據(jù)方程(5)建立參數(shù)迭代格式

Xk+1=CXkCT+Y0,k=0,1,2,…

(6)

用拉直算子將式(6)轉(zhuǎn)化為列向量形式

vec(Xk+1)=(C?C)vec(Xk)+vec(Y0),k=0,1,2,…

(7)

式(7)是一階線性定常迭代格式,它的迭代矩陣為C?C,且有ρ(C?C)=ρ2(C)。

有如下的收斂性定理。

定理1[3]對于任意的初始矩陣X(0),格式(6)收斂的充要條件是ρ(C)<1。

當(dāng)A,B是可逆矩陣時,設(shè)C的任一特征值為λc,對應(yīng)的特征向量為x,則有Cx=λcx,即 (αA-1B-I)-1(αA-1B+I)x=λcx,于是α(λc-1)A-1Bx=(λc+1)x。因為A-1B可逆,所以λc≠-1。于是有,

(8)

(9)

(10)

即C的特征值可由B-1A的特征值表示出來。

定理2[3]設(shè)A,B可逆,對任意的初始矩陣X0,格式(6)收斂的充要條件是

(11)

推論1 設(shè)A可逆,且B-1A的特征值的實部大(小)于零,則對任何負(正)數(shù)α,格式(6)收斂。

定理1表明,只要選取適當(dāng)實數(shù)α，使得ρ(C)<1,格式(6)就能夠收斂。需要研究的問題是：如何選取實數(shù)α,才能使格式(6)收斂最快,即使ρ(C)最小。

2 最優(yōu)參數(shù)的選取

證明 因為B-1A=B-1/2(B-1/2AB-1/2)B1/2,即B-1相似于B-1/2AB-1/2,而B-1/2AB-1/2合同于A,所以B-1A的特征值為正數(shù)。根據(jù)定理1,對于任何負數(shù)α,格式(6)收斂。B-1A的特征值排序為η1≥η2≥…ηn>0,由式(10)可得

(12)

(13)

ρ(C)的最小值只能在g1(α)和gn(α)的交點處取得。令g1(α)=gn(α),可求得

(14)

確定格式(6)的最優(yōu)參數(shù)時,需要計算B-1A的最大和最小特征值。為了避免計算B-1,下面討論格式(6)的近似最優(yōu)參數(shù)的選取問題。

證明 設(shè)A和B的特征值分別為

λ1≥λ2≥…≥λn,u1≥u2≥…≥un。

根據(jù)廣義特征值的極值原理,由式(8)可得

(15)

(16)

(17)

易見ρ(C)≤q(α),欲使ρ(C)最小,只要q(α)最小即可。類似于(14)的推導(dǎo),可得

(18)

3 迭代校正法

迭代法的收斂性可以通過迭代校正得到提高?？紤]矩陣方程(1),記g(X)=AXBT+BXAT。此外,設(shè)矩陣方程的唯一解為X*,X1,X2,…,Xm(m>1)為X*的不同近似矩陣,且g(Xi)≠F,(i=1,2,…,m)。利用諸Xi提供了的信息構(gòu)造矩陣X,使?jié)M足

(19)

稱矩陣X為方程(1)關(guān)于X1,X2,…,Xm(m>1)的校正解,確定X的過程為校正過程。整體校正模型(WMC)[4]如下

(20)

所謂整體校正,是指Xi的每一個元素對X的對應(yīng)元素的貢獻比例相同。

將校正過程施加于任何一種求解矩陣方程(1)的迭代格式,可以改善原格式的收斂狀態(tài)。設(shè)求解矩陣方程(1)的單步迭代格式為

Xi=φ(X(i-1)),i=1,2,3,…

(21)

改進格式(21)為

(22)

一方面,當(dāng)格式(21)收斂時,格式(22)必定收斂。進一步,當(dāng)?shù)趇次校正過程的整體縮減系數(shù)αi<1時,格式(22)比(21)收斂的快。另一方面,當(dāng)格式(22)不收斂時,由于αi<1幾乎對所有i成立,所以格式(22)也能收斂。

4 參數(shù)迭代加速算法

下面討論關(guān)于迭代格式(6)的加速方法。

令X0=Y0,則根據(jù)數(shù)學(xué)歸納法,迭代格式(6)可寫為

(23)

當(dāng)ρ(C)<1,則

(24)

其中X*是原方程的精確解。

采用級數(shù)(24)的部分和逼近X*時,為了加速收斂過程,可取該級數(shù)的前2{k+1}項構(gòu)成的部分和進行迭代計算,即

于是,可得部分和迭代格式

(25)

5 數(shù)值實驗

給出兩個數(shù)值例子來驗證對其最優(yōu)參數(shù)選取以及對應(yīng)迭代校正算法和加速算法的有效性。所有的實驗都是在如下環(huán)境中實現(xiàn)的：Windows 10系統(tǒng),Matlab R2015a。

例1 考慮矩陣方程AXBT+BXAT=I,其中

通過計算得到矩陣A,B的特征值

λ1=0.4903,λ2=1.9128,λ3=9.5969
u1=2.1981,u2=3.5550,u3=5.2470

B-1A的特征值為η1=2.0825,η2=0.2029,η3=0.5195。

由(14)和(18),可得αopt=-0.6501,αaopt=-0.6387。

根據(jù)以上的算法,得到矩陣方程的解和停止準則為10-10。對α取不同的值,則數(shù)值結(jié)果如下。從圖1和圖2可以看出參數(shù)迭代法使用最優(yōu)參數(shù)(αopt=-0.6501)和近似最優(yōu)參數(shù)(αaopt=-0.6387)與加速算法的比較。從表1，可以看到參數(shù)迭代法(6)和迭代校正法(22)的比較。 同時,可以從表2看到參數(shù)迭代法(6)和加速算法(25)的比較。在經(jīng)過迭代校正和加速后，易見(6)的收斂性大大地提高了。

圖1 例1中的迭代次數(shù)與誤差(αopt=-0.6501)Fig.1 The iterations and errors for Case 1

圖2 例1中的迭代次數(shù)與誤差(αaopt=-0.6387)Fig.2 The iterations and errors for Case 1

表1 不同α對于例1的數(shù)值結(jié)果

表2 不同α對于例1的數(shù)值結(jié)果

其中

通過計算得到矩陣方程AXBT+BXAT=I的解(m=100)。其中停止準則為‖xk+1-xk‖≤10-10令初始矩陣X0=Y0,對α取不同的值,則數(shù)值結(jié)果如下。根據(jù)定理4,求得αopt=-1.9916和αaopt=-1.8614。從圖3和圖4可以看出參數(shù)迭代法使用最優(yōu)參數(shù)(αopt=-1.9916)和近似最優(yōu)參數(shù)(αaopt=-1.8614)與加速算法的比較。同時,可以從表3 和表4看到參數(shù)迭代法(6)分別與迭代校正法(22)加速算法(25)的比較。數(shù)值結(jié)果見表3和表4。

圖3 例2中的迭代次數(shù)與誤差(αopt=-1.9916)Fig.3 The iterations and errors for Case 2

圖4 例2中的迭代次數(shù)與誤差(αaopt=-1.8614)Fig.4 The iterations and errors for Case 2

表3 不同α對于例2的數(shù)值結(jié)果

表4 不同α對于例2的數(shù)值結(jié)果

續(xù)表