预条件平方Smith法求解连续Lyapunov方程

蔡兆克, 鲍 亮, 初 鲁

(华东理工大学数学系,上海 200237)

预条件平方Smith法求解连续Lyapunov方程

蔡兆克, 鲍 亮, 初 鲁

(华东理工大学数学系,上海 200237)

探讨了如何数值求解连续时间的Lyapunov矩阵方程AX+XAT+BBT=0,给出了一种预条件的平方Smith算法,该算法首先利用交替方向隐式法即ADI法处理连续Lyapunov方程,构造出含ADI参数的对称Stein方程;然后利用平方Smith法迭代产生Krylov子空间中的低秩逼近形式。得到一些数值实验,这些例子表明预条件平方Smith法是非常有效的。

连续Lyapunov方程; 预条件; 平方Smith算法;ADI;Krylov空间

Lyapunov方程是矩阵理论研究中一类非常重要的矩阵方程,在许多工程理论领域中都发挥着重要作用,如稳定性分析[1]、模型降阶[2-4]等问题均需要求解该类方程。

本文考虑连续时间的Lyapunov矩阵方程:

AX+XAT+BBT=0

(1)

当方程(1)的系数矩阵规模较小时,由Bartels-Stewart、Golub、Hammarling等[1,3,5-6]提出的以Schur分解为基础的直接法求解效果非常好,但由于该算法需要o(n3) 的运算量和o(n2) 的存储量,因此对于较大规模系数矩阵的方程,该方法不适用。而对于较大规模的矩阵方程,一般采用迭代法如Jacobi迭代法、ADI迭代法、Smith迭代法等[11-13,15-16]进行求解,特别是基于Galerkin投影[9]的迭代法成为了近年来的主要研究方向,如全局GMRES、块GMRES等Krylov子空间法[7-13]是求解此类大型矩阵方程的有效方法;该类迭代方法主要是利用Kronecker积形成一个n2×n2的大型线性系统,结合Arnoldi过程构造出Krylov子空间中的低秩逼近解,其中ADI预条件处理方式也成为加速数值解收敛的有效方式,如Jbilou在文献[21]中利用低秩ADI预条件子将大型Lyapunov矩阵方程转换成Stein方程进而通过全局Arnoldi过程求解;此外,Sadkane在文献[22]中给出了求解离散时间Lyapunov方程的平方Smith算法。本文成功地将其进一步推广到连续时间Lyapunov矩阵方程的数值求解过程中,并用数值实验例子说明本文的预条件平方Smith算法相较于经典Krylov子空间算法是非常有效的。

本文基于ADI预条件处理,介绍了一种平方Smith方法以求解连续时间Lyapunov矩阵方程在Krylov子空间中的低秩逼近解,并在这一求解过程的具体实施中提出了一种非常重要的改进奇异值分解定理的处理方式;同时给出预条件平方Smith算法,以及误差和残量的估计;进而给出了数值实验例子及结果。

1 ADI预条件与平方Smith法

1.1 ADI预条件

Lyapunov矩阵方程(1) 的解Xi可由交替方向隐式迭代法即ADI迭代法[14-15]产生,即有线性系统

(2)

(3)

其中X0=0,{μ1,μ2,…,μi,…}∈-。显然,联立式(2)和式(3)有

(4)

(5)

AXAT-X+BBT=0

(6)

且参数μ产生于

(7)

本文中,始终假设矩阵A是稳定矩阵(也即,矩阵A的特征值位于单位圆盘内),则ρ(A)<1,那么方程(6)有如下形式的唯一解[19]:

(8)

1.2 平方Smith算法[22]

(9)

而由式(9)易得到:

(10)

因此,可以利用Krylov子空间k(A,B)=range(B,AB,…,A2k-1B)中的低秩矩阵Zk逼近。使用块Arnoldi过程构造Krylov子空间k(A,B)中的低秩矩阵Zk,块Arnoldi算法如下:

算法1:(1)对B作QR分解:B=Q1R1,其中Q1∈n×p,R1∈p×p

(2)对j=1,…,2m作循环,且2m≪n

Wj=AQj

对i=1,…,j作循环

对Wj作QR分解:Wj=Qj+1Hj+1,j

1.3Krylov子空间中奇异值分解的应用

首先,根据平方Smith迭代公式(9)和(10)知:X0=BBT

那么,当k=1,X1=(B,AB)(B,AB)T,由块Arnoldi算法的结论知(B,AB)=(Q1R1,AQ1R1)=

根据该定理可有奇异值分解:

类似k=1情况的讨论有:

