邵新慧, 彭程
东北大学 理学院, 辽宁 沈阳 110819
收稿日期：2015-12-30
基金项目：国家自然科学基金资助项目(11071033)。
作者简介：邵新慧(1970-)，女，山东青岛人，东北大学副教授，博士。

摘要：针对Sylvester矩阵方程给出了一种基于梯度的迭代解法.通过引入一个松弛参数和应用层次识别原理, 构建了一种新型的迭代方法求解一类Sylvester矩阵方程.收敛分析表明, 在一定的假设条件下对于任意初始值, 迭代解都收敛到精确解.数值算例也表明了所给方法的有效性和优越性.
关键词：Sylvester矩阵方程迭代解法梯度收敛松弛参数
Iterative Solutions to Sylvester Matrix Equations
SHAO Xin-hui, PENG Cheng
School of Sciences, Northeastern University, Shenyang 110819, China
Corresponding author: SHAO Xin-hui, E-mail: xinhui1002@126.com
Abstract: The gradient-based iterative solutions to Sylvester matrix equations are given. By introducing a relaxation parameter and applying the hierarchical identification principle, an iterative algorithm is constructed to solve Sylvester matrix equations. Convergence analysis indicates that the iterative solutions converge with the exact solutions to any initial value under certain assumptions. Numerical examples are given to testify the efficiency of the proposed method.
Key Words: Sylvester matrix equationiterative solutiongradientconvergencerelaxation parameter
Sylvester矩阵方程在很多领域中都有涉及, 例如矩形域上的椭圆边值问题离散后会得到一个Sylvester形式的矩形方程；在常微分方程定性理论研究及数值求解的隐式Runge-Kwutta方法与块方法中常会碰到这类方程的求解问题；在线性统计领域中也会有Sylvester矩阵方程；在控制理论及应用中, 极点配置、观测器及构造Sylvester函数中都涉及上述方程的求解.除此之外, 特殊形式的Sylvester方程, 即Lyapunov方程, 在实际应用中, 特别是在连续和离散稳定性分析中具有重要意义, 在系统与控制领域中经常会遇到.与此同时, Sylvester矩阵方程在图像恢复、模型降阶、神经网络、信号处理等领域中也有涉及.Sylvester矩阵方程的求解问题通常包括精确解和数值解的研究.尽管在系统稳定性分析等众多应用中, 精确解的研究非常重要, 但是, 当系统的参数未知时, 不可能得到精确解, 此时数值解的研究就非常必要.近些年, 国内外很多学者对各种类型的Sylvester矩阵方程的数值解进行研究, 得到了一些有效的数值求解算法^[1-15].
本文从Sylvester矩阵方程AXB=F的迭代解法入手, 先给出了它的迭代解的收敛条件, 然后介绍了计算Sylvester矩阵方程AXB+CXD=F数值解的分层迭代法，以及解决Sylvester矩阵方程AX+XB=C的松弛迭代法.
1 矩阵方程AXB+CXD=F的松弛迭代法1.1 相关引理对于任意的两个整数i < j, I[i, j]表示集合{i, i+1, …, j}, 对任意方阵A, , tr (A), ‖A‖分别表示矩阵共轭, 矩阵的迹和矩阵的F范数.对X=[x₁x₂…x_n]∈C^m×n, 设vec (X)=[x₁^Tx₂^T…x_n^T]^T, ‖A‖=‖vec (A)‖.对于任意矩阵M和N, M∈C^r×m, N∈C^n×s, MXN表示克罗内克积, 有vec(MXN)=(N^T?M)vec (X).
引理1^[2] ??对于AXB=F, 其中A∈C^p×m, B∈C^n×q, F∈C^p×q是已知矩阵, X∈C^m×n是未知矩阵.若满足或, 且求解方程AXB=F的迭代算法表示为X(k+1)=X(k)+μA^T(F-AX(k)B)B^T.若方程AXB=F有唯一解X_*, 则迭代解X(k)收敛到唯一解, 即.
引理2^[2] ??若A∈C^p×m是列满秩矩阵, B∈C^n×q是行满秩矩阵(p≥m, n≤q), 则方程AXB=F有唯一解, 且X=(A^TA)^-1A^TFB^T(BB^T)^-1, 所对应的齐次方程AXB=O有唯一解X=O.
1.2 迭代格式的建立考虑Sylvester矩阵方程

(1)

其中:A, C∈C^m×r；B, D∈C^s×n；F∈C^m×n是已知的常数矩阵；X∈C^r×s是待求矩阵.
假设X_*表示方程(1) 的精确解, X(k)表示X_*的第k步迭代的近似解.应用分层辨识原理^{[2-5, 7-8]}, 定义矩阵

(2)

(3)

引入一个松弛因子, 建立如下松弛递归格式:

(4)

(5)

