蔡兆克, 鮑 亮, 初 魯

(華東理工大學(xué)數(shù)學(xué)系,上海 200237)

預(yù)條件平方Smith法求解連續(xù)Lyapunov方程

蔡兆克, 鮑 亮, 初 魯

(華東理工大學(xué)數(shù)學(xué)系,上海 200237)

探討了如何數(shù)值求解連續(xù)時(shí)間的Lyapunov矩陣方程AX+XAT+BBT=0,給出了一種預(yù)條件的平方Smith算法,該算法首先利用交替方向隱式法即ADI法處理連續(xù)Lyapunov方程,構(gòu)造出含ADI參數(shù)的對(duì)稱Stein方程;然后利用平方Smith法迭代產(chǎn)生Krylov子空間中的低秩逼近形式。得到一些數(shù)值實(shí)驗(yàn),這些例子表明預(yù)條件平方Smith法是非常有效的。

連續(xù)Lyapunov方程; 預(yù)條件; 平方Smith算法;ADI;Krylov空間

Lyapunov方程是矩陣?yán)碚撗芯恐幸活惙浅Ｖ匾木仃嚪匠?在許多工程理論領(lǐng)域中都發(fā)揮著重要作用,如穩(wěn)定性分析[1]、模型降階[2-4]等問(wèn)題均需要求解該類方程。

本文考慮連續(xù)時(shí)間的Lyapunov矩陣方程:

AX+XAT+BBT=0

(1)

當(dāng)方程(1)的系數(shù)矩陣規(guī)模較小時(shí),由Bartels-Stewart、Golub、Hammarling等[1,3,5-6]提出的以Schur分解為基礎(chǔ)的直接法求解效果非常好,但由于該算法需要o(n3) 的運(yùn)算量和o(n2) 的存儲(chǔ)量,因此對(duì)于較大規(guī)模系數(shù)矩陣的方程,該方法不適用。而對(duì)于較大規(guī)模的矩陣方程,一般采用迭代法如Jacobi迭代法、ADI迭代法、Smith迭代法等[11-13,15-16]進(jìn)行求解,特別是基于Galerkin投影[9]的迭代法成為了近年來(lái)的主要研究方向,如全局GMRES、塊GMRES等Krylov子空間法[7-13]是求解此類大型矩陣方程的有效方法;該類迭代方法主要是利用Kronecker積形成一個(gè)n2×n2的大型線性系統(tǒng),結(jié)合Arnoldi過(guò)程構(gòu)造出Krylov子空間中的低秩逼近解,其中ADI預(yù)條件處理方式也成為加速數(shù)值解收斂的有效方式,如Jbilou在文獻(xiàn)[21]中利用低秩ADI預(yù)條件子將大型Lyapunov矩陣方程轉(zhuǎn)換成Stein方程進(jìn)而通過(guò)全局Arnoldi過(guò)程求解;此外,Sadkane在文獻(xiàn)[22]中給出了求解離散時(shí)間Lyapunov方程的平方Smith算法。本文成功地將其進(jìn)一步推廣到連續(xù)時(shí)間Lyapunov矩陣方程的數(shù)值求解過(guò)程中,并用數(shù)值實(shí)驗(yàn)例子說(shuō)明本文的預(yù)條件平方Smith算法相較于經(jīng)典Krylov子空間算法是非常有效的。

本文基于ADI預(yù)條件處理,介紹了一種平方Smith方法以求解連續(xù)時(shí)間Lyapunov矩陣方程在Krylov子空間中的低秩逼近解,并在這一求解過(guò)程的具體實(shí)施中提出了一種非常重要的改進(jìn)奇異值分解定理的處理方式;同時(shí)給出預(yù)條件平方Smith算法,以及誤差和殘量的估計(jì);進(jìn)而給出了數(shù)值實(shí)驗(yàn)例子及結(jié)果。

1 ADI預(yù)條件與平方Smith法

1.1 ADI預(yù)條件

Lyapunov矩陣方程(1) 的解Xi可由交替方向隱式迭代法即ADI迭代法[14-15]產(chǎn)生,即有線性系統(tǒng)

(2)

(3)

其中X0=0,{μ1,μ2,…,μi,…}∈-。顯然,聯(lián)立式(2)和式(3)有

(4)

(5)

AXAT-X+BBT=0

(6)

且參數(shù)μ產(chǎn)生于

(7)

本文中,始終假設(shè)矩陣A是穩(wěn)定矩陣(也即,矩陣A的特征值位于單位圓盤內(nèi)),則ρ(A)<1,那么方程(6)有如下形式的唯一解[19]:

(8)

1.2 平方Smith算法[22]

(9)

而由式(9)易得到:

(10)

因此,可以利用Krylov子空間k(A,B)=range(B,AB,…,A2k-1B)中的低秩矩陣Zk逼近。使用塊Arnoldi過(guò)程構(gòu)造Krylov子空間k(A,B)中的低秩矩陣Zk,塊Arnoldi算法如下:

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

(2)對(duì)j=1,…,2m作循環(huán),且2m?n

Wj=AQj

對(duì)i=1,…,j作循環(huán)

對(duì)Wj作QR分解:Wj=Qj+1Hj+1,j

1.3Krylov子空間中奇異值分解的應(yīng)用

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

那么,當(dāng)k=1,X1=(B,AB)(B,AB)T,由塊Arnoldi算法的結(jié)論知(B,AB)=(Q1R1,AQ1R1)=

