段曉敏， 趙新玉， 孫華飛

(1.大連交通大學(xué) 材料科學(xué)與工程學(xué)院，遼寧，大連 116028;2.大連交通大學(xué) 理學(xué)院，遼寧，大連 116028;3.北京理工大學(xué) 數(shù)學(xué)學(xué)院，北京 100081)

?

Stein方程數(shù)值解的黎曼梯度算法

段曉敏1，2， 趙新玉1， 孫華飛3

(1.大連交通大學(xué) 材料科學(xué)與工程學(xué)院，遼寧，大連 116028;2.大連交通大學(xué) 理學(xué)院，遼寧，大連 116028;3.北京理工大學(xué) 數(shù)學(xué)學(xué)院，北京 100081)

基于正定矩陣流形的信息幾何結(jié)構(gòu)， 使用黎曼梯度算法來(lái)獲得Stein方程的數(shù)值解. 利用彎曲的黎曼流形上的測(cè)地距離作為算法的目標(biāo)函數(shù)，并將流形上的測(cè)地線作為算法的收斂路徑. 通過(guò)數(shù)值實(shí)驗(yàn)討論了算法的步長(zhǎng)和收斂速度的關(guān)系，從而得到算法的最優(yōu)步長(zhǎng).

Stein方程;黎曼梯度算法;數(shù)值模擬

考慮Stein矩陣方程

(1)

式中：A為任意的n×n實(shí)矩陣；Q為給定的n×n正定矩陣；X為待求的未知矩陣. 它在許多領(lǐng)域中起著重要的作用，如系統(tǒng)與控制理論、動(dòng)態(tài)規(guī)劃等科學(xué)與工程計(jì)算領(lǐng)域[1-3]. 在過(guò)去幾十年里，人們對(duì)這個(gè)方程的求解問(wèn)題進(jìn)行了各種研究，例如，Bartels-Stewart方法，Hessenberg-Schur方法，Schur方法和QR分解法[1]. 根據(jù)偏序集上的Kronecker內(nèi)積和固定點(diǎn)原理，文獻(xiàn)[4]討論了方程存在唯一正定矩陣解的充分條件，即，如果I-ATA是一個(gè)正定矩陣， 那么方程(1)存在唯一的解， 并且這個(gè)解是對(duì)稱(chēng)正定的矩陣， 其中I為n階單位矩陣.

信息幾何是在黎曼流形上采用現(xiàn)代微分幾何方法來(lái)研究統(tǒng)計(jì)學(xué)和信息領(lǐng)域問(wèn)題而提出的一套新的理論體系，已經(jīng)逐步應(yīng)用于神經(jīng)網(wǎng)絡(luò)、信號(hào)處理、系統(tǒng)與控制等領(lǐng)域[5]. 而將信息幾何理論應(yīng)用到研究矩陣方程的數(shù)值解法卻很少[6]. 本文將應(yīng)用信息幾何方法求解方程(1)的數(shù)值解問(wèn)題. 注意到解是一個(gè)對(duì)稱(chēng)正定矩陣， 并且所有的對(duì)稱(chēng)正定矩陣的集合構(gòu)成了一個(gè)流形，因此借助于流形的幾何結(jié)構(gòu)來(lái)研究這個(gè)矩陣的數(shù)值解更為方便. 為了滿(mǎn)足這樣的需要，提出解方程(1)的黎曼梯度下降法， 使用彎曲的黎曼流形上的測(cè)地距離作為目標(biāo)函數(shù)， 并且把測(cè)地線作為算法的收斂路徑. 最后， 用兩個(gè)數(shù)值例子演示黎曼梯度下降法的有效性.本文中， 用X>0表示X是對(duì)稱(chēng)正定的矩陣. 如果X-Y>0， 將會(huì)用X>Y表示.

1 預(yù)備知識(shí)

本節(jié)簡(jiǎn)要介紹黎曼流形上的幾何概念、黎曼梯度下降算法以及正定矩陣流形的幾何結(jié)構(gòu)，詳細(xì)內(nèi)容見(jiàn)文獻(xiàn)[7-8]．

1.1 黎曼流形上的幾何概念

給定一個(gè)正則函數(shù)f:M→R， 以向量V∈TXM為方向的黎曼梯度Xf代表f沿著V的變化率. 也就是說(shuō)， 給定任何光滑曲線φX，V:[0，1]→M， 使其滿(mǎn)足φX，V(0)=X∈M和X，V(0)=V， 那么黎曼梯度Xf是TXM上唯一滿(mǎn)足下式的切向量

一些光滑流形上的優(yōu)化問(wèn)題可以利用黎曼梯度下降法來(lái)解決[9]. 因此， 考慮流形M上的微分方程

