修正共軛梯度算法求解四元數(shù)Sylvester張量方程

胡晶晶 ,柯藝芬 ,馬昌鳳

(1.福建師范大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院,福建 福州 350117;2.福建師范大學(xué)分析數(shù)學(xué)及應(yīng)用教育部重點實驗室,福建 福州 350117;3.福建省應(yīng)用數(shù)學(xué)中心,福建 福州 350117)

1.引言

四元數(shù)經(jīng)常出現(xiàn)在許多領(lǐng)域,如量子物理、信號處理、彩色圖像處理和計算機科學(xué)等[1-4].由于四元數(shù)域比復(fù)數(shù)域提供了更多的自由度,因此在處理多維問題時,四元數(shù)比實數(shù)和復(fù)數(shù)具有更多的通用性和靈活性.因此,越來越多的學(xué)者對四元數(shù)的理論性質(zhì)和計算問題感興趣,并取得了許多有價值的成果.例如,文[5-6]通過復(fù)表示方法獲得了四元數(shù)矩陣方程AXB+CY D=E和AHXA+BHY B=C的極小范數(shù)η-(反)-Hermitian最小二乘解的顯示形式.文[7]通過實表示方法推導(dǎo)出了四元數(shù)矩陣方程AXB+CY D=E的η-Hermitian和η-反-Hermitian解.文[8]建立了一種有效的迭代方法去獲得四元數(shù)矩陣方程AXB+CY D=E的極小范數(shù)η-Hermitian和η-反-Hermitian最小二乘解.文[9]考慮了四元數(shù)Sylvester張量方程的最小二乘問題,通過求解其等價形式從而得到了四元數(shù)Sylvester張量方程的最小二乘解.

本文采用以下符號: R表示實數(shù)集,Rn表示n維實向量空間,Rn×n表示n×n階實矩陣集合,表示N階I1×I2×···×IN維實張量的集合,Q表示四元數(shù)體,Qn×n表示n×n階四元數(shù)矩陣集合,表示N階I1×I2×···×IN維四元數(shù)張量集合;AT表示矩陣A的轉(zhuǎn)置,R(A)表示矩陣A的列空間.

2.預(yù)備知識

定義2.1四元數(shù)a可以表示為

其中as ∈R,s=1,2,3,4,且i,j,k滿足

對于兩個四元數(shù)a=a1+a2i+a3j+a4k∈Q和b=b1+b2i+b3j+b4k∈Q,它們的乘積定義為

定義2.2四元數(shù)張量A可以表示為

其中As ∈,s=1,2,3,4.四元數(shù)張量A的共軛可被定義為

類似地,對于四元數(shù)矩陣A,有著相似的表示.

例2.1[10]張量A ∈R3×4×2的兩個正面切片分別為

張量A沿模-1展開成矩陣A(1),其中

定義2.4[11]對于張量A ∈,算子‘vec’將張量的列疊加成一個向量,定義vec(A)=vec(A(1)),其中A(1)是張量A沿模-1的展開矩陣.

例2.2若張量A ∈R3×4×2如例2.1所示,則vec(A)=[1,2,···,24]T.

證四元數(shù)Sylvester張量方程(2.2)可寫作

將式(2.7)和式(2.8)代入式(2.6),根據(jù)定義2.7,四元數(shù)張量相等意味著實部和虛部分別對應(yīng)相等,從而可得方程組(2.3).

進(jìn)一步,根據(jù)Kronecker積和向量化算子,從而可得線性方程組(2.4).

注2.1值得注意的是當(dāng)線性方程組(2.4)的系數(shù)矩陣規(guī)模非常大時,尋找其數(shù)值解將面臨很大的計算挑戰(zhàn).文[12]表明,基于張量格式的算法總體上比經(jīng)典形式的算法效率更高.下面我們將構(gòu)造一類基于張量形式的高效迭代算法求解方程組(2.3),并建立迭代算法的收斂結(jié)果.

注2.2記

3.全局收斂性

本節(jié)主要考慮四元數(shù)張量方程(2.2)的求解.利用定理2.1,可以將四元數(shù)張量方程轉(zhuǎn)化為實張量方程組(2.3),考慮用張量形式的修正共軛梯度(Modify Conjugate Gradient,MCG)算法求解(2.3),從而得到方程(2.2)的解.

為了方便起見,我們引入下列線性算子:

從而四元數(shù)Sylvester張量方程(2.2)可寫成

步4 置k:=k+1,返回步2.

定理3.1假設(shè)序列{Rs(k)}和{Qt(k)}(s=1,2,3,4,t=1,2,3,4)由算法3.1生成,則下列等式成立

證首先,我們證明式(3.3)和式(3.4)對于0≤u

從而,當(dāng)k=1時,式(3.3)和式(3.4)成立.

現(xiàn)假設(shè)式(3.3)和式(3.4)對于k=l成立.對于k=l+1,根據(jù)算法3.1和式(3.2),我們有

最后,我們證明當(dāng)u

從而,當(dāng)k=l+1時,式(3.3)和式(3.4)也成立,則對所有的0≤uv時,式(3.3)和式(3.4)也成立.這樣對所有的uv,結(jié)論得證.

定理3.2記(t=1,2,3,4)為張量方程組(3.1)的解,則對于任給的初始張量Xt(0)(t=1,2,3,4),由算法3.1生成的序列{Rs(k)}和{Qt(k)}(s=1,2,3,4,t=1,2,3,4)滿足

證用數(shù)學(xué)歸納法.當(dāng)k=0時,有

假設(shè)當(dāng)k=l(l ≥1)時,式(3.5)成立.對于k=l+1,我們可得