其中:μ是收敛因子；ω是一个松弛因子, 满足0 < ω < 1.
将方程(2) 和(3) 代入到方程(4) 和(5) 中, 得到如下格式:

(6)

(7)

由于上述两个等式的右端包括未知矩阵X, 算法无法实现.根据分层辨识原理^{[2-5, 7-8]}, 方程(6) 和(7) 中X可以被它的迭代近似值X(k-1) 代替, 得到

(8)

(9)

由于求的是近似解X(k)而不是X₁(k)和X₂(k), 遵循平衡策略形成第k步迭代的近似解:

(10)

(11)

(12)

于是迭代算法改成如下格式:

(13)

1.3 收敛性分析定理1 ??设方程AXB+CXD=F有唯一解X_*, 其中A, C∈C^m×r, B, D∈C^s×n, F∈C^m×n, X∈C^r×s.设λ_max(AA^T)=λ₁, λ_max(BB^T)=λ₂, λ_max(CC^T)=λ₃, λ_max(DD^T)=λ₄.若收敛因子μ满足, 则对于任意初值X(0), 算法(13) 给出的迭代解收敛于唯一解X_*, 即.
证明??定义误差矩阵

(14)

(15)

(16)

设

(17)

将式(15) 代入到式(10), (11)：
将式(16), (17) 代入到上式, 可得X₁(k)=X(k-1)-ωμA^T(ξ(k-1)+η(k-1))B^T,
X₂(k)=X(k-1)-(1-ω)μC^T(ξ(k-1)+η(k-1))D^T.
由于|X|²=tr (X^TX), tr (AB)=tr (BA), tr (A)=tr (A^T), 得
.设λ_max(AA^T)=λ₁, λ_max(BB^T)=λ₂, λ_max(CC^T)=λ₃, λ_max(DD^T)=λ₄.可得‖A^T(ξ(k-1)+η(k-1))B^T‖²≤λ₁λ₂‖(k-1)+η(k-1)², 则 tr [(ξ(k-1)+η(k-1))ξ^T(k-1)+(ξ(k-1)+η(k-1))^Tξ(k-1)].则+ω²μ²λ₁λ₂‖ξ(k-1)+η(k-1)‖²-μωtr [(ξ(k-1)+η(k-1))ξ^T(k-1)+(ξ(k-1)+η(k-1))^Tξ(k-1)].
同样(1-ω)²μ²λ₃λ₄‖ξ(k-1)+η(k-1)‖²-μ(1-ω)tr [(ξ(k-1)+η(k-1))η^T(k-1)+(ξ(k-1)+η(k-1))^Tη(k-1)].
将等式(12) 两边同时取范数, 有
可得.若μ满足, 有, 则‖ξ(k)+η(k)‖²→0.进而得ξ(k)+η(k)→0, 结合引理1, .
2 数值模拟例??考虑方程AXB+CXD=F, 其中:设方程AXB+CXD=F的精确解为
利用算法(10)~(12), 取初始矩阵, 定义相对误差.根据定理1计算得到, λ_max(AA^T)=λ₁=2, λ_max(BB^T)=λ₂=2, λ_max(CC^T)=λ₃=5, λ_max(DD^T)=λ₄=2, 故有
通过大量实验, 当ω=0.49及μ==0.283 29时, 有更小的相对误差和较快的迭代速度, 如图 1~图 3所示.
图 1(Fig. 1)

图 1 当ω=0.49时算法(10)~(12) 的收敛性Fig.1 Convergence performance of algorithm (10)~(12) for different μ when ω=0.49

图 2(Fig. 2)

图 2 ω=0.49时算法(13) 的收敛性Fig.2 Convergence performance of algorithm (13) for different μ when ω=0.49

图 3(Fig. 3)

图 3 当μ=100/353时算法(13) 的收敛性Fig.3 Convergence performance of algorithm (13) for different ω when μ=100/353

当ω=0.49, μ=100/353时, 表 1给出用算法(13) 和文献[2]中算法求得解与相对误差δ₁.
表 1(Table 1)

表 1 本文算法与文献[2]中算法比较Table 1 Comparison of algorithm (13) and algorithm in ref. [2] k x11 x12 x21 x22 δ/% δ1/% 1 2.407 0 1.132 7 2.407 0 5.238 7 28.34 28.57 2 0.944 9 1.824 4 2.749 1 4.603 7 8.04 8.16 3 1.111 7 1.927 4 2.948 0 5.012 1 2.28 2.33 4 0.997 6 1.984 6 2.979 0 4.968 5 0.65 0.66 5 1.008 9 1.993 9 2.995 5 5.000 4 0.18 0.19 6 1.000 0 1.998 7 2.998 2 4.997 5 0.052 574 0.054 399 7 1.000 7 1.999 5 2.999 6 5.000 0 0.014 986 0.015 543 8 1.000 0 1.999 9 2.999 9 4.999 8 0.004 275 8 0.004 440 7 9 1.000 1 2.000 0 3.000 0 5.000 0 0.001 221 1 0.001 268 8 10 1.000 0 2.000 0 3.000 0 5.000 0 0.000 349 02 0.000 362 51 11 1.000 0 2.000 0 3.000 0 5.000 0 0.000 099 848 0.000 103 57