在考慮黎曼流形M上的優(yōu)化算法時(shí)要考慮到流形的幾何特征. 具體地， 從初始值X0開(kāi)始， 設(shè)在迭代算法的步驟t∈N時(shí)刻， 算法的迭代點(diǎn)為Xt∈M. 利用文獻(xiàn)[9]的方法， 可以將Xt沿著測(cè)地線，向目標(biāo)函數(shù)J(Xt)的黎曼梯度的反方向移動(dòng)，則得到下一個(gè)點(diǎn)Xt+1. 也就是說(shuō)，如果設(shè)目標(biāo)函數(shù)J(Xt)沿著Xt方向的黎曼梯度為

那么黎曼梯度下降算法可以表示為

式中η為能夠使得算法收斂的任何合適的步長(zhǎng)．

1.2 正定矩陣流形的幾何結(jié)構(gòu)

所有的n×n對(duì)稱(chēng)正定矩陣可以構(gòu)成一個(gè)流形，記為P(n). 流形P(n)不是一般線性群GL(n)的李子群，而是GL(n)的一個(gè)子流形. 在流形P(n)上，可以定義不同的度量，這將賦予它不同的幾何結(jié)構(gòu). 設(shè)所有n×n對(duì)稱(chēng)矩陣組成的集合為S(n)，那么P(n)的切空間為S(n). 即任何對(duì)稱(chēng)矩陣的指數(shù)映射都是一個(gè)正定矩陣，并且任何正定矩陣的對(duì)數(shù)映射都是一個(gè)對(duì)稱(chēng)矩陣.設(shè)所有n×n實(shí)矩陣組成的集合為M(n)，把流形P(n)視為M(n)的一個(gè)子集，那么Frobenius內(nèi)積定義為

式中：B1，B2∈P(n).Frobenius內(nèi)積是歐氏的，不依賴(lài)于點(diǎn)的選擇. 這個(gè)內(nèi)積導(dǎo)致了P(n)是一個(gè)平坦的黎曼流形. 而且測(cè)地線是直線，還有一定的限制條件. 在流形P(n)上改進(jìn)的Frobenius內(nèi)積在B處的定義為

(2)

式中X，Y∈S(n). 于是， 帶有黎曼度量(2)的流形P(n)成為了黎曼流形， 而且， 它使得黎曼流形P(n)是彎曲的.

2 黎曼梯度算法

如果A滿(mǎn)足

(3)

那么，方程(1)的解是唯一正定的[2]. 設(shè)L(X):=Q+ATXA，顯然，對(duì)任何的X∈P(n)，L(X)是一個(gè)正定矩陣. 本文目標(biāo)是尋找一個(gè)矩陣X∈P(n)使得L(X)與給定的矩陣Q越接近越好.注意到黎曼流形P(n)是測(cè)地完備的， 因此用Q和L(X)的測(cè)地距離作為目標(biāo)函數(shù)J(X)

那么，方程(1)的精確解可以被定義為

(4)

為了求得J(X)的最小值點(diǎn)， 利用微分方程， 顯然最小值在X*處實(shí)現(xiàn)， 它滿(mǎn)足

(5)

那么黎曼梯度下降算法可以表示為

(6)

式中η為能夠使得算法(6)收斂的任何合適的步長(zhǎng)．

為了使用算法式(5)(6)，首先需要計(jì)算目標(biāo)函數(shù)J(Xt)的黎曼梯度. 根據(jù)黎曼梯度的定義， 獲得如下定理．

定理1 目標(biāo)函數(shù)J(X)在X處的黎曼梯度為

證明 計(jì)算實(shí)值函數(shù)

根據(jù)文獻(xiàn)[4]的性質(zhì)2.1以及一些矩陣跡的性質(zhì)有

事實(shí)上， 因?yàn)?/p>

對(duì)?A∈GL(n)和?B∈P(n)成立， 那么XJ(X)可以被重寫(xiě)為

證畢.

綜上， 黎曼梯度下降算法(RGA)為

式中初值X0∈P(n). 為了使RGA算法具有快速的收斂速度， 步長(zhǎng)η可以通過(guò)實(shí)驗(yàn)得到， 細(xì)節(jié)將在模擬實(shí)例中給出．

3 例 證

本章給出兩個(gè)例子說(shuō)明黎曼梯度算法(RGA)的計(jì)算有效性. 模擬中的誤差設(shè)ρ(L(X))=Q<10-15，其中ρ(·)表示矩陣的譜半徑．

例1 首先考慮當(dāng)n=5時(shí)，Stein方程

Q=X-ATXA

的數(shù)值解， 式中：

Q=

A=

它們滿(mǎn)足前提條件式(3)．

