基于多項(xiàng)式預(yù)處理的特殊雙變量矩陣方程異類約束解算法

周咸富，段復(fù)建

(桂林電子科技大學(xué) 數(shù)學(xué)與計(jì)算科學(xué)學(xué)院，廣西 桂林 541004)

討論一類特殊的雙變量矩陣方程

AXB+AYB=F

(1)

的異類約束解加速問題，其中A，B，F(xiàn)∈Rn×n，X∈Rn×n，Y∈SRn×n，且A、B是對(duì)稱正定矩陣。對(duì)于方程(1)的研究,不同約束條件可用不同的方法，一般采用共軛梯度法求解。劉莉[1]利用共軛梯度法求矩陣方程AXB+CYD=E的中心對(duì)稱最小二乘解；方玲等[2]利用極小殘差法求解AXB+CYD=E的對(duì)稱最小二乘解。但這些方法迭代過程中的收斂性不確定，即便收斂，收斂速度也較慢。針對(duì)收斂速度慢的問題，進(jìn)行預(yù)處理是提高收斂速度的主要途徑之一。Zhao等[3]提出了一種基于參數(shù)迭代預(yù)處理方法和共軛梯度修正形式(MCG)相結(jié)合的并行預(yù)條件求解算法；溫瑞萍等[4]提出了不完全SAOR預(yù)條件共軛梯度法；Jia等[5]提出了一種基于殘差的稀疏近似反預(yù)處理方法。對(duì)于矩陣方程異類約束解問題，李書連等[6]通過修改系數(shù)矩陣和其方向，提出了修正共軛梯度算法，解決了雙變量矩陣方程組最小二乘問題；牛婷婷等[7]將多變量異類約束問題推廣到一類離散時(shí)間代數(shù)Riccati矩陣方程，建立了雙迭代算法；宋衛(wèi)紅等[8]利用逆矩陣的Neumann級(jí)數(shù)形式，將多變量異類約束解推廣到離散對(duì)偶代數(shù)Riccati矩陣方程，構(gòu)造了雙迭代算法。

預(yù)處理方法大多應(yīng)用于解同類線性方程組，用于雙變量矩陣方程異類約束條件的研究不多，且對(duì)于異類約束的矩陣方程收斂速度的研究也不多見。因此，將多項(xiàng)式預(yù)處理[9]和共軛梯度法相結(jié)合，提出多項(xiàng)式預(yù)處理共軛梯度法。

1 多項(xiàng)式預(yù)處理技術(shù)

定義同階對(duì)稱正定矩陣A與B的內(nèi)積為[A,B]=tr(ATB)，由此導(dǎo)出的Frobenius范數(shù)為

設(shè)Rn×n為n×n階矩陣集合，SRn×n為n×n階對(duì)稱矩陣集合，σ為A的任意奇異值，S=[a,b]為σ的取值區(qū)間，a和b分別為σ的最小與最大奇異值,τ為B的任意奇異值，L=[d,e]為τ的取值區(qū)間，d和e分別為λ的最小、最大奇異值。

迭代法求解矩陣方程，收斂速度與系數(shù)矩陣奇異值的分布有關(guān)。奇異值分布越集中，條件數(shù)越接近于1，收斂速度越快，反之，則會(huì)很慢。利用多項(xiàng)式對(duì)矩陣方程進(jìn)行預(yù)處理，可使其奇異值的分布相對(duì)集中，條件數(shù)接近于1，從而提高收斂速度。

對(duì)矩陣方程AXB+AYB=F作等價(jià)變換

C(A)AXBC(B)+C(A)AYBC(B)=C(A)FC(B)，

(2)

其中C(A)、C(B)為次數(shù)小于n的實(shí)系數(shù)多項(xiàng)式矩陣。若C(A)A和C(B)B的譜條件遠(yuǎn)遠(yuǎn)小于A、B的譜條件數(shù)，則利用共軛梯度法求解式(2)的收斂速度將大大快于直接求式(1)的收斂速度。所以選取系數(shù)多項(xiàng)式時(shí)，應(yīng)使C(A)和C(B)在某種意義下接近A-1和B-1，這樣經(jīng)過處理后的迭代次數(shù)遠(yuǎn)比之前少。若A或B的最大奇異值和最小奇異值相差不大，則不需要進(jìn)行預(yù)處理。令

F(A)=C(A)A,G(B)=BC(B),

Q=C(A)FC(B)，

則式(1)的求解問題等價(jià)于

F(A)XG(B)+F(A)YG(B)=Q

(3)

的求解。

由多項(xiàng)式插值原理可得A的預(yù)處理多項(xiàng)式：

(4)

(5)

(6)

a≤σi≤b，i=1，2，…，n，

d≤τi≤e，i=1,2,…,n，

