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

胡晶晶 ,柯藝芬 ,馬昌鳳

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

1.引言

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

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

2.預備知識

定義2.1四元數a可以表示為

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

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

定義2.2四元數張量A可以表示為

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

類似地,對于四元數矩陣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.

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

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

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

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

注2.2記

3.全局收斂性

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

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

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

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

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

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

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

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

最后,我們證明當u

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

定理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)滿足

證用數學歸納法.當k=0時,有

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

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

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

通過算法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的展開矩陣,根據定理3.1,有

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

根據定理3.1,有

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

4.極小范數解

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

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

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

5.數值實驗

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

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

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

圖1 例5.1的數值結果

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

令精確解

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

圖2 例5.2的數值結果

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

例5.3考慮四元數Sylvester張量方程

其中系數矩陣滿足

與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].數據來源于HB(Harwell-Boeing) group.在這個例子中出現的這四種測試矩陣(pivtol,gre-343,mesh3em5和sphere3)經常出現在統計、有向加權圖和結構問題中.這些矩陣的稀疏結構如圖5.1-5.4所示,測試矩陣的性質見表5.1,其中‘nz’表示非零項的個數.需要注意的是,這里的稀疏矩陣的階數較大,因此四元數Sylvester張量方程(5.3)的解的維數非常大.令

表5.1 例5.3測試矩陣的性質

圖5.1 第一個測試矩陣的稀疏結構

圖5.2 第二個測試矩陣的稀疏結構

圖5.3 第三個測試矩陣的稀疏結構

圖5.4 第四個測試矩陣的稀疏結構

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

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

圖5.5 例5.3的收斂結果

6.結論

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