預條件平方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