陳世軍， 余勝斌

(陽光學院 基礎教研部，福建 福州 350015)

在工程學、控制理論等領(lǐng)域中很多問題都可以歸結(jié)為求矩陣方程約束解.近年來，諸多學者對單變量二次矩陣方程進行了許多研究，取得了豐富的成果并建立了一些有效的算法[1—4].例如，崔學蓮[1]提出了二次矩陣方程X2+AX-B=O在矩陣區(qū)間[L,U]的約束解的迭代算法，并證明了該算法的收斂性；金繼承[2]提出了保結(jié)構(gòu)加倍算法,用于阻尼系統(tǒng)的二次方程的數(shù)值求解，并給出了算法的迭代格式和收斂性證明；李媛媛[4]利用固定點理論和矩陣范數(shù)的相關(guān)知識, 給出了二次矩陣方程X2+2AX+B=O有解的充要條件和有唯一解的判定定理.共軛梯度算法是求矩陣方程約束解的有效算法之一.近幾十年來，很多學者研究了基于共軛梯度算法的矩陣方程(組)的約束解.對于非線性矩陣方程求解，牛頓算法和共軛梯度算法是一種有效的算法[5—17].孫穎異[5]通過修正共軛梯度參數(shù)，構(gòu)造新的搜索方向, 提出兩類修正的WYL共軛梯度法，證明了算法的全局收斂性； J.H.Long[6]等建立了求解二次矩陣方程的改進牛頓算法，且給出了算法的收斂性證明；張凱院[10]建立了雙矩陣變量Riccati矩陣方程AX1B+CX2D+X1E1X2+X2E2X1=E3對稱解的牛頓-MCG算法.對于多變量二次矩陣方程異類約束解，目前還沒有得到研究，因此本文基于牛頓算法和修正共軛梯度算法，建立了求多變量二次矩陣方程的異類約束解，這里的約束解是指對稱解、中心對稱解和自反解.

本文研究三變量二次矩陣方程的異類約束解，矩陣方程形式如下：

M1X1N1+M2X2N2+M3X3N3+

X1H1X2+X2H2X3+X1H3X3+Q=O.

(1)

1 求矩陣方程(1)異類約束解的牛頓算法

定義1劃分n階單位矩陣I=(e1,e2,…,en)，稱S=(en,en-1,…,e1)為n階次單位矩陣．若矩陣X∈Rn×n滿足SXS=X，則稱X為中心對稱矩陣.

定義2設P∈Rn×n滿足PT=P，P2=I，若X∈Rn×n滿足PXP=X，則稱X為關(guān)于矩陣P的自反矩陣.

為了使表達式更簡潔，引進以下矩陣記號：

對矩陣函數(shù)ψ(X+tY)在t=0處求導，得

M1Y1N1+M2Y2N2+M3Y3N3+

X1H1Y2+Y1H1X2+X2H2Y3+

Y2H2X3+X1H3Y3+Y1H3X3.

記

φX(Y)=φX(Y1,Y2,Y3)=

M1Y1N1+M2Y2N2+M3Y3N3+

X1H1Y2+Y1H1X2+X2H2Y3+

Y2H2X3+X1H3Y3+Y1H3X3.

γ(Y)=Y1H1Y2+Y2H2Y3+Y1H3Y3.

因此可得

ψ(X+Y)=ψ(X1+Y1,X2+Y2,X3+Y3)=

ψ(X)+φX(Y)+γ(Y)，

(2)

其中φX(Y)是指矩陣函數(shù)ψ(X)在“點X”處沿著“方向Y”的Frechet導數(shù).

不妨假設(X1,X2,X3)∈Ω1-3-7是二次矩陣方程(1)的一組近似解，則求二次矩陣方程(1)的一組解(X1,X2,X3)等同于求一組校正值(Y1,Y2,Y3)∈Ω1-3-7使之滿足方程ψ(X+Y)=O，當(Y1,Y2,Y3)譜半徑均小于1時，由牛頓算法原理可以忽略矩陣(Y1,Y2,Y3)的高次項γ(Y)，因而方程(2)可以近似為ψ(X+Y)≈ψ(X)+φX(Y).因此求解ψ(X+Y)=O可轉(zhuǎn)化為求(Y1,Y2,Y3)∈Ω1-3-7使之滿足線性矩陣方程

φX(Y)=-ψ(X).

(3)

結(jié)合牛頓算法原理，建立求矩陣方程(1)異類約束1-3-7解的牛頓算法：

φX(k)(Y(k))=-ψ(X(k)).

Step 3 計算X(k+1)=X(k)+Y(k)，令k∶=k+1，轉(zhuǎn)到Step 2.

牛頓算法中有收斂性結(jié)論：假設矩陣X*是二次矩陣方程(1)的一個單根，并且初始矩陣X(1)能充分接近X*，則牛頓算法迭代所得矩陣序列{X(k)}能二次收斂于X*.

2 問題的提出