根據(jù)該定理可有奇異值分解:

類似k=1情況的討論有:

在一般情況下,

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

2 誤差與殘量

2.1 概述

本節(jié)討論第k步迭代時(shí)的誤差Λk和殘量Γk。

2.2 誤差Λk

由命題1的證明過(guò)程可知,當(dāng)且僅當(dāng)ρa(bǔ)ρb<1時(shí)該命題成立,本文總是假設(shè)該條件成立;且當(dāng)k→∞時(shí),‖X-Xk‖→0;但是,如果ρa(bǔ)ρb非常接近于1或者常數(shù)Ca和Cb至少一個(gè)數(shù)值十分大,那么收斂速度就會(huì)大大變小。第1種情形時(shí),文獻(xiàn)中ADI迭代法[18]能夠極小化譜半徑;第2種情形時(shí),不存在低秩逼近解,與本文的假設(shè)條件不相符合,故不作考慮。

命題2 對(duì)任意的k≥0,成立

當(dāng)k≥1時(shí),由Zk表達(dá)式的理論推導(dǎo)過(guò)程可知:

那么,

(11)

其中

(12)

推論1 對(duì)任意的k≥0,成立

證明 由于

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

2.3 殘量Γk

關(guān)于殘量Γk=BBT+AXkAT-Xk,在命題3中對(duì)其進(jìn)行了估計(jì),并給出一種計(jì)算殘量的方式,該命題作為數(shù)值實(shí)驗(yàn)部分迭代停機(jī)準(zhǔn)則中殘量計(jì)算的理論基礎(chǔ)。

命題3 對(duì)任意的k≥0,殘量Γk滿足

證明 (1)根據(jù)殘量Γk的定義有

當(dāng)k≥1時(shí),

等式兩邊同時(shí)取范數(shù)即可。

2.4 預(yù)條件平方Smith算法

本節(jié)給出以命題3中殘量Γk的范數(shù)為迭代收斂判斷準(zhǔn)則的預(yù)條件平方Smith算法,算法如下:

算法2 PLRSS

(1)對(duì)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 數(shù)值實(shí)驗(yàn)

數(shù)值實(shí)驗(yàn)采用Matlab編程語(yǔ)言,針對(duì)Plrss算法、塊Krylov子空間法、全局Kryloy子空間法、Jacobi迭代法進(jìn)行編程實(shí)現(xiàn),討論分析不同規(guī)模系數(shù)矩陣情況下這些算法的收斂情況。在Windows8.1(64 bit)系統(tǒng)環(huán)境下的MATLAB R2014a中運(yùn)行,處理器為Intel Core2 @ 2.20 GHz,內(nèi)存為6.00 GB。

例1:在取上述隨機(jī)系數(shù)矩陣時(shí),若n=100,p=2,比較Plrss算法與塊Krylov子空間法、全局Krylov子空間法、Jacobi法的收斂情況,如圖1所示(本文所有數(shù)值實(shí)驗(yàn)結(jié)果圖表中:Plrss為Plrss算法,Block為塊Krylov子空間法,Global為全局Krylov子空間法,Jacobi為Jacobi迭代法),結(jié)果表明Plrss算法在求解連續(xù)時(shí)間Lyapunov方程時(shí)收斂效果更好。

圖1 當(dāng)n=100,p=2時(shí)各數(shù)值算法的結(jié)果比較Fig.1 Comparison of the algorithms with n=100,p=2

例2:為了進(jìn)一步說(shuō)明本文研究的Plrss算法適合于B為瘦長(zhǎng)形矩陣情況,取n=100,p=20與例1進(jìn)行比較。實(shí)驗(yàn)結(jié)果見(jiàn)圖2。結(jié)果表明當(dāng)n/p 小即弱化了p?n 這一條件時(shí),Plrss算法比塊Krylov子空間算法效果要差一些。

圖2 當(dāng)n=100,p=20時(shí),各數(shù)值算法結(jié)果比較Fig.2 Comparison of the algorithms with n=100,p=20

例3:為了進(jìn)一步說(shuō)明本文所研究的Plrss算法對(duì)于系數(shù)矩陣規(guī)模較大的連續(xù)時(shí)間Lyapunov方程也具有較好的求解效果,取n=500,p=2。實(shí)驗(yàn)結(jié)果如圖3所示。結(jié)果表明隨著系數(shù)矩陣規(guī)模的增大,相比塊Krylov子空間法、全局Krylov子空間法、Jacobi法,Plrss算法仍然具有更好的收斂效果。

4 結(jié)束語(yǔ)

本文提出了一種求解連續(xù)時(shí)間Lyapunov矩陣方程的預(yù)條件平方Smith算法,并在數(shù)值實(shí)驗(yàn)中,通過(guò)與塊Krylov子空間、全局Krylov子空間、Jacobi算法收斂結(jié)果的比較,直觀地說(shuō)明了在求解連續(xù)時(shí)間Lyapunov矩陣方程時(shí),Plrss算法具有很好的收斂效果。

圖3 當(dāng)n=500,p=2時(shí),各數(shù)值算法結(jié)果比較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

中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金

蔡兆克(1990-),男,安徽人,碩士生,研究方向?yàn)閿?shù)值代數(shù)及其應(yīng)用。E-mail:zhaoke_cai@163.com

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

O242.2

A