矩陣方程的一個二次收斂迭代算法*

劉 健，孫洪春

(1.魯南技師學院信息工程系，山東 臨沂 276021；2. 臨沂大學數學與統計學院，山東 臨沂 276005)

矩陣方程的一個二次收斂迭代算法*

劉 健1，孫洪春2

(1.魯南技師學院信息工程系，山東 臨沂 276021；2. 臨沂大學數學與統計學院，山東 臨沂 276005)

提出了求解一類矩陣方程的一個收斂算法, 在較寬松的條件下, 證明了所給算法的全局收斂性和二次收斂率,也給出了算法的數值試驗,試驗表明算法是有效的.

矩陣方程；全局收斂性；二次收斂率

引言

考察的矩陣方程是：

(1)

其中Ai∈Cm×p,Bi∈Cq×n,C∈Cm×n為已知矩陣,X∈Cp×q為未知矩陣.為敘述方便，令

(2)

F(X):Cpq→Cmn.在本文中,假設方程(1)的解集非空,記為Ω*.

矩陣方程是代數學中重要的問題之一, 廣泛應用于計算數學、量子力學、計算機安全、自動控制理論、系統識別、振動理論、圖像處理等.近年來，無論在實際背景應用還是在學術前景研究方面，矩陣方程解的存在性和求解算法的設計已成為數值代數研究的一個熱門課題，引起了國內外許多學者濃厚的研究興趣.Xu等[1]給出了矩陣方程最小二乘解的表達式和解的存在唯一性的充分必要條件. Sheng等[2]提出了求解矩陣方程的對稱解和斜對稱解的兩種有效的迭代方法. Liao等[3]給出了矩陣方程的一種最小二乘解的最小范數解析表達式，并給出了求解該類解的算法和數值實驗. Peng[4]提出了兩種求解矩陣方程的迭代方法，相比較該迭代方法的收斂速度更快、精度更高. Liao[5]基于廣義奇異值分解探討了半正定矩陣集上的矩陣方程最小二乘解. 彭亞新[6]用迭代法研究了矩陣方程AXB=F的一般解、對稱解、(反)中心對稱解、雙對稱解與對稱次反對稱解等問題. Li等[7]基于Moore-Penrose逆和Kronecker積，提出了求解矩陣方程解的一般表達式和相應的數值算法. 本文提出了求解矩陣方程的一個新矩陣迭代算法, 在較寬松的條件下, 證明了所給算法的全局收斂性和二次收斂率，也給出了算法的數值試驗,試驗表明算法是有效的.

本文所用符號含義：用‖·‖ 表示向量的2范數，‖·‖F表示矩陣的Frobenius范數，dist(X,Ω*)表示矩陣X到解集Ω*的最近距離，AT表示矩陣A的轉置.

1 預備知識

為分析所給下面迭代算法， 回顧如下的一些知識.

定義1 設矩陣A=(aij)m×n∈Rm×n,B=(bij)s×t∈Rs×t, 則A和B的Kronecker積定義為

定義2 設矩陣A=(aij)m×n∈Rm×n， 則A的拉直就是按矩陣A的列的先后順序進行上下疊放，記為vec(A)，即

對于以上兩種運算，有以下性質.

命題1 對可乘矩陣A,B,C,D,下述結論成立

(A1) (A?B)T=AT?BT； (A2) vec(ABC)=(CT?A)vec(B)

由式(2)和命題1的(A2)知:

(3)

vecF(X)=AvecX-vecC=0

(4)

因此，Ω*和Ω是等價的，其中Ω={vecX∈Rpq|AvecX=vecC}.

設矩陣M,N∈Rm×n,定義矩陣的內積為:

因此，可定義矩陣M∈Rm×n的Frobenius范數：

(5)

命題2 對于集合W={x∈Rn|Dx=d}，D∈Rl×n,d∈Rl. 則

PW(x)=(I-D+D)x+D+d, ?x∈Rn,

式中PW(x)是向量x∈Rn在集合W上的投影，D+是矩陣D的廣義逆.

由命題2，結合(4)可得:

(6)

根據式(5)的定義，可得:

dist(X,Ω*)≤‖A+‖‖F(X)‖F