定義1設(shè)A∈Rn×n，A為非奇異矩陣，‖·‖為定義在Rn×n上的范數(shù)矩陣，稱

cond(A)v=‖A-1‖v‖A‖v，v=1,2,or,∞

為A關(guān)于范數(shù)‖·‖的條件數(shù)。

通常情況下，若矩陣A的條件數(shù)大，則稱A為病態(tài)矩陣，反之稱為良態(tài)矩陣，條件數(shù)越趨近于1，則越穩(wěn)定。

采用最常見的Frobenius范數(shù)反映矩陣病態(tài)性，

當(dāng)A為實(shí)對(duì)稱矩陣時(shí)，有

其中λmax、λmin分別為矩陣A最大、最小特征值。

定理1設(shè)式(1)中系數(shù)的譜條件數(shù)較大，經(jīng)過一次多項(xiàng)式預(yù)處理后，得到新的矩陣方程最小與最大奇異值比值比原方程最小與最大奇異值比值提高16倍。

證明設(shè)經(jīng)過一次插值多項(xiàng)式預(yù)條件，得到C(A)A的最小、最大特征值分別為a′、b′，C(B)B的最小、最大特征值分別為d′、e′。記

2 多項(xiàng)式預(yù)處理共軛梯度算法

對(duì)于方程(1)用共軛梯度法求解，若條件數(shù)遠(yuǎn)大于1，則收斂速度較慢。為了提高收斂速度，引入多項(xiàng)式預(yù)處理技術(shù)構(gòu)造新的算法。

算法1定義1+ε1,1+ε2,(ε1>0,ε2>0)分別為系數(shù)矩陣A、B的非零最大、最小奇異值比值的滿意的臨界值。

1)給出初始解X0∈Rn×n、Y0∈SRn×n和初始邊界a0、b0、d0、e0。

2)對(duì)于i=0,1,…,n，求矩陣多項(xiàng)式：

3)若bi/ai≤1+ε,ei/di≤1+ε2，則停止預(yù)處理，轉(zhuǎn)步驟5)，否則，令

Ai+1=F(Ai)，

Bi+1=G(Bi)，

Qi+1=C(Ai)FC(Bi)。

4)計(jì)算ai+1=1，bi+1=(ai+bi)2/4aibi，di+1=1，ei+1=(di+ei)2/4diei，轉(zhuǎn)步驟2)，i=i+1。

5)令

F(A)=F(Ai)，G(B)=G(Bi)，Q=Qi+1。

6)計(jì)算

R0=Q-(F(A)X0G(B)+F(A)Y0G(B))，

P0=FT(A)R0GT(B)，

U0=P0，

V0=Q0，k=0。

7)若Rk=0，或Rk≠0，而‖Uk‖2+‖Vk‖2=0，停止，否則，k=k+1。

8)計(jì)算

Xk=Xk-1+αk-1Uk-1，

Yk=Yk-1+αk-1Vk-1，

Rk=Q-(F(A)XkG(B)+F(A)YkG(B))。

9)計(jì)算

10)轉(zhuǎn)步驟8)。

需要注意的是，在實(shí)際求解中，系數(shù)矩陣的最大與最小奇異值很難求出，因此，利用上述算法在求解C(A)、C(B)時(shí)，對(duì)于a0、b0、d0、e0可利用合適的估計(jì)式，

對(duì)于算法1，一般都有以下2個(gè)性質(zhì)。

性質(zhì)1若矩陣方程(3)是相容的，算法中{Ri}、{Ui}、{Vi}滿足以下關(guān)系：

〈Ri,Rj〉=0, 〈Ui,Uj〉+〈Vi,Vj〉=0，

其中i,j=1，2，…，k，i≠j。

性質(zhì)2若線性方程(1)是相容的，且[X*,Y*]是其一個(gè)解，則對(duì)于任意初始矩陣對(duì)[X0,Y0]，算法中的迭代序列{Xi}、{Yi}、{Ri}、{Pi}、{Qi}、{Ui}、{Vi}滿足如下關(guān)系：

其中i=1,2,…，k-1。

定理2若矩陣方程(1)有解，則對(duì)任意選定的初始迭代矩陣對(duì)[X1,Y1]，應(yīng)用算法1，矩陣方程(1)異類約束解可經(jīng)過有限步迭代得到。

證明若Ri≠0，i=1,2,…,pq，由性質(zhì)2可得‖Ui‖2+‖Vi‖2≠0。根據(jù)算法1，由迭代計(jì)算可得Rpq+1、Upq+1、Vpq+1，且由性質(zhì)1知

〈Rpq+1,Ri〉=0， 〈Ri,Rj〉=0，

