M-矩阵代数Riccati方程的一类简单迭代法.pdf

M-矩阵代数Riccati方程的一类简单迭代法关晋瑞1，任孚鲛1，邵荣侠2（1.太原师范学院 数学与统计学院，山西 晋中 030619；2.新疆财经大学 统计与数据科学学院，新疆 乌鲁木齐 830012）摘要：文章研究了M-矩阵代数Riccati方程的数值解法。当方程的系数矩阵为正则奇异M-矩阵时，现有的一些数值方法在计算中存在一定程度的困难。为此提出了一类简单迭代法以求解方程，该方法在每步迭代中只用到矩阵乘法，运算量小且易于实现。理论分析和数值实验表明该方法是可行的，而且在一定情况下有效。关键词：代数Riccati方程；正则M-矩阵；牛顿法；迭代法中图分类号：O241.6文献标识码：A文章编号：1008-9659（2023）03-0001-05文章研究如下形式代数Riccati方程的数值解法XCX-AX-XD+B=0（1）其中A Rm m，B Rm n，C Rn m，D Rn n是已知矩阵，而X Rm n是未知矩阵，且由系数矩阵构成的K=()D-C-BA（2）是M-矩阵。这种类型的代数Riccati方程在文献中一般被称为非对称代数Riccati方程或者M-矩阵代数Riccati方程，文章沿用M-矩阵代数Riccati方程的名称。M-矩阵代数Riccati方程在中子运输理论、排队论、应用概率论和控制论等领域中具有广泛的应用1-4。式（2）中的K是非奇异的M-矩阵或者不可约奇异的M-矩阵时，人们对式（1）已进行了深入的研究，并取得了大量的成果5-13。而对于K是正则可约奇异M-矩阵的情形，相关的研究则较少，一方面是由于问题的复杂性，另一方面是由于许多经典的数值方法如牛顿法、保结构加倍算法等的收敛性很难得到保证。但在应用中这种情形也很重要，例如在Markov链的研究中，往往会遇到K是正则奇异M-矩阵的情形。GUO等人首先研究了方程最小非负解的存在性以及一些简单性质，又对保结构加倍算法的收敛性进行了深入的分析，指出在一定条件下该方法仍具有二次收敛率14-15。GUAN等人给出最小非负解存在性的一类较为简单的证明，并进行了推广16。最近，MA等人对于K是正则可约奇异M-矩阵的情形进行了研究，并提出了两类数值方法：修正牛顿法和Jacobi型不动点迭代法17。修正牛顿法的迭代格式如下：(I+A-XkC)Xk+1+Xk+1(I+D-CXk)=B-XkCXk+(+)Xk，k=0,1,2,其中X0=0，且 0，0是两个给定的参数。Jacobi型不动点迭代法的迭代格式如下：(I+A1)Xk+1+Xk+1(I+D1)=B+XkCXk+(+)Xk+A2Xk+XkD2，k=0,1,2,其中X0=0，A=A1-A2，D=D1-D2，且A1与D1分别是A与D的对角元，而 0，0是两个给定的参数。有关理论分析表明，修正牛顿法和Jacobi型不动点迭代法计算所得序列Xk都是适定的，且单调递增收敛于方程的最小非负解S17.在一般情况下，这两类方法都具有线性收敛率。尽管修正牛顿法收敛较快，但是在Vol.42，No.3Sept.2023第42卷 第3期2023年9月新疆师范大学学报（自然科学版）Journal of Xinjiang Normal University（Natural Sciences Edition）收稿日期 2023-04-11 基金项目 国家自然科学基金项目（12001395）；山西省科技创新人才团队专项资助项目（202204051002018）。作者简介 关晋瑞（1985-），男，山西襄汾人，副教授，主要从事数值代数及其应用方面研究。1新疆师范大学学报（自然科学版）2023年每步迭代中需要求解Sylvester方程，因此运算量较大；而Jacobi型不动点迭代法每步迭代运算量小，但是收敛较慢。综上所述，对于K是正则可约奇异M-矩阵的情形，目前的研究还不够深入，因此仍需要深入探讨。文章对于K是正则奇异M-矩阵的情形进行研究，并提出一类简单迭代法以求解方程。该方法在每步迭代中只用到矩阵乘法，运算量小、易于实现，而且对于K是正则奇异M-矩阵的情形具有较好的数值效果。理论分析和数值实验表明文章的方法可行，而且在一定情况下有效。1 预备知识用Rm n表示实数域上全体m n的矩阵，用A B表示矩阵的Kronecher乘积，用vec(A)表示矩阵A的拉直向量，用|表示一个给定的矩阵范数，用(A)表示矩阵A的谱半径，用I表示n阶单位矩阵。设A=(aij)Rm n，若对于任意的i j都有aij 0（aij 0）成立，则称A为非负矩阵（正矩阵），记为A 0（A 0）.设A=(aij)Rn n，若对于任意的i j都有aij 0成立，则称A为Z-矩阵。显然，任意Z-矩阵都可以表示为A=sI-B的形式，其中B为非负矩阵。若进一步有s (B)成立，则称该Z-矩阵为M-矩阵。特别地，当s (B)时称A为非奇异的M-矩阵，当s=(B)时称A为奇异的M-矩阵。定义1.114 设A Rn n为M-矩阵，若存在正向量u 0，使得Au 0，则称A为正则M-矩阵。引理1.118 设A Rn n为Z-矩阵，则下面几个结论等价：（1）A是非奇异的M-矩阵；（2）A-1 0；（3）存在正向量v 0使得Av 0；（4）A的全部特征值都具有正实部。引理1.214 设式（2）中的K是正则M-矩阵，则方程（1）存在最小非负解S，且D-CS和A-SC也是正则M-矩阵。2 一类简单迭代法对于方程（1），令A=s1I-N1,D=s2I-N2（3）其中N1和N2是非负矩阵。将上式代入方程（1），化简可得(s1+s2)X=B+XCX+N1X+XN2对上式进行迭代计算，可得一类简单迭代法Xk+1=1s1+s2()B+XkCXk+N1Xk+XkN2,k=0,1,2,（4）其中X0=0.上述方法在每步迭代中只用到矩阵加法和乘法，因此迭代简单且易于实现。当m=n时，方法的运算量为8n3，这与Jacobi型不动点迭代法的运算量相同，因此每步迭代的运算量很小。引理2.1 对于方程（1），设K为正则M-矩阵，且K 0，则由迭代格式（4）得到的序列Xk是适定的，且对任意的k 0满足0 Xk S，0 Xk Xk+1（5）其中S为方程的最小非负解。证明 由于K为正则M-矩阵，且K 0，显然有s1+s2 0.因此，序列Xk是适定的。利用数学归纳法来证明结论（5）.当k=0时，显然有X0=0 S.另外，由X1=1s1+s2B 0=X0可知结论（5）成立。2关晋瑞，等：M-矩阵代数Riccati方程的一类简单迭代法假设对于k-1结论成立，则由Xk=1s1+s2()B+Xk-1CXk-1+N1Xk-1+Xk-1N2S=1s1+s2()B+SCS+N1S+SN2可得Xk-S=1s1+s2()Xk-1CXk-1-SCS+N1Xk-1-N1S+Xk-1N2-SN2 =1s1+s2(Xk-1-S)CS+Xk-1C(Xk-1-S)+N1(Xk-1-S)+(Xk-1-S)N2 0以及Xk-Xk+1=1s1+s2()Xk-1CXk-1-XkCXk+N1Xk-1-N1Xk+Xk-1N2-XkN2 =1s1+s2(Xk-1-Xk)CXk+Xk-1C(Xk-1-Xk)+N1(Xk-1-Xk)+(Xk-1-Xk)N2 0于是结论（5）对于k也成立，从而对任意的k 0，结论（5）都成立。证毕。定理2.1 对于方程（1），设K为正则M-矩阵，K 0，S为方程的最小非负解，则由迭代格式（4）得到的序列Xk是适定的，且单调递增收敛于S.证明 根据引理2.1可知，由迭代格式（4）得到的序列Xk是适定的，单调递增而且有上界。于是Xk存在极限，记为limk Xk=X.在迭代格式（4）两端取极限，可知X为方程（1）的解。由引理2.1可知，X是非负的，且满足X S.另一方面，由于S为方程的最小非负解，于是S X.因此，有X=S.证毕。定理2.2 迭代法（4）的收敛率为limk|vec(S-Xk)|ks1+s2-s1+s2其中为A-SC的最小特征值，为D-CS的最小特征值。证明 在引理2.1的证明中，已经得到S-Xk=1s1+s2(S-Xk-1)CS+Xk-1C(S-Xk-1)+N1(S-Xk-1)+(S-Xk-1)N2于是有S-Xk=1s1+s2(S-Xk-1)(N2+CS)+(N1+Xk-1C)(S-Xk-1)1s1+s2(S-Xk-1)(N2+CS)+(N1+SC)(S-Xk-1)从而vec(S-Xk)1s1+s2I (N1+SC)+(N2+CS)T I vec(S-Xk-1).1(s1+s2)kI (N1+SC)+(N2+CS)T Ikvec(S-X0)在上式两端取范数，然后开k次方，再取极限，可得limk|vec(S-Xk)|k1s1+s2I (N1+SC)+(N2+CS)T I另外，注意到N1+SC=s1I-A+SC，N2+CS=s2I-D+CS3新疆师范大学学报（自然科学版）2023年而由引理1.2知A-SC，D-CS都是M-矩阵，于是可得limk|vec(S-Xk)|ks1+s2-s1+s2其中为A-SC的最小特征值，为D-CS的最小特征值。证毕。推论2.1 对于非临界情形，迭代法具有线性收敛率，而对于临界情形，迭代法具有半线性收敛率。证明 对于非临界情形，A-SC和D-CS中至少有一个是非奇异M-矩阵，根据引理1.1可知和至少有一个是正的，从而迭代法的收敛因子r=s1+s2-s1+s2 1，此时迭代法具有线性收敛率。对于临界情形，A-SC和D-CS都是奇异M-矩阵，这时=0，从而迭代法的收敛因子r=1，此时迭代法只具有半线性收敛率。证毕。推论2.2 在分裂（3）中，当s1=max aii，s2=max dii时，迭代法（4）具有最小的收敛因子。证明 注意到收敛因子r=s1+s2-s1+s2的表达式中当s1+s2越小时，r的值也越小。于是当s1=max aii，s2=max dii时，s1+s2最小，且满足N1=s1I-A 0，N2=s2I-D 0.证毕。3 数值例子本节通过两个数值例子对文章所提出的方法进行验证，并与文献 17 中的Jacobi型不动点迭代法进行比较。从迭代步数（IT）、计算时间（CPU）以及相对误差（RES）等三方面对两个算法的数值效果进行比较，其中相对误差定义为RES=XCX-AX-XD+BXCX+AX+XD+B实验中，所有的程序都用Matlab（R2012a）编写，并在一台个人笔记本电脑（2G CPU和8G内存）上运行，当迭代满足RES 0是参数。实验取n=256，对于参数的不同取值，利用两个方法进行计算，其中Jacobi型不动点迭代法中的参数取为1，数值实验结果如表2所示。从表中可见，两个方法都可以较快收敛，而文章的方法略快一些。4关晋瑞，等：M-矩阵代数Riccati方程的一类简单迭代法表2 例3.2的数值结果0.511.52迭代法Jacobi-FPSIMJacobi-FPSIMJacobi-FPSIMJacobi-FPSIMIT2114231526172920CPU0.38440.25890.39680.26150.41560.28690.49440.3342RES7.7860e-074.8890e-078.1495e-078.8685e-076.7125e-078.9365e-078.4290e-077.0971e-07通过上述两个例子可见，文章所提出的方法可行，而且在一定条件下与Jacobi型不动点迭代法相比，有效。参考文献：1 ZHOU K M，KHARGONEKAR P P.An Algebraic Riccati Equation Approach to H Optimization J.Systems Control Let，1988，（11）：85-91.2 ROGERS L C G.Fluid Models in Queueing Theory and Wiener-Hopf Factorization of Markov Chains J.Ann Appl Probab，1994，（04）：390-413.3 JUAN G J，LIN W W.Nonsymmetric Algebraic Riccati Equations and Hamilton-like MatricesJ.SIAM J Matrix Anal Appl，1999，（20）：228-243.4 BEAN N G，OREILLY M M，TAYLOR P G.Algorithms for Return Probabilities for Stochastic Fluid Flows J.Stochastic Models，2005，（21）：149-184.5 BINI D A，IANNAZZO B，MEINI B.Numerical Solution of Algebraic Riccati Equations M.Philadelphia：SIAM，2012：33-194.6 GUO C H.Nonsymmetric Algebraic Riccati Equations and Wiener-Hopf Factorization for M-matrices J.SIAM J Matrix Anal Appl，2001，23（01）：225-242.7 BAI Z Z，GUO X X，XU S F.Alternately Linearized Implicit Iteration Methods for the Minimal Nonnegative Solutions of the Nonsymmetric Algebraic Riccati Equations J.Numer Linear Algebra Appl，2006，13（08）：655-674.8 GUO X X，LIN W W，XU S F.A Structure-preserving Doubling Algorithm for Nonsymmetric Algebraic Riccati Equation J.Numer Math，2006，103（03）：393-412.9 WANG W G，WANG W C，LI R C.Alternating-directional Doubling Algorithm for M-matrix Algebraic Riccati Equations J.SIAM J Matrix Anal Appl，2012，33（01）：170-194.10 XUE J G，LI R C.Highly Accurate Doubling Algorithms for M-matrix Algebraic Riccati Equations J.Numer Math，2017，（135）：733-767.11 GUAN J R.Modified Alternately Linearized Implicit Iteration Method for M-matrix Algebraic Riccati Equations J.Appl Math Comput，2019，（347）：442-448.12 DONG L Q，LI J C.A Class of Complex Nonsymmetric Algebraic Riccati Equations Associated with H-matrix J.J Comput Appl Math，2020，（368）：112567.13 LIU C L，WANG W G，XUE J G，et al.Accurate Numerical Solution for Structured M-matrix Algebraic Riccati Equations J.J.Comput Appl Math，2021，（396）：113614.14 GUO C H.On Algebraic Riccati Equations Associated with M-matrices J.Linear Algebra Appl，2013，（439）：2800-2814.15 GUO C H，Lu D.On Algebraic Riccati Equations Associated with Regular Singular M-matrices J.Linear Algebra Appl，2016，（493）：108-119.16 GUAN J R，REN F J.A Note on the Minimal Nonnegative Solution for Regular M-matrix Algebraic Riccati Equations J.Mathematical Theory and Applications，2021，41（04）：77-82.17 MA C F，LU H Z.Numerical Study on Nonsymmetric Algebraic Riccati Equations J.Mediterr J Math，2016，（13）：4961-4973.18 BERMAN A，PLEMMONS R J.Nonnegative Matrices in the Mathematical SciencesM.New York：Academic Press，1994，132-158.（下转第12页）5新疆师范大学学报（自然科学版）2023年3 李椿，章立源，钱尚武.热学 M.北京：高等教育出版社，2015：16-17.4 秦允豪.热学 M.北京：高等教育出版社，2011：270-271.5 赵凯华，罗蔚茵.热学 M.北京：高等教育出版社，1998：3-5.6 高炳坤，李复.引入热力学温标的一种方案 J.大学物理，2006，25（02）：33-14.7 李椿，章立源，钱尚武.热学 M.北京：高等教育出版社，2015：8-9.8 李椿，章立源，钱尚武.热学 M.北京：高等教育出版社，2015：25-26.9 汪志诚.热力学 统计物理 M.北京：高等教育出版社，2011：270-271.10 李椿，章立源，钱尚武.热学 M.北京：高等教育出版社，2015：148-150.11 黄淑清，聂宜如，申先甲.热学教程 M.北京：高等教育出版社，2011：96-99.12 赵凯华，罗蔚茵.热学 M.北京：高等教育出版社，1998：186-188.13 马晓栋.关于热力学第二定律的定量讲述 J.乌鲁木齐成人教育学院学报，1999，7（01）：55-57.A Deep Understanding of the Concept of Temperature based on Teaching PracticeHAN Ji-ning，WEI Wei*，MA Xiao-dong（Xinjiang Key Laboratory for Luminescence Minerals and Optical Functional Materials，School of Physics and Electronic Engineering，Xinjiang Normal University，Urumqi，Xinjiang，830054，China）Abstract：Based on human cognitive process to the concept of temperature which is complex and abstract，this paper firstly elaborates the definition methods，argumentation process and deficiency of the qualitative and quantitative definitions of temperature in most thermology textbooks from the qualitative and quantitative perspectives and points out the shortcomings in the definition process.Secondly，the definitions methods and argumentation process of the concise qualitative definition and the concise quantitative definition of temperature in other literatures are given，and their completeness and superiority are discussed in detail.It is clearly pointed out that the law contents make up for the loopholes of the theorem.Finally，it makes necessary supplementary explanation to the details of the quantitative definition of temperature.Keywords：Qualitative definition；Quantitative definition；Thermodynamic temperature scale；Ideal gas temperature scaleA Simple Iteration Method for M-matrix Algebraic Riccati EquationsGUAN Jin-rui1，REN Fu-jiao1，SHAO Rong-xia2（1.School of Mathematics and Statistics，Taiyuan Normal University，Jinzhong，Shanxi，030619，China；2.School of Statistics and Data Science，Xinjiang University of Finance and Economics，Urumqi，Xinjiang，830012，China）Abstract：This paper studies numerical solution of M-matrix algebraic Riccati equations.When the coefficient matrix of the equation is a regular singular M-matrix，some existing numerical methods encounter certain difficulties.Therefore，a simple iteration method is proposed to solve the equation.This method only uses matrix multiplication in each iteration，with small computational complexity and easy implementation.Theoretical analysis and numerical experiments show that the method in this paper is feasible and effective under certain circumstances.Keywords：Algebraic Riccati equation；Regular M-matrix；Newton method；Iterative method（上接第5页）12