(7)

還可得到如下性質. 基于(3)，容易得vecF′(X)=A. 對于(4)，有:

(8)

式中a>0為常數.

?X,Y∈Rpq.由(3)得：

(9)

再根據(5)的定義方式，?X,Y∈Rpq.進一步可得：

(10)

2 算法和收斂性

本節，給出解(1)的一個新矩陣迭代算法，基于預備知識中的分析，證明所給算法的全局收斂性和二次收斂率.

算法1

第一步:選擇參數ε≥0，初始點X0∈Rp q. 令k=0.

第二步:若‖F(Xk‖)≤ε，停止；否則，轉第三步 3.

第三步:求方程組的解Dk作為在點Xk處的一個搜索方向

(ATA+μkI)vecDk=-ATvecF(Xk)

(11)

下面，給出算法1的收斂性和收斂率分析. 首先分析算法1的可行性.

對于方程(11)，由μk>0，則矩陣(ATA+μkI)是正定的，所以(11)有唯一解Dk. 而且vecDk也是下面無約束優化的一個下降方向，

(12)

▽φ(Xk)TvecDk=[vec′F(Xk)TvecF(Xk)]TvecDk=-[ATvecF(Xk)]T(ATA+μkI)-1[ATvecF(Xk)]<0.

定義θk(D)=‖AvecD+vecF(Xk)‖2+μk‖vecD‖2，函數θk為嚴格凸函數,因此,無約束極小化問題

minθk(D)

(13)

等價于(11). 事實上，

▽θk(Dk)=2AT[AvecDk+vecF(Xk)]+2μkvecDk=2[ATA+2μkI]vecDk+2ATvecF(Xk)=0.

引理1 對于算法1產生的Xk,X*∈Ω*，則由(11)求得的解Dk滿足

‖Dk‖≤dist(Xk,Ω*)

(14)

‖AvecDk+vecF(Xk)‖≤c1dist2(Xk,Ω*)

(15)

證明 因為vecDk是最小化問題(13)的解,從而有：

(16)

根據θk的定義及(16)、(8)有：

同理可證第二個不等式

3.1.4 豬的室間隔缺損模型 這類模型非常重要，由于這類疾病模型屬于自然發生，能夠有效模擬人類收縮功能正常的心力衰竭模型[18-19]。

引理2 對于算法1產生的Xk-1,Xk,有:

dist(Xk,Ω*)≤c2dist2(Xk-1,Ω*)

(17)

式中c2=c1‖A+‖.

證明 由式(8)得

0 = ‖vecF'(Xk-1)vecDk-1-(vecF(Xk-1+Dk-1)-vecF(Xk-1))‖≥

‖vecF(Xk-1+Dk-1)‖-‖vecF'(Xk-1)vecDk-1+ vecF(Xk-1)‖

即‖vecF(Xk-1+Dk-1)‖≤‖vecF'(Xk-1)vecDk-1+ vecF(Xk-1)‖.再由Xk=Xk-1+Dk-1知：

dist(Xk,Ω*) = dist(Xk-1+Dk-1,Ω*)≤‖A+‖‖F(Xk-1+Dk-1)‖F=

‖A+‖‖vecF(Xk-1+Dk-1)‖≤‖A+‖‖AvecDk-1+ vecF(Xk-1)‖≤c1‖A+‖dist2(Xk-1,Ω*).

式中第一個不等式是根據式(7)而得，最后不等式成立是因為式(15).令c2=c1(A+),有式(17)成立.

證明 首先證明l=1時結論成立.

‖X1-X*‖≤‖X0+D0-X*‖ ≤‖X0-X*‖ + ‖D0‖ ≤

‖X0-X*‖ + dist(X0,Ω*) ≤‖X0-X*‖ + ‖X0-X*‖ ≤2r.

(18)

再由(13)式得：

(19)

于是有：

(20)

基于以上引理,下面給出點列{Xk}的收斂定理.

{dist(Xk,Ω*)}二次收斂于0.

從而{Xk}為柯西點列, 因此點列{Xk}收斂.

算法2

第二步: 若‖F(Xk)‖≤ε，停止；否則，轉第三步.

