修正共轭梯度算法求解四元数Sylvester张量方程.pdf

应用数学MATHEMATICA APPLICATA2024,37(1):238-250修正共轭梯度算法求解四元数Sylvester张量方程胡晶晶1,2,3,柯艺芬1,2,3,马昌凤1,2,3(1.福建师范大学数学与统计学院,福建 福州 350117;2.福建师范大学分析数学及应用教育部重点实验室,福建 福州 350117;3.福建省应用数学中心,福建 福州 350117)摘要：本文提出一类张量形式的修正共轭梯度算法求解四元数Sylvester张量方程.证明在不计舍入误差的情况下,所提方法可在有限迭代步内获得张量方程组的解.进一步,通过选择特殊类型的初始张量,可获得方程组的唯一极小Frobenius范数解.通过数值算例验证了所提出算法的可行性和有效性.关键词：四元数;Sylvester张量方程;修正共轭梯度算法中图分类号：O241.6AMS(2010)主题分类：15A69;15B33;65J10文献标识码：A文章编号：1001-9847(2024)01-0238-131.引言四元数经常出现在许多领域,如量子物理、信号处理、彩色图像处理和计算机科学等14.由于四元数域比复数域提供了更多的自由度,因此在处理多维问题时,四元数比实数和复数具有更多的通用性和灵活性.因此,越来越多的学者对四元数的理论性质和计算问题感兴趣,并取得了许多有价值的成果.例如,文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维实向量空间,Rnn表示n n阶实矩阵集合,RI1I2IN表示N阶I1 I2 IN维实张量的集合,Q表示四元数体,Qnn表示 n n阶四元数矩阵集合,QI1I2IN表示N阶I1 I2 IN维四元数张量集合;AT表示矩阵A的转置,R(A)表示矩阵A的列空间.2.预备知识定义2.1四元数a可以表示为a=a1+a2i+a3j+a4k Q,收稿日期：2023-02-01基金项目：国家自然科学基金资助项目(12371378,11901098);福建省自然科学基金(2023J011127,2020J05034)作者简介：胡晶晶,女,汉族,河南人,研究方向:控制论中的矩阵计算.通讯作者：柯艺芬.第 1 期胡晶晶等：修正共轭梯度算法求解四元数Sylvester张量方程239其中as R,s=1,2,3,4,且i,j,k满足i2=j2=k2=1,ij=ji=k,jk=kj=i,ki=ik=j.对于两个四元数a=a1+a2i+a3j+a4k Q和b=b1+b2i+b3j+b4k Q,它们的乘积定义为ab=a1b1 a2b2 a3b3 a4b4+(a1b2+a2b1+a3b4 a4b3)i+(a1b3 a2b4+a3b1+a4b2)j+(a1b4+a2b3 a3b2+a4b1)k.定义2.2四元数张量A可以表示为A=A1+A2i+A3j+A4k QI1I2IN,其中As RI1I2IN,s=1,2,3,4.四元数张量A的共轭可被定义为A=A1 A2i A3j A4k QI1I2IN.类似地,对于四元数矩阵A,有着相似的表示.定义2.310给定一个N阶实张量A RI1I2IN.将N阶张量元素ai1i2iN映射为模-n矩阵A(n)RInI1In1In+1IN的元素ain,j,其中j=1+Nk=1,k=n(ik 1)k1m=1,m=nIm.例2.110张量A R342的两个正面切片分别为A(:,:,1)=147102581136912,A(:,:,2)=131619221417202315182124,张量A沿模-1展开成矩阵A(1),其中A(1)=147101316192225811141720233691215182124.定义2.411对于张量A RI1I2IN,算子vec将张量的列叠加成一个向量,定义vec(A)=vec(A(1),其中A(1)是张量A沿模-1的展开矩阵.例2.2若张量A R342如例2.1所示,则vec(A)=1,2,24T.定义2.510一个N阶实张量X RI1I2IN与一个实矩阵A(n)RJnIn的n-模式积用符号X nA(n)表示,这是一个定义为I1 In1 Jn In+1 IN的实张量,其元素定义为(X nA(n)i1in1jnin+1iN=Inin=1xi1iniNajnin.定义2.610对于实张量X,Y RI1I2IN,张量X和Y的内积定义为X,Y=vec(X)Tvec(Y)=i1,i2,iNxi1i2iNyi1i2iN.性质2.1对于实张量X,Y,Z RI1I2IN,由定义2.6所定义的内积满足以下性质:1)对称性:X,Y=Y,X;2)线性性:k1X+k2Y,Z=k1X,Z+k2Y,Z,其中k1,k2 R;3)正定性:对于所有的X=O,X,X 0,其中O为RI1I2IN上的零张量,即O的所有元素为0.240应用数学2024张量X RI1I2IN的Frobenius范数 定义为X:=X,X=i1,i2,iN|xi1i2iN|2=vec(X),vec(X).性质2.210若X RI1In1JIn+1IN和Y RI1In1KIn+1IN,且A RJK,则有X,Y nA=X nAT,Y.引理2.19线性映射L:RI1I2IN RI1I2IN定义为:L(X)=X 1A(1)+X 2A(2)+X NA(N),(2.1)其中A(n)RInIn(n=1,2,N),X RI1I2IN.对于任意的张量X,Y RI1I2IN,根据性质2.1和性质2.2,则L与其共轭线性算子LT满足L(X),Y=X,LT(Y),这里LT(Y)=Y 1(A(1)T+Y 2(A(2)T+Y N(A(N)T.定义2.7 四元数张量A=A1+A2i+A3j+A4k QI1I2IN,B=B1+B2i+B3j+B4k QI1I2IN,A=B的充要条件为As=Bs,其中As,Bs RI1I2IN(s=1,2,3,4).本文主要研究如下四元数域上的Sylvester张量方程:X 1A(1)+X 2A(2)+X NA(N)+X 1B(1)+X 2B(2)+X NB(N)=C,(2.2)其中四元数矩阵A(n)QInIn,B(n)QInIn(n=1,2,N),C QI1I2IN为已知四元数张量,X QI1I2IN为未知四元数张量.定理2.1对于四元数矩阵A(n)=A(n)1+A(n)2i+A(n)3j+A(n)4k QInIn,B(n)=B(n)1+B(n)2i+B(n)3j+B(n)4k QInIn(n=1,2,N),四元数张量C=C1+C2i+C3j+C4k QI1I2IN和X=X1+X2i+X3j+X4k QI1I2IN.四元数Sylvester张量方程(2.2)等价于如下张量方程组:LA(n)1+B(n)1(X1)LA(n)2B(n)2(X2)LA(n)3B(n)3(X3)LA(n)4B(n)4(X4)=C1,LA(n)2+B(n)2(X1)+LA(n)1B(n)1(X2)+LA(n)4B(n)4(X3)LA(n)3B(n)3(X4)=C2,LA(n)3+B(n)3(X1)LA(n)4B(n)4(X2)+LA(n)1B(n)1(X3)+LA(n)2B(n)2(X4)=C3,LA(n)4+B(n)4(X1)+LA(n)3B(n)3(X2)LA(n)2B(n)2(X3)+LA(n)1B(n)1(X4)=C4,(2.3)其中LA(n)sB(n)s(X)=X1(A(1)sB(1)s)+X2(A(2)sB(2)s)+XN(A(N)sB(N)s),s=1,2,3,4.此外,实张量方程组(2.3)等价于下列线性方程组x=b,(2.4)=Kro(LA(n)1+B(n)1)Kro(LA(n)2B(n)2)Kro(LA(n)3B(n)3)Kro(LA(n)4B(n)4)Kro(LA(n)2+B(n)2)+Kro(LA(n)1B(n)1)+Kro(LA(n)4B(n)4)Kro(LA(n)3B(n)3)Kro(LA(n)3+B(n)3)Kro(LA(n)4B(n)4)+Kro(LA(n)1B(n)1)+Kro(LA(n)2B(n)2)Kro(LA(n)4+B(n)4)+Kro(LA(n)3B(n)3)Kro(LA(n)2B(n)2)+Kro(LA(n)1B(n)1),(2.5)x=vec(X1)vec(X2)vec(X3)vec(X4),b=vec(C1)vec(C2)vec(C3)vec(C4),第 1 期胡晶晶等：修正共轭梯度算法求解四元数Sylvester张量方程241Kro(LA(n)sB(n)s)=Nn=1I(IN)I(In+1)(A(n)sB(n)s)I(In1)I(I1),s=1,2,3,4,其中I(m)代表m阶单位矩阵.证四元数Sylvester张量方程(2.2)可写作Nn=1(X nA(n)+X nB(n)=C,(2.6)其中X nA(n)=X1nA(n)1 X2nA(n)2 X3nA(n)3 X4nA(n)4+X1nA(n)2+X2nA(n)1+X3nA(n)4 X4nA(n)3i+X1nA(n)3 X2nA(n)4+X3nA(n)1+X4nA(n)2j+X1nA(n)4+X2nA(n)3 X3nA(n)2+X4nA(n)1k,(2.7)X nB(n)=X1nB(n)1+X2nB(n)2+X3nB(n)3+X4nB(n)4+X1nB(n)2 X2nB(n)1 X3nB(n)4+X4nB(n)3i+X1nB(n)3+X2nB(n)4 X3nB(n)1 X4nB(n)2j+X1nB(n)4 X2nB(n)3+X3nB(n)2 X4nB(n)1k.(2.8)将式(2.7)和式(2.8)代入式(2.6),根据定义2.7,四元数张量相等意味着实部和虚部分别对应相等,从而可得方程组(2.3).进一步,根据Kronecker积和向量化算子,从而可得线性方程组(2.4).注2.1值得注意的是当线性方程组(2.4)的系数矩阵规模非常大时,寻找其数值解将面临很大的计算挑战.文12表明,基于张量格式的算法总体上比经典形式的算法效率更高.下面我们将构造一类基于张量形式的高效迭代算法求解方程组(2.3),并建立迭代算法的收敛结果.注2.2记LTA(n)sB(n)s(X)=X 1(A(1)s B(1)s)T+X 2(A(2)s B(2)s)T+X N(A(N)s B(N)s)T,其中s=1,2,3,4.根据Kronecer积的性质,可得Kro(LTA(n)sB(n)s)=(Kro(LA(n)sB(n)s)T.3.全局收敛性本节主要考虑四元数张量方程(2.2)的求解.利用定理2.1,可以将四元数张量方程转化为实张量方程组(2.3),考虑用张量形式的修正共轭梯度(Modify Conjugate Gradient,MCG)算法求解(2.3),从而得到方程(2.2)的解.为了方便起见,我们引入下列线性算子:M11=LA(n)1+B(n)1,M12=LA(n)2B(n)2,M13=LA(n)3B(n)3,M14=LA(n)4B(n)4,M21=LA(n)2+B(n)2,M22=LA(n)1B(n)1,M23=LA(n)4B(n)4,M24=LA(n)3B(n)3,M31=LA(n)3+B(n)3,M32=LA(n)4B(n)4,M33=LA(n)1B(n)1,M34=LA(n)2B(n)2,M41=LA(n)4+B(n)4,M42=LA(n)3B(n)3,M43=LA(n)2B(n)2,M44=LA(n)1B(n)1.从而四元数Sylvester张量方程(2.2)可写成4t=1Mst(Xt)=Cs,s=1,2,3,4,(3.1)242应用数学2024其中Xt RI1I2IN,Cs RI1I2IN.算法3.1MCG算法求解张量方程组(3.1)步0输入矩阵A(n)s RInIn,B(n)s RInIn(n=1,2,N,s=1,2,3,4),张量CsRI1I2IN(s=1,2,3,4).选择初始张量Xt(0)RI1I2IN(t=1,2,3,4).步1对于s=1,2,3,4,t=1,2,3,4,计算Rs(0)=Cs4t=1Mst(Xt(0),Qt(0)=4s=1MTst(Rs(0).置k:=0.步2若4s=1Rs(k)2=0 或者4t=1Qt(k)2=0,则停止;否则,转步3.步3对于s=1,2,3,4,t=1,2,3,4,计算k=4s=1Rs(k)24t=1Qt(k)2,Xt(k+1)=Xt(k)+kQt(k),Rs(k+1)=Rs(k)k4t=1Mst(Qt(k),k=4s=1Rs(k+1)24s=1Rs(k)2,Qt(k+1)=4s=1MTst(Rs(k+1)+kQt(k).步4置k:=k+1,返回步2.在算法3.1中,MTst是线性算子Mst的共轭线性算子.根据引理2.1可得,对任给的张量XtRI1I2IN,Ys RI1I2IN,满足Mst(Xt),Ys=Xt,MTst(Ys).(3.2)定理3.1假设序列Rs(k)和Qt(k)(s=1,2,3,4,t=1,2,3,4)由算法3.1生成,则下列等式成立4s=1Rs(u),Rs(v)=0,u,v=0,1,2,u=v,(3.3)4t=1Qt(u),Qt(v)=0,u,v=0,1,2,u=v.(3.4)证首先,我们证明式(3.3)和式(3.4)对于0 u v k 成立,对k用数学归纳法,当k=1时,根据算法3.1和式(3.2),我们有4s=1Rs(0),Rs(1)=4s=1Rs(0),Rs(0)04t=1Mst(Qt(0)=4s=1Rs(0)2 04s=14t=1Rs(0),Mst(Qt(0)第 1 期胡晶晶等：修正共轭梯度算法求解四元数Sylvester张量方程243=4s=1Rs(0)2 04t=14s=1MTst(Rs(0),Qt(0)=4s=1Rs(0)2 04t=1Qt(0)2=0和4t=1Qt(0),Qt(1)=4t=1Qt(0),4s=1MTst(Rs(1)+1Qt(0)=4s=14t=1Mst(Qt(0),Rs(1)+14t=1Qt(0)2=4s=1Rs(0)Rs(1)0,Rs(1)+14t=1Qt(0)2=104s=1Rs(1)2+14t=1Qt(0)2=0.从而,当k=1时,式(3.3)和式(3.4)成立.现假设式(3.3)和式(3.4)对于k=l成立.对于k=l+1,根据算法3.1和式(3.2),我们有4s=1Rs(l),Rs(l+1)=4s=1Rs(l),Rs(l)l4t=1Mst(Qt(l)=4s=1Rs(l)2 l4s=14t=1Rs(l),Mst(Qt(l)=4s=1Rs(l)2 l4t=14s=1MTst(Rs(l),Qt(l)=4s=1Rs(l)2 l4t=1Qt(l)l1Qt(l 1),Qt(l)=4s=1Rs(l)2 l4t=1Qt(l)2=0和4t=1Qt(l),Qt(l+1)=4t=1Qt(l),4s=1MTst(Rs(l+1)+lQt(l)=4s=14t=1Mst(Qt(l),Rs(l+1)+l4t=1Qt(l)2=4s=1Rs(l)Rs(l+1)l,Rs(l+1)+l4t=1Qt(l)2=1l4s=1Rs(l+1)2+l4t=1Qt(l)2=0.最后,我们证明当u l时,式(3.3)和式(3.4)成立.事实上,我们有4s=1Rs(u),Rs(l+1)=4s=1Rs(u),Rs(l)l4t=1Mst(Qt(l)244应用数学2024=4s=1Rs(u),Rs(l)l4s=14t=1Rs(u),Mst(Qt(l)=l4t=14s=1MTst(Rs(u),Qt(l)=l4t=1Qt(u)u1Qt(u 1),Qt(l)=0和4t=1Qt(u),Qt(l+1)=4t=1Qt(u),4s=1MTst(Rs(l+1)+lQt(l)=4s=14t=1Mst(Qt(u),Rs(l+1)=4s=1Rs(u)Rs(u+1)u,Rs(l+1)=0.从而,当k=l+1时,式(3.3)和式(3.4)也成立,则对所有的0 u v时,式(3.3)和式(3.4)也成立.这样对所有的u=v,结论得证.定理3.2记Xt(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)满足4t=1Qt(k),Xt Xt(k)=4s=1Rs(k)2,k=0,1,2,.(3.5)证用数学归纳法.当k=0时,有4t=1Qt(0),Xt Xt(0)=4s=14t=1MTst(Rs(0),Xt Xt(0)=4s=1Rs(0),4t=1Mst(Xt Xt(0)=4s=1Rs(0),Cs4t=1Mst(Xt(0)=4s=1Rs(0)2.假设当k=l(l 1)时,式(3.5)成立.对于k=l+1,我们可得4t=1Qt(l+1),Xt Xt(l+1)=4t=14s=1MTst(Rs(l+1)+lQt(l),Xt Xt(l+1)=4s=1Rs(l+1),4t=1Mst(Xt Xt(l+1)+l4t=1Qt(l),Xt Xt(l)lQt(l)=4s=1Rs(l+1),Cs4t=1Mst(Xt(l+1)+l4t=1Qt(l),Xt Xt(l)ll4t=1Qt(l),Qt(l)第 1 期胡晶晶等：修正共轭梯度算法求解四元数Sylvester张量方程245=4s=1Rs(l+1)2.从而当k=l+1时,式(3.5)成立.由数学归纳法可得,式(3.5)对所有的k=0,1,2,都成立.2定理3.3若实张量方程组(3.1)是相容的,对任意的初始张量X1(0),X2(0),X3(0),X4(0),在不计舍入误差的情况下,其精确解可通过算法3.1在有限迭代步内获得.证若4s=1Rs(k)2=0,k=0,1,2,4I 1(I=I1I2IN),根据定理3.2,有4t=1Qt(k)2=0.通过算法3.1可得到Xt(4I)和Rs(4I),令R(k)=diagR1(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(u),R(v)=4s=1Rs(u),Rs(v)=0,u,v=0,1,2,4I 1,u=v,则R(0),R(1),R(4I 1)可构成如下线性空间的一组正交基,H=E|E=E10000E20000E30000E4,其中Es RI1I2I3IN(s=1,2,3,4).根据定理3.1,有R(k),R(4I)=4s=1Rs(k),Rs(4I)=0,k=0,1,2,4I 1.故Rs(4I)=O.则X1(4I),X2(4I),X3(4I),X4(4I)是实张量方程组(3.1)的解.4.极小范数解引理4.113若A Rmn,b Rm,线性方程组Ax=b有解x且x R(AT),则x是Ax=b的唯一极小范数解.定理4.1若实张量方程组(3.1)是相容的,如果选择初始张量为Xt(0)=4s=1MTst(Hs(0),t=1,2,3,4,(4.1)其中Hs(0)RI1I2IN(s=1,2,3,4)是任意张量(特别地,我们选择Xt(0)=O,t=1,2,3,4),则由算法3.1生成的解X1,X2,X3,X4是实张量方程组(3.1)的唯一极小Frobenius范数解.证若选择初始张量如式(4.1)所定义,根据算法3.1,存在张量Hs(k)RI1I2IN(s=1,2,3,4)满足Xt(k)=4s=1MTst(Hs(k),根据定理3.3,存在张量Hs RI1I2IN(s=1,2,3,4)满足Xt=4s=1MTst(Hs).(4.2)246应用数学2024下面我们证明解Xt(t=1,2,3,4)是唯一极小Frobenius范数解.根据Kronecker积和注2.2,我们有x=vec(X1)vec(X2)vec(X3)vec(X4)=Kro(MT11)Kro(MT21)Kro(MT31)Kro(MT41)Kro(MT12)Kro(MT22)Kro(MT32)Kro(MT42)Kro(MT13)Kro(MT23)Kro(MT33)Kro(MT43)Kro(MT14)Kro(MT24)Kro(MT34)Kro(MT44)vec(H1)vec(H2)vec(H3)vec(H4)=+Kro(LTA(n)1+B(n)1)+Kro(LTA(n)2+B(n)2)+Kro(LTA(n)3+B(n)3)+Kro(LTA(n)4+B(n)4)Kro(LTA(n)2B(n)2)+Kro(LTA(n)1B(n)1)Kro(LTA(n)4B(n)4)+Kro(LTA(n)3B(n)3)Kro(LTA(n)3B(n)3)+Kro(LTA(n)4B(n)4)+Kro(LTA(n)1B(n)1)Kro(LTA(n)2B(n)2)Kro(LTA(n)4B(n)4)Kro(LTA(n)3B(n)3)+Kro(LTA(n)2B(n)2)+Kro(LTA(n)1B(n)1)vec(H1)vec(H2)vec(H3)vec(H4)=Kro(LA(n)1+B(n)1)Kro(LA(n)2B(n)2)Kro(LA(n)3B(n)3)Kro(LA(n)4B(n)4)Kro(LA(n)2+B(n)2)+Kro(LA(n)1B(n)1)+Kro(LA(n)4B(n)4)Kro(LA(n)3B(n)3)Kro(LA(n)3+B(n)3)Kro(LA(n)4B(n)4)+Kro(LA(n)1B(n)1)+Kro(LA(n)2B(n)2)Kro(LA(n)4+B(n)4)+Kro(LA(n)3B(n)3)Kro(LA(n)2B(n)2)+Kro(LA(n)1B(n)1)Tvec(H1)vec(H2)vec(H3)vec(H4)R(T).由于线性方程组(2.4)和实张量方程组(3.1)等价,根据引理4.1可得x为(2.4)的唯一极小范数解,从而由算法3.1生成的解Xt(t=1,2,3,4)是实张量方程组(3.1)的唯一极小Frobenius范数解.5.数值实验为了验证本文算法的可行性,本节通过三个例子来说明所提出算法的效果.所有的程序都是基于Bader和Kolda14开发的MATLAB张量工具箱,实验中,选取N=3,当满足4s=1Rs(k)2 106或4t=1Qt(k)2 106,(5.1)迭代终止.所有算法均在带有MATLAB R2018a的个人PC机上运行.此外,定义残差范数RES(k):=4s=1Rs(k)2.例5.1考虑四元数Sylvester张量方程X 1A(1)+X 2A(2)+X 3A(3)+X 1B(1)+X 2B(2)+X 3B(3)=C,(5.2)其中A(1)=i+j+k2 k1 i+jii+j+2k1 k2i+j1+ki+j+4k,B(1)=1+i+jk1+i+k4i+j+2k1+i+j+k1+k1+k1+i+j1+i+k,A(2)=2i+j+k1+ij1 j+2ki+jkik1+i j,B(2)=i1+i1 j+k2j3+i+2k1+j+k2+i+j+k2k1+i+j,A(3)=1+i+k2+jkji+j+2k1 ki+j+k01+i,B(3)=1+j+2k2+ii+j+k1+j1+i+2j+kj kj+k1+k1+i,C(:,:,1)=1+i+j+k2+i+2j+k1+j+2k1+i+j+k3kii+j1+i+2j1,第 1 期胡晶晶等：修正共轭梯度算法求解四元数Sylvester张量方程247C(:,:,2)=3+j+k1+2i1+i+j+ki1+j+ki+kj+2k1+i+jj+k,C(:,:,3)=1+ji1+i+j+k1+j+k1+i+k2i+j1+j+k1+i1+i+k.考虑初始张量Xt(0)=O(t=1,2,3,4),当满足式(5.1)时迭代终止.根据算法3.1,我们可获得式(5.2)的极小Frobenius范数解,数值结果如图1所示.050100150200Iterative number10-610-410-2100102104RES(k)图1例5.1的数值结果020406080100120140160Iterative number10-610-410-2100102104106108RES(k)图2例5.2的数值结果例5.2考虑四元数Sylvester张量方程(5.2),设A(1)=A11+A12i+A13j+A14k,B(1)=B11+B12i+B13j+B14k,A(2)=A21+A22i+A23j+A24k,B(2)=B21+B22i+B23j+B24k,A(3)=A31+A32i+A33j+A34k,B(3)=B31+B32i+B33j+B34k,这里A11=triu(hilb(n),A12=triu(ones(n,n),A13=eye(n),A14=zeros(n,n),A21=tridiag(n,1,2,1),A22=eye(n),A23=zeros(n,n),A24=tridiag(n,0.5,6,0.5),A31=A32=A33=A34=ones(n,n),B11=hankel(1:n),B12=ones(n,n),B13=B14=zeros(n,n),B21=ones(n,n),B22=tridiag(n,1,2,1),B23=B24=eye(n),B31=eye(n,n),B32=B33=B34=zeros(n,n).令精确解X=tenones(n,n,n)+tenones(n,n,n)i+tenones(n,n,n)j+tenones(n,n,n)k.从而可得式(5.2)的右端项张量C.在这个例子中,张量方程(5.2)是相容的,通过选择初始张量Xt(0)=O(t=1,2,3,4),对于n=10时,执行算法3.1得到的收敛曲线如图2所示.下面的例子通过求解稀疏四元数Sylvester张量方程来检测算法3.1的有效性.例5.3考虑四元数Sylvester张量方程X 1A(1)+X 2A(2)+X 3A(3)=C,(5.3)其中系数矩阵满足A(1)=A11+A12i+A13j+A14k,A(2)=A21+A22i+A23j+A24k,248应用数学2024A(3)=A31+A32i+A33j+A34k,A11=A24=A33,A13=A22=A31,A14=A23=A32,A21=A34,S1=tridiag(1,2,1),S2=tridiag(1,4,1),S3=tridiag(1,3,1),S4=tridiag(1,4,1).与Ahmadi-Asl和Beik8工作稍有不同,我们考虑以下四种情况:情况1A11=S1,A12=pivtol,A21=S2,A13=S3,A14=S4;情况2A11=S1,A12=gre343,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.050100nz=306050100pivtol图5.1第一个测试矩阵的稀疏结构0100200300nz=13100100200300gre_343图5.2第二个测试矩阵的稀疏结构0100200300nz=13770100200300mesh3em5图5.3第三个测试矩阵的稀疏结构050100150200250nz=1794050100150200250sphere3图5.4第四个测试矩阵的稀疏结构123456789Iterative number10-810-610-410-2100102104106RES(k)Case 1Case 2Case 3Case 4图5.5例5.3的收敛结果第 1 期胡晶晶等：修正共轭梯度算法求解四元数Sylvester张量方程249测试矩阵来自于Davis收集的矩阵15.数据来源于HB(Harwell-Boeing)group.在这个例子中出现的这四种测试矩阵(pivtol,gre343,mesh3em5和sphere3)经常出现在统计、有向加权图和结构问题中.这些矩阵的稀疏结构如图5.1-5.4所示,测试矩阵的性质见表5.1,其中nz表示非零项的个数.需要注意的是,这里的稀疏矩阵的阶数较大,因此四元数Sylvester张量方程(5.3)的解的维数非常大.令X=tendiag(ones(1,n),n,n,n)+tendiag(ones(1,n),n,n,n)i+tendiag(ones(1,n),n,n,n)j+tendiag(ones(1,n),n,n,n)k,其中对于情况1-4,n分别取为102,343,289和258作为其精确解.通过选择初始张量Xt(0)=O,应用算法3.1去求解张量方程(5.3),当满足式(5.1)时迭代终止,相应的收敛结果如图5.5所示,结果证实了算法3.1的收敛性,可以看出残差范数是显著下降的.因此,算法3.1对于求解稀疏张量方程是有效的.表5.1例5.3测试矩阵的性质矩阵行数列数结构nz条件数pivtol102102非对称306109.6070gre343343343非对称1310111.9763mesh3em5289289对称13774.9660sphere3258258对称17946.4894e+166.结论本文通过将四元数Sylvester张量方程等价地转化为实数域上的张量方程组并引入线性算子,构造了基于张量形式的修正共轭梯度算法求其等价形式.理论分析表明对任给的初始张量,由算法3.1生成的迭代序列在不计舍入误差的情况下,可在有限迭代步内收敛到方程组的精确解.进一步,通过选择特殊类型的初始张量,可获得唯一极小Frobenius范数解.最后,数值例子证明了该算法的有效性.参考文献：1 ADLER S L.Quaternionic Quantum Mechanics and Quantum FieldsM.New York:Oxford Univer-sity Press,1995.2 BIHAN N L,MARS J.Singular value decomposition of matrices:a new tool for vector-sensor signalprocessingJ.Signal Processing,2004,84:1177-1199.3 ELL T A,SANGWINE S J.Hypercomplex Fourier transforms of color imagesJ.IEEE TransactionsImage Processing,2007,16(1):22-35.4 VINCE J.Quaternions for Computer GraphicsM.London:Springer Publishing Company,Incorpo-rated,2011.5 YUAN S F,WANG Q W,ZHANG X.Least-squares problem for the quaternion matrix equationAXB+CY D=E over different constrained matricesJ.International Journal of Computer Math-ematics,2013,90:565-576.6 YUAN S F,WANG Q W,XIONG Z P.The least squares-Hermitian problems of quaternion matrixequation AHXA+BHY B=CJ.Filomat,2014,28:1153-1165.7 ZHANG F X,WEI M S,LI Y,et al.An efficient real representation method for least squares problem ofthe quaternion constrained matrix equation AXB+CY D=EJ.International Journal of ComputerMathematics,2021,98:1408-1419.8 BEIK F P A,AHMADI-ASL S.An iterative algorithm for-(anti)-Hermitian least-squares solutionsof quaternion matrix equationsJ.Electronic Journal of Linear Algebra,2015,80:1081-3810.250应用数学20249 WANG Q W,XU X J,DUAN X F.Least squares solution of the quaternion Sylvester tensor equa-tionJ.Linear Multilinear Algebra,2021,69:104-130.10 BADER B W.Tensor decompositions and applicationsJ.SIAM Review,2009,51(3):455-500.11 XU X J,WANG Q W.Extending BiCG and BiCR methods to solve the Stein tensor equationJ.Computers and Mathematics with Applications,2019,77:3117-3127.12 BEIK F P A,SABERI-MOVAHED F,AHMADI-ASL S.On the Krylov subspace methods basedon tensor format for positive definite Sylvester tensor equationsJ.Numerical Linear Algebra withApplications,2016,23:444-466.13 PENG Y X,HU X Y,ZHANG L.An iteration method for the symmetric solutions and the optimalapproximation solution of the matrix equation AXB=CJ.Journal of Computational and AppliedMathematics,2005,160:763-777.14 BADER B W,KOLDA T G.MATLAB Tensor Toolbox,Version 2.6EB/OL.(2015-02-06)2023-11-03.https:/www.tensortoolbox.org.15 UniversityofFloridaSparseMatrixCollectionwebpageEB/OL.(2015-06-10)2023-11-03.https:/www.cise.ufl.edu/research/sparse/matrices/listbyid.html.Modified Conjugate Gradient Algorithm for Solving theQuaternion Sylvester Tensor EquationHU Jingjing1,2,3,KE Yifen1,2,3,MA Changfeng1,2,3