表 1 本文算法与文献[2]中算法比较 Table 1 Comparison of algorithm (13) and algorithm in ref. [2]

通过图表可以发现, 算法(13) 明显优于文献[2]中的算法.
3 结论本文讨论了实数域下Sylvester矩阵方程的迭代解法, 其主要思想是基于梯度迭代法与分层迭代原理, 通过引入松弛因子, 建立新的迭代格式, 从而获得Sylvester矩阵方程的数值解.还给出了相应的收敛性分析, 进行了理论证明, 得出如下结论:在特定条件下, 对于任意初始值, 方程的迭代解均收敛于精确解.最后通过数值算例说明了所给新算法的有效性、可行性以及优越性.从讨论可以看出, 同等情况下, 带有松弛因子的迭代算法与之前的方法相比, 可以带来较小的迭代误差与较快的迭代速度.另外, 本文给出的迭代算法更具有普遍性, 可以解多种Sylvester矩阵方程.
参考文献

[1]	Evans D J, Martins M M. AOR method for AX-XB=C[J].International Journal of Computer Mathematics, 1994, 52(2): 75–82.
[2]	Ding F, Chen T. Hierarchical least squares identification methods for multivariable systems[J].IEEE Transactions on Automatic Control, 2005, 50(3): 397–402.DOI:10.1109/TAC.2005.843856
[3]	徐树方. 控制论中的矩阵计算[M]. 北京: 高等教育出版社, 2011: 27-29. ( Xu Shu-fang. Matrix computation in cybernetics[M]. Beijing: Higher Education Press, 2011: 27-29.)
[4]	Ding F, Chen T. Gradient based iterative algorithms for solving a class of matrix equations[J].IEEE Transactions on Automatic Control, 2005, 50(8): 1216–1221.DOI:10.1109/TAC.2005.852558
[5]	Ding F, Liu P X, Ding J. Iterative solutions of the generalized Sylvester matrix equations by using the hierarchical identification principle[J].Applied Mathematics and Computation, 2008, 197(1): 41–50.DOI:10.1016/j.amc.2007.07.040
[6]	Xie L, Ding J, Ding F. Gradient based iterative solutions for general linear matrix equations[J].Computers & Mathematics with Applications, 2009, 58(7): 1441–1448.
[7]	Ding J, Liu Y J, Ding F. Iterative solutions to matrix equations of form AXB=F[J].Computers & Mathematics with Applications, 2010, 59(11): 3500–3507.
[8]	Kilicman A, Zhour Z A. Vector least-squares solutions for coupled singular matrix equations[J].Journal of Computational and Applied Mathematics, 2007, 206(2): 1051–1069.DOI:10.1016/j.cam.2006.09.009
[9]	Wu A G, Lu L L, Duan G R. Iterative algorithms for solving a class of complex conjugate and transpose matrix equations[J].Applied Mathematics and Computation, 2011, 217(21): 8343–8353.DOI:10.1016/j.amc.2011.02.113
[10]	Song C Q, Chen G L. An efficient algorithm for solving extended Sylvester conjugate transpose matrix equations[J].Arab Journal of Mathematical Sciences, 2011, 17(21): 115–134.
[11]	Song C Q, Chen G L, Zhao L L. Iterative solutions to coupled Sylvester transpose matrix equations[J].Applied Mathematical Modelling, 2011, 35(10): 4675–4683.DOI:10.1016/j.apm.2011.03.038
[12]	Niu Q, Wang X, Lu L Z. A relaxed gradient based algorithm for solving Sylvester equations[J].Asian Journal of Control, 2011, 13(3): 461–464.DOI:10.1002/asjc.v13.3
[13]	Dehghan M, Hajarian M. An iterative algorithm for the reflexive solutions of the generalized coupled Sylvester matrix equations and its optimal approximation[J].Applied Mathematics and Computation, 2008, 202(2): 571–588.DOI:10.1016/j.amc.2008.02.035
[14]	Tian Z L, Gu C Q. A numerical algorithm for Lyapunov equations[J].Applied Mathematics and Computation, 2008, 202(1): 44–53.DOI:10.1016/j.amc.2007.12.057
[15]	Hou J J, Peng Z Y, Zhang X L. An iterative method for the least squares symmetric solution of matrix equation AXB=C[J].Numerical Algorithms, 2006, 42(2): 181–192.DOI:10.1007/s11075-006-9037-3