第三步: 求下面方程組的解Dk使得(ATA+μkI)vecDk=-ATvecF(Xk).

若

‖F(Xk+Dk)‖F≤γ‖F(Xk)‖F

(21)

第四步: 令mk是使得下面不等式成立的最小非負整數m,

φ(Xk+βmDk)-φ(Xk)≤αβm▽φ(Xk)TDk

(22)

下面，給出算法2的全局收斂性和二次收斂率.

定理2設{Xk}是算法2產生的序列,則序列{Xk}的任一聚點都是方程(12)的解,即是方程(1)的解. 進一步有{dist(Xk,Ω*)}二次收斂于0.

證明假設方程(12)在Xk點的穩點▽φ(Xk)≠0, 則ATvecF(Xk)≠0, 又由(11)知Dk≠0, 且

▽φ(Xk)TvecDk=[vec′F(Xk)TvecF(Xk)]TvecDk=-[vecDk]T(ATA+μkI)-1[vecDk]<0.

▽φ(Xk)TvecDk=[vec′F(Xk)TvecF(Xk)]TvecDk

(1-βmk)‖(ATA+μkI)-1AT‖≤c2‖A‖‖A+‖

(21)

(22)

3 數值試驗

下面給出算法1的數值試驗,實驗表明所提算法是有效的.

例2 令F(X)=AXB-C=0, 其中:

此時‖F(X*)‖F=5.9669e-014.

[1]Xu Guiping, Wei Musheng, Zheng Daosheng. On solutions of matrix equation.AXB+CYD=F[J]. Linear Algebra and its Applications, 1998, 279: 93-109.

[2]Sheng Xingping, Chen Guoliang. An iterative method for the symmetric and skew symmetric solutions of a linear matrix equationAXB+CYD=E[J]. Journal of Computational and Applied Mathematics, 2010, 233: 3030-3040.

[3]Liao Anping, Lei Yuan. Least-Squares Solution with the Minimum-Norm for the Matrix Equation (AXB,GXH)=(C,D)[J]. Computers and Mathematms with Apphcatmns, 2005 ,50: 539-549.

[4]Peng Zhenyun. New matrix iterative methods for constraint solutions of the matrix equationAXB=C[J]. Journal of Computational and Applied Mathematics, 2010, 235:726-735.

[5]Liao A P, Bai Z Z. Least squares solutions ofAXB=Dover symmetric positive semidefinte mareices [J]. Journal of Computational Mathematics. 2003,21(2):175-182.

[6]彭亞新.求解約束矩陣方程及其最佳逼近的迭代法研究[D]. 長沙: 湖南大學數學與計量經濟學院, 2004.

[7]Li Hongyi, Gao Zongsheng, Zhao Di. Least squares solutions of the matrix equationAXB+CYD=Ewith the least norm for symmetric arrowhead matrices[J]. Applied Mathematics and Computation, 2014, 226(1): 719-724.

[8]程云鵬，矩陣論[M].西安：西北工業大學出版社，1999.

[9]王宜舉，修乃華.非線性最優化理論與算法[M].北京:科學出版社, 2016.

[10]Facchinei F. Minimization of SC1 functions and the Maratos effect [ J ]. Operations Research Letters, 1995, 17(3):131-137.

AQuadraticConvergenceIterativeAlgorithmforMatrixEquations

LIU Jian1, SUN Hong-chun2

(1.Department of Information Engineering, Lunan Technician University, Linyi Shandong 276021, China;2. School of Mathematics and Statistics, Linyi University, Linyi Shandong 276005, China)

This paper presents a convergence algorithm to solve matrix equations. Under relaxed conditions, the global convergence and quadratic convergence rate of the given algorithm are proved. The numerical experiment is also given to show that the algorithm is effective.

matrix equations; global convergence; quadratic convergence rate

1673-2103(2017)05-0001-06

2017-08-05

國家自然科學基金項目(11671228)；中國物流學會與中國物流與采購聯合會計劃項目 (2015CSLKT3-199)

劉健(1968-)，男，山東臨沂人，副教授，研究方向：最優化算法及其應用.

孫洪春(1968-), 男, 山東費縣人, 教授, 碩士，研究方向：最優化理論與算法.

O 221

A