從而當(dāng)k=l+1時,式(3.5)成立.由數(shù)學(xué)歸納法可得,式(3.5)對所有的k=0,1,2,···都成立.□

定理3.3若實張量方程組(3.1)是相容的,對任意的初始張量[X1(0),X2(0),X3(0),X4(0)],在不計舍入誤差的情況下,其精確解可通過算法3.1在有限迭代步內(nèi)獲得.

通過算法3.1可得到Xt(4I)和Rs(4I),令R(k)=diag{R1(k)(1),R2(k)(1),R3(k)(1),R4(k)(1)},k=0,1,2,···,4I-1,其中Rs(k)(1)表示張量Rs(k)(s=1,2,3,4)沿模-1的展開矩陣,根據(jù)定理3.1,有

則R(0),R(1),···,R(4I-1)可構(gòu)成如下線性空間的一組正交基,

根據(jù)定理3.1,有

故Rs(4I)=O.則[X1(4I),X2(4I),X3(4I),X4(4I)]是實張量方程組(3.1)的解.

4.極小范數(shù)解

引理4.1[13]若A ∈Rm×n,b ∈Rm,線性方程組Ax=b有解x?且x?∈R(AT),則x?是Ax=b的唯一極小范數(shù)解.

定理4.1若實張量方程組(3.1)是相容的,如果選擇初始張量為

由于線性方程組(2.4)和實張量方程組(3.1)等價,根據(jù)引理4.1可得x?為(2.4)的唯一極小范數(shù)解,從而由算法3.1生成的解(t=1,2,3,4)是實張量方程組(3.1)的唯一極小Frobenius范數(shù)解.

5.數(shù)值實驗

為了驗證本文算法的可行性,本節(jié)通過三個例子來說明所提出算法的效果.所有的程序都是基于Bader和Kolda[14]開發(fā)的MATLAB張量工具箱,實驗中,選取N=3,當(dāng)滿足

迭代終止.所有算法均在帶有MATLAB R2018a的個人PC機上運行.此外,定義殘差范數(shù)

考慮初始張量Xt(0)=O(t=1,2,3,4),當(dāng)滿足式(5.1)時迭代終止.根據(jù)算法3.1,我們可獲得式(5.2)的極小Frobenius范數(shù)解,數(shù)值結(jié)果如圖1所示.

圖1 例5.1的數(shù)值結(jié)果

例5.2考慮四元數(shù)Sylvester張量方程(5.2),設(shè)

令精確解

從而可得式(5.2)的右端項張量C.在這個例子中,張量方程(5.2)是相容的,通過選擇初始張量Xt(0)=O(t=1,2,3,4),對于n=10時,執(zhí)行算法3.1得到的收斂曲線如圖2所示.

圖2 例5.2的數(shù)值結(jié)果

下面的例子通過求解稀疏四元數(shù)Sylvester張量方程來檢測算法3.1的有效性.

例5.3考慮四元數(shù)Sylvester張量方程

其中系數(shù)矩陣滿足

與Ahmadi-Asl和Beik[8]工作稍有不同,我們考慮以下四種情況:

情況1A11=S1,A12=pivtol,A21=S2,A13=S3,A14=S4;

情況2A11=S1,A12=gre-343,A21=S2,A13=S3,A14=S4;

情況3A11=S1,A12=mesh3em5,A21=S2,A13=S3,A14=S4;

情況4A11=S1,A12=sphere3,A21=S2,A13=S3,A14=S4.

測試矩陣來自于Davis收集的矩陣[15].數(shù)據(jù)來源于HB(Harwell-Boeing) group.在這個例子中出現(xiàn)的這四種測試矩陣(pivtol,gre-343,mesh3em5和sphere3)經(jīng)常出現(xiàn)在統(tǒng)計、有向加權(quán)圖和結(jié)構(gòu)問題中.這些矩陣的稀疏結(jié)構(gòu)如圖5.1-5.4所示,測試矩陣的性質(zhì)見表5.1,其中‘nz’表示非零項的個數(shù).需要注意的是,這里的稀疏矩陣的階數(shù)較大,因此四元數(shù)Sylvester張量方程(5.3)的解的維數(shù)非常大.令

表5.1 例5.3測試矩陣的性質(zhì)

圖5.1 第一個測試矩陣的稀疏結(jié)構(gòu)

圖5.2 第二個測試矩陣的稀疏結(jié)構(gòu)

圖5.3 第三個測試矩陣的稀疏結(jié)構(gòu)

圖5.4 第四個測試矩陣的稀疏結(jié)構(gòu)

其中對于情況1-4,n分別取為102,343,289和258作為其精確解.

通過選擇初始張量Xt(0)=O,應(yīng)用算法3.1去求解張量方程(5.3),當(dāng)滿足式(5.1)時迭代終止,相應(yīng)的收斂結(jié)果如圖5.5所示,結(jié)果證實了算法3.1的收斂性,可以看出殘差范數(shù)是顯著下降的.因此,算法3.1對于求解稀疏張量方程是有效的.

圖5.5 例5.3的收斂結(jié)果

6.結(jié)論

本文通過將四元數(shù)Sylvester張量方程等價地轉(zhuǎn)化為實數(shù)域上的張量方程組并引入線性算子,構(gòu)造了基于張量形式的修正共軛梯度算法求其等價形式.理論分析表明對任給的初始張量,由算法3.1生成的迭代序列在不計舍入誤差的情況下,可在有限迭代步內(nèi)收斂到方程組的精確解.進(jìn)一步,通過選擇特殊類型的初始張量,可獲得唯一極小Frobenius范數(shù)解.最后,數(shù)值例子證明了該算法的有效性.