其中i,j=1,2,…,pq,i≠j。序列{Ri}構(gòu)成了矩陣空間Rp×q的一個(gè)正交基，故Rpq+1=0，即Xpq+1、Ypq+1為方程的一個(gè)解。

引理1設(shè)算法產(chǎn)生的迭代序列{Xk}在某種范數(shù)意義下收斂,

若存在實(shí)數(shù)α>0和一個(gè)與迭代次數(shù)無關(guān)的常數(shù)0

‖Xk+1-X*‖=q‖Xk-X*‖α,

則稱算法產(chǎn)生的迭代序列{Xk}具有Q-α階收斂速度。特別地，當(dāng)α=1,0

為證明多項(xiàng)式預(yù)處理共軛梯度法產(chǎn)生的迭代序列{Xk}、{Yk}具有Q-線性收斂速度，作如下變換。

記Mk=Xk-X*,Nk=Yk-Y*，則Mk+1=Mk+αkUk,Nk+1=Nk+αkVk,即

且由性質(zhì)2可知，對(duì)任意的i≠k，有

(7)

(8)

結(jié)合算法1，經(jīng)過計(jì)算證明，產(chǎn)生的序列{Mk}、{Nk}具有如下結(jié)論。

性質(zhì)3對(duì)k≥1,有

(9)

(10)

證明由式(7)、(8)可得

tr((Mk+αkUk)T(Mk+αkUk))+

tr((Nk+αkVk)T(Nk+αkVk))=

2αk‖Rk‖2=‖Mk‖2+‖Nk‖2+

(‖Uk‖2+‖Vk‖2)=

定理3設(shè)迭代序列{Xk}、{Yk}由多項(xiàng)式預(yù)處理共軛梯度算法經(jīng)過0次預(yù)處理產(chǎn)生，且A、B為對(duì)稱正定矩陣，則有

(11)

其中X*、Y*是方程(1)的解，σn=a,σ1=b,τq=d,τ1=e分別是A、B的最小、最大奇異值。

證明假定A、B的奇異值分解為

A=UΣVT,B=WDTT。

(12)

其中：U、V、W、T皆為正交矩陣；

且

b=σ1≥σ2≥…≥σn=a，

e=τ1≥τ2≥…≥τq=d。

由式(10)、(12)可得，

‖Mk+1‖2+‖Nk+1‖2+

‖Mk+1‖2+‖Nk+1‖2+

‖Mk+1‖2+‖Nk+1‖2+

‖Mk+1‖2+‖Nk+1‖2+

‖Mk+1‖2+‖Nk+1‖2+

‖Mk+1‖2+‖Nk+1‖2+

‖Mk+1‖2+‖Nk+1‖2≤

即有

定理4多項(xiàng)式預(yù)處理共軛梯度算法經(jīng)過一次預(yù)處理后的Q-收斂速度為

證明由定理1知，經(jīng)過一次預(yù)處理，最小與最大奇異值比值提高了16倍，結(jié)合定理3可得結(jié)論，證畢。

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

已知矩陣A,B,F∈R5×5，且A、B是對(duì)稱正定矩陣，求X∈R5×5,Y∈SR5×5，使

AXB+AYB=F。

X0=zeros(5)， Y0=zeros(5)。

其中A的最大和最小奇異值分別為b=37.436 0,a=0.024 0，B的最大和最小奇異值分別為e=97.310 1,d=0.025 1。

使用共軛梯度迭代法進(jìn)行求解，在迭代過程中，取ε=1×10-8作為迭代終止條件，即當(dāng)‖Rk‖<ε，迭代停止。對(duì)于本例，可以求得如下結(jié)果：

迭代次數(shù)k=3 510，迭代時(shí)間t=0.287 921 s，誤差‖R3 510‖=6.637 3×10-9<ε。

利用多項(xiàng)式預(yù)處理共軛梯度算法，求解AXB+AYB=F。同樣的迭代終止條件，本算法的迭代結(jié)果如表1所示。從表1可看出，預(yù)處理次數(shù)i=10，迭代次數(shù)k=12，迭代時(shí)間t=0.005 191 s，本算法比共軛梯度法迭代次數(shù)減少了3 498，迭代時(shí)間縮短了0.282 730 s。

表1 本算法的迭代結(jié)果

4 結(jié)束語

基于共軛梯度算法，引入多項(xiàng)式預(yù)處理技術(shù)，構(gòu)造了多項(xiàng)式預(yù)處理共軛梯度算法，并用于求解特殊雙變量矩陣方程異類約束的問題，同時(shí)證明了該算法是收斂，且具有Q-線性收斂速度。數(shù)值實(shí)驗(yàn)表明,多項(xiàng)式預(yù)處理共軛梯度算法比共軛梯度算法收斂速度更快，需要的迭代時(shí)間更短，且算法是有效、可行的。