記A1=M1,B1=N1,C1=M2,D1=N2,E1=M3,F1=N3，A2=I,B2=H1X2,C2=X1H1,D2=I,E2=X2H2,F2=I，A3=I,B3=H3X3,C3=I,D3=H2X3,E3=X1H3,F3=I,G=-ψ(X).建立求線性矩陣方程(3)異類約束解的MCG算法，給出方程(3)的一般形式：

(4)

其中Ai,Bi,Ci,Di,Ei,Fi,G∈Rn×n(i=1,2,3).主要討論下面2個問題：

問題1假設關(guān)于Y1,Y2,Y3的方程(4)有解，則求(Y1,Y2,Y3)∈Ω1-3-7，使方程(4)成立.

問題2假如方程(4)無異類約束1-3-7解，則求(Y1,Y2,Y3)∈Ω1-3-7，使得

3 求解問題1的MCG算法

基于共軛梯度算法的基本思路，結(jié)合方程(4)中矩陣特點，修改有關(guān)矩陣的類型，建立求解問題1的MCG算法.先引入記號：

算法1 (問題1的MCG算法)

Step 2 若Rk=O，或者Rk≠O且Qk=O，則停止計算；否則，繼續(xù)計算

Step 3 計算

Rk+1=G-μ(Y(k+1))，

Qk+1=

Step 4 令k∶=k+1，轉(zhuǎn)Step 2計算.

性質(zhì)2若k≥2，則算法1中得到的矩陣Ri和Qj滿足

i≠j；i,j=1,2,…,k.

4 求解問題2的MCG算法

由定理1可知當Rk≠O而Qk=O時算法1中斷，即方程(4)無異類約束1-3-7解.因此，需要求方程(4)的最小二乘異類約束解，即求矩陣(Y1,Y2,Y3)∈Ω1-3-7使得

下面通過將矩陣方程組按行拉直，構(gòu)造方程(4)的等價正規(guī)矩陣方程，將求解問題2轉(zhuǎn)化為求解等價矩陣方程異類約束1-3-7解的問題.然后參照算法1，建立求解問題2的MCG算法.

引進記號：

u(Y1,Y2,Y3)=

v(Y1,Y2,Y3)=

M=

定理2求解問題2等價于求約束正規(guī)矩陣方程f(Y)=K的異類約束1-3-7解，且該正規(guī)矩陣方程一定有異類約束1-3-7解.

證明求解問題2等價于求(Y1,Y2,Y3)∈Ω1-3-7，使得

(5)

下面證明求解極小值問題(5)等價于求矩陣方程f(Y)=K的異類約束解.將矩陣方程組

(6)

參照算法1以及MCG算法原理，建立求矩陣方程f(Y)=K異類約束解，即求解問題2的MCG算法如下：

算法2(問題2的MCG算法)

Step 2 若Rk=O，則停止算法2計算；否則，繼續(xù)計算

Step 3 計算

Step 4 令k∶=k+1，轉(zhuǎn)Step 2.

5 數(shù)值算例

下面給出算例1，在算例1中通過改變矩陣的階數(shù)，說明文中建立的MCG算法1和MCG算法2是可行的，具有很高的效率，且能求出矩陣方程(1)異類約束1-3-7解或異類約束最小二乘解.這里計算時間(s)、系數(shù)矩陣的階數(shù)、MCG算法1、MCG算法2、MCG算法1中斷次數(shù)、實際誤差的F-范數(shù)分別用time、n、k1、k2、k12和error表示.結(jié)合牛頓算法以及MCG算法1、MCG算法2，采用如下步驟計算：

Step 3 計算X(k+1)=X(k)+Y(k)，置k∶=k+1，轉(zhuǎn)Step 2.

例1用文中建立的牛頓-MCG雙迭代算法求矩陣方程(1)的一組異類約束1-3-7解和最小二乘解，方程(1)系數(shù)矩陣均為n階方陣，系數(shù)矩陣、終止準則ε如下：

Mi=Ni=Hi=I(i=1,2,3)，

選取初始矩陣X1=1.3I,X2=0.94I,X3=1.05I，Y1=Y2=Y3=O，按計算步驟求得矩陣方程(1)的異類約束1-3-7解.結(jié)果見表1(MATLAB軟件2014版-PIV3.0 GHZ微機).

表1 方程(1)的異類約束解計算結(jié)果

當n=5時，可得矩陣方程(1)的異類約束1-3-7解為

若選取初始矩陣為X1=2I,X2=0.2I,X3=3I，Y1=Y2=Y3=O，n=5，按計算步驟得到矩陣方程(1)的異類約束1-3-7解為

從表1可以看出，本文建立的牛頓-MCG雙迭代算法可以求出矩陣方程(1)的一組異類約束1-3-7解，且有很高的計算效率.當n=5時，取不同的初始矩陣，得到了不同的異類約束1-3-7解，這說明矩陣方程(1)異類約束1-3-7解不唯一，這也是牛頓-MCG算法的特點，不要求方程解唯一，只要求方程(1)有解即可.