在一般情况下,

Xk≈(Zk-1,Ak-1Zk-1)(Zk-1,Ak-1Zk-1)T;同样也有:

2 误差与残量

2.1 概述

本节讨论第k步迭代时的误差Λk和残量Γk。

2.2 误差Λk

由命题1的证明过程可知,当且仅当ρaρb<1时该命题成立,本文总是假设该条件成立;且当k→∞时,‖X-Xk‖→0;但是,如果ρaρb非常接近于1或者常数Ca和Cb至少一个数值十分大,那么收敛速度就会大大变小。第1种情形时,文献中ADI迭代法[18]能够极小化谱半径;第2种情形时,不存在低秩逼近解,与本文的假设条件不相符合,故不作考虑。

命题2 对任意的k≥0,成立

当k≥1时,由Zk表达式的理论推导过程可知:

那么,

(11)

其中

(12)

推论1 对任意的k≥0,成立

证明 由于

而又由命题1和命题2,那么有

2.3 残量Γk

关于残量Γk=BBT+AXkAT-Xk,在命题3中对其进行了估计,并给出一种计算残量的方式,该命题作为数值实验部分迭代停机准则中残量计算的理论基础。

命题3 对任意的k≥0,残量Γk满足

证明 (1)根据残量Γk的定义有

当k≥1时,

等式两边同时取范数即可。

2.4 预条件平方Smith算法

本节给出以命题3中残量Γk的范数为迭代收敛判断准则的预条件平方Smith算法,算法如下:

算法2 PLRSS

(1)对B作QR分解:B=Q1R1;

(2)令 U=I,V=I,S=R1,Zs=[],p=dim(Q1),j=0,iter=0,rst=0;

(4)若j=2k,则

否则,约化奇异值分解消除中小于tolsvd的奇异值:

Z=jUS;

iter=iter+1;

(6)若dim(j+1)>mmax,则

(Zs,Z)列正交化后,并赋值Zs=(Zs,Z);

约化奇异值分解:

3 数值实验

数值实验采用Matlab编程语言,针对Plrss算法、块Krylov子空间法、全局Kryloy子空间法、Jacobi迭代法进行编程实现,讨论分析不同规模系数矩阵情况下这些算法的收敛情况。在Windows8.1(64 bit)系统环境下的MATLAB R2014a中运行,处理器为Intel Core2 @ 2.20 GHz,内存为6.00 GB。

例1:在取上述随机系数矩阵时,若n=100,p=2,比较Plrss算法与块Krylov子空间法、全局Krylov子空间法、Jacobi法的收敛情况,如图1所示(本文所有数值实验结果图表中:Plrss为Plrss算法,Block为块Krylov子空间法,Global为全局Krylov子空间法,Jacobi为Jacobi迭代法),结果表明Plrss算法在求解连续时间Lyapunov方程时收敛效果更好。

图1 当n=100,p=2时各数值算法的结果比较Fig.1 Comparison of the algorithms with n=100,p=2

例2:为了进一步说明本文研究的Plrss算法适合于B为瘦长形矩阵情况,取n=100,p=20与例1进行比较。实验结果见图2。结果表明当n/p 小即弱化了p≪n 这一条件时,Plrss算法比块Krylov子空间算法效果要差一些。

图2 当n=100,p=20时,各数值算法结果比较Fig.2 Comparison of the algorithms with n=100,p=20

例3:为了进一步说明本文所研究的Plrss算法对于系数矩阵规模较大的连续时间Lyapunov方程也具有较好的求解效果,取n=500,p=2。实验结果如图3所示。结果表明随着系数矩阵规模的增大,相比块Krylov子空间法、全局Krylov子空间法、Jacobi法,Plrss算法仍然具有更好的收敛效果。

4 结束语

本文提出了一种求解连续时间Lyapunov矩阵方程的预条件平方Smith算法,并在数值实验中,通过与块Krylov子空间、全局Krylov子空间、Jacobi算法收敛结果的比较,直观地说明了在求解连续时间Lyapunov矩阵方程时,Plrss算法具有很好的收敛效果。

图3 当n=500,p=2时,各数值算法结果比较Fig.3 Comparison of the algorithms with n=500,p=2

[1] LASALLE J P,LEFSCHETZ S.Stability by Lyapunov’s Direct Method:With Applications [M].New York:Academic Press,1961.

[2] MOORE B C.Principal component analysis in linear systems:Controllability,observability,and model reduction [J].IEEE Transactions on Automatic Control,1981,26(1):17-32.