在下面的模擬中， 使用RGA算法取5階單位陣作為初值X0. 將步長(zhǎng)大小和收斂速度的關(guān)系用圖1表示出來(lái). 可以發(fā)現(xiàn)， 當(dāng)步長(zhǎng)從0變化到0.54時(shí)， 迭代次數(shù)逐漸減少， 步長(zhǎng)比0.54大時(shí)算法收斂速度又變慢了. 因此， 通過(guò)實(shí)驗(yàn)， 臨界值0.54被選為RGA算法的步長(zhǎng)．

經(jīng)過(guò)Matlab程序計(jì)算， 需要32次迭代得到方程的數(shù)值解為

針對(duì)這組數(shù)據(jù)算法的收斂速度請(qǐng)見(jiàn)圖2．

例2 考慮10階的Stein矩陣方程， 其中Q∈P(10)，A滿(mǎn)足條件(3). 類(lèi)似上節(jié)的計(jì)算過(guò)程， 選擇最優(yōu)步長(zhǎng)為0.52，模擬結(jié)果顯示RGA算法需要11次迭代， 計(jì)算時(shí)間為0.002 5 s.

4 結(jié) 論

借助于正定矩陣流形的信息幾何結(jié)構(gòu)，利用黎曼梯度算法求解Stein方程的數(shù)值解.這個(gè)算法使用彎曲的黎曼流形上的測(cè)地距離作為目標(biāo)函數(shù)， 并將流形上的測(cè)地線作為算法的收斂路徑.為了給出最快的收斂速度， 優(yōu)化的步長(zhǎng)被給出. 最后通過(guò)兩個(gè)數(shù)值例子表明黎曼梯度下降法求解Stein方程數(shù)值解的有效性．

[1] Gajic Z， Qureshi M T J. Lyapunov matrix equation in system stability and control[M]. [S.l.]: Courier Corporation， 2008.

[2] Heinen J. A technique for solving the extended discrete Lyapunov matrix equation[J]. IEEE Transactions on Automatic Control， 1972，17(2):156-157.

[3] Ran A C M， Reurings M C B. The symmetric linear matrix equation[J]. Electronic Journal of Linear Algebra， 2002，9(16):93-107.

[4] Ran A C M， Reurings M C B. A fixed point theorem in partially ordered sets andsome applications to matrix equations[J]. Proceedings of the American Mathematical Society， 2003，132:1435-1443．

[5] Amari S， Nagaoka H. Methods of information geometry[M]. Oxford: Oxford University Press， 2000．

[6] Duan X， Sun H， Peng L， et al. A natural gradient algorithm for the solution of discrete algebraic Lyapunov equations based on the geodesic distance[J]. Applied Mathematics and Computation， 2013，219:9899-9905．

[7] Moakher M. A differential geometric approach to the geometric mean of symmetric positive-definite matrices[J]. SIAM Journal on Matrix Analysis and Applications， 2005，26(3):735-747.

[8] Spivak M. A comprehensive introduction to differential geometry[M]. 3rd ed. Berkeley: Publish or Perish， 1999.

[9] Udriste C. Convex functions and optimization methods on Riemannian manifolds[M]. Berlin: Springer Science & Business Media， 1994.

(責(zé)任編輯：李兵)

Riemannian Gradient Algorithm for the Numerical Solution of Stein Equations

DUAN Xiao-min1，2， ZHAO Xin-yu1， SUN Hua-fei3

(1.School of Science， Dalian Jiaotong University， Dalian，Liaoning 116028， China;2.School of Materials Science and Engineering， Dalian Jiaotong University， Dalian，Liaoning 116028， China;3.Mathematical School， Beijing Institute of Technology， Beijing 100081， China)

A Riemannian gradient algorithm based on information geometric structures of a manifold consisting of all symmetric positive-definite matrices was proposed to calculate the numerical solution of Stein equations. In this algorithm， the geodesic distance on the curved Riemannian manifoldis taken as an objective function and the geodesic curve was treated as the convergence path. Also the optimal variable step-sizes corresponding to the minimum value of the objective function were provided in order to improve the convergence speed.

Stein equation; Riemannian gradient algorithm;simulation

2015-07-03

國(guó)家自然科學(xué)基金資助項(xiàng)目(61401058， 61179031);中國(guó)博士后科學(xué)基金面上資助項(xiàng)目(2015M581323)

段曉敏(1979—)，女，博士，講師，E-mail：dxmhope@djtu.edu.cn；趙新玉(1979-)，男，博士，副教授，E-mail：xyz@djtu.edu.cn.

孫華飛(1958—)，男，教授，博士生導(dǎo)師，E-mail：huafeisun@bit.edu.cn.

O 186

A

1001-0645(2016)02-0201-04

10.15918/j.tbit1001-0645.2016.02.018