[3] SAFONOV M G,CHIANG R Y.A Schur method for balanced-truncation model reduction [J].IEEE Transactions on Automatic Control,1989,34(7):729-733.

[4] TOMBS M S,POSTLETHWAITE I.Truncated balanced realization of a stable non-minimal state-space system [J].International Journal of Control,1987,46(4):1319-1330.

[5] BARTELS R H,STEWART G W.Solution of the matrix equationAX+XB=C[J].Communications of the ACM,1972,15(9):820-826.

[6] GOLUB G H,NASH S,VAN LOAN C.A Hessenberg-Schur method for the problemAX+XB=C[J].IEEE Transactions on Automatic Control,1979,24(6):909-913.

[7] HU D Y,REICHEL L.Krylov-subspace methods for the Sylvester equation [J].Linear Algebra and Its Applications,1992,172:283-313.

[8] JAIMOUKHA I M,KASENALLY E M.Krylov subspace methods for solving large Lyapunov equations [J].SIAM Journal on Numerical Analysis,1994,31(1):227-251.

[9] JBILOU K,RIQUET A J.Projection methods for large Lyapunov matrix equations [J].Linear Algebra and Its Applications,2006,415(2):344-358.

[10] SAAD Y.Numerical solution of large Lyapunov equations [J].Signal Processing,Scattering & Operator Theory & Numerical Methods,Proc Mtns,2010,47(2):503-511.

[11] SIMONCINI V.A new iterative method for solving large-scale Lyapunov matrix equations [J].SIAM Journal on Scientific Computing,2007,29(3):1268-1288.

[12] LU A,Wachspress E L.Solution of Lyapunov equations by alternating direction implicit iteration [J].Computers & Mathematics with Applications,1991,21(9):43-58.

[13] PENZL T.A cyclic low-rank Smith method for large sparse Lyapunov equations [J].SIAM Journal on Scientific Computing,1999,21(4):1401-1418.

[14] YOUNG D.The numerical solution of elliptic and parabolic partial differential equations [J].Survey of Numerical Analysis,1961:380-438.

[15] WACHSPRESS E L.Iterative solution of the Lyapunov matrix equation [J].Applied Mathematics Letters,1988,1(1):87-90.

[16] YOUSEF S.Iterative Methods for Sparse Linear Systems [M].Philadelphia:SIAM,2003.

[17] HOM R A,JOHNSON C R.Topics in Matrix Analysis [M].New York:Cambridge UP,1991.

[18] BENNER P,KHOURY G E,SADKANE M.On the squared Smith method for large-scale Stein equations [J].Numerical Linear Algebra with Applications,2014,21(5):645-665.

[19] LANCASTER P,RODMAN L.Algebraic Riccati Equations [M].Oxford:Clarendon Press,1995.

[20] GOLUB G H,VAN LOAN C F.Matrix Computations [M].Baltimore:JHU Press,2012.

[21] JBILOU K.ADI preconditioned Krylov methods for large Lyapunov matrix equations [J].Linear Algebra and Its Applications,2010,432(10):2473-2485.

[22] SADKANE M.A low-rank Krylov squared Smith method for large-scale discrete-time Lyapunov equations [J].Linear Algebra and Its Applications,2012,436(8):2807-2827.

A Preconditioned Squared Smith Method forContinuous-Time Lyapunov Equations

CAI Zhao-ke, BAO Liang, CHU Lu

(Department of Mathematics,East China University of Science and Technology,Shanghai 200237,China)

This paper proposes a preconditioned squared Smith algorithm to solve the continuous-time Lyapunov matrix equations AX+XAT+BBT=0numerically.Themethodfirstusesthealternatingdirectionalimplicit(ADI)methodandtransformstheoriginalequationstotheequivalentsymmetricSteinmatrixequationswithsomeADIparameters.ThenweadoptthesquaredSmithalgorithmtoseeksolutionsoftheSteinequationsbygeneratingthesquaredSmithiterationsinsomelow-rankformswiththeKrylovsubspaces.Andwegivesomenumericalexperimentstoshowtheefficiencyofthisalgorithmfinally.

continuous-time Lyapunov equations; preconditioned; squared Smith algorithm; ADI; Krylov subspace

1006-3080(2016)06-0881-06

10.14135/j.cnki.1006-3080.2016.06.021

2016-02-01

中央高校基本科研业务费专项资金

蔡兆克(1990-),男,安徽人,硕士生,研究方向为数值代数及其应用。E-mail:zhaoke_cai@163.com

鲍 亮,E-mail:lbao@ecust.edu.cn

O242.2

A