非線性矩陣方程數(shù)值解的廣義哈密頓算法

段曉敏, 趙新玉, 孫華飛, 紀(jì)雪婷

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

考慮一類廣義非線性矩陣方程

(1)

注意到方程(1)的解是對(duì)稱正定矩陣, 而對(duì)稱正定矩陣全體可以構(gòu)建成一個(gè)流形，所以基于流形的幾何結(jié)構(gòu)提出求解方程(1)的算法具有很大的優(yōu)勢(shì). 著名的學(xué)者Fiori[8]提出了適用于流形上某些問(wèn)題的廣義哈密頓算法(EHA). 這個(gè)算法屬于漸進(jìn)二階收斂的, 比流形上常用的黎曼梯度算法更加高效. 本文基于正定矩陣流形的幾何結(jié)構(gòu), 提出求解方程(1)的EHA算法. 將流形上2個(gè)點(diǎn)之間的測(cè)地距離作為目標(biāo)函數(shù)，并計(jì)算目標(biāo)函數(shù)在流形上的黎曼梯度. 最后通過(guò)數(shù)值模擬比較EHA算法和經(jīng)典的MSIM算法的計(jì)算行為.

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

本節(jié)介紹黎曼流形上的幾何概念和正定矩陣流形的幾何結(jié)構(gòu)[9-10].

1.1 黎曼流形上的幾何概念

給定X∈M,V∈TXM, 在黎曼流形M上連接2點(diǎn)的最短線被稱之為測(cè)地線, 記γX,V:[0,1]→M. 流形上2點(diǎn)之間測(cè)地距離定義為測(cè)地線的長(zhǎng)度l(γX,V). 如果流形上的任意2點(diǎn)之間都有一條測(cè)地線連接, 那么就能夠測(cè)量這個(gè)流形上任意2點(diǎn)之間的距離.

(2)

1.2 正定矩陣流形上的幾何概念

所有n×n正定對(duì)稱矩陣可以構(gòu)成一個(gè)流形P(n)，稱之為正定矩陣流形. 所有n×n對(duì)稱矩陣組成的集合為S(n)，那么P(n)的切空間為S(n). 即, 任何一個(gè)對(duì)稱矩陣的指數(shù)映射是一個(gè)正定矩陣, 反之，任何一個(gè)正定矩陣的對(duì)數(shù)映射都是一個(gè)對(duì)稱矩陣.

定義P(n)上的Frobenius內(nèi)積

式中：矩陣B1,B2∈P(n)； tr為矩陣的跡. 這個(gè)內(nèi)積賦予P(n)一個(gè)平坦的幾何結(jié)構(gòu), 而且測(cè)地線退化為直線, 并具有一定的局限性.

改進(jìn)的Frobenius內(nèi)積在B處定義為

gB(X,Y)=gI(B-1XB-1Y)=tr(B-1XB-1Y),

(3)

式中X,Y∈S(n). 定義了黎曼度量(3)的P(n)成為一個(gè)黎曼流形, 并且使得P(n)是彎曲的. 在B點(diǎn)以X為方向的測(cè)地線為

顯然這條測(cè)地線無(wú)條件的被完全包含在P(n)里.P(n)上2點(diǎn)之間的測(cè)地距離為

2 廣義哈密頓算法

從而方程(1)的精確解X*在目標(biāo)函數(shù)J(X)的最小值點(diǎn)實(shí)現(xiàn)

首先需要計(jì)算目標(biāo)函數(shù)J(X)的黎曼梯度.

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

(4)

式中a:=ln(X-1L(X))L(X)-1.

根據(jù)式(2)中函數(shù)的黎曼梯度的計(jì)算方法, 將目標(biāo)函數(shù)J(s(t))對(duì)參數(shù)t求導(dǎo)得

根據(jù)文獻(xiàn)[8],P(n)上的廣義哈密頓算法為

(5)

式中μ>0為摩擦因數(shù). 將系統(tǒng)(5)離散化可以表示為

(6)

式中η為學(xué)習(xí)效率. 2個(gè)參數(shù)η和μ滿足

(7)

式中λmax為目標(biāo)函數(shù)J(X)的黑森矩陣的最大特征值. 將式(4)代入式(5), 得到計(jì)算方程(1)數(shù)值解的廣義哈密頓算法(EHA)表示為

離散化之后的EHA算法為

(8)

因此計(jì)算方程(1)數(shù)值解的EHA算法可以表示為

①X0為輸入矩陣X的初始值,V0為是V的初始方向，選取ε>0為容忍系數(shù);

② 計(jì)算ρ(L(X)-X)，式中ρ為矩陣的譜半徑;

③ 如果ρ(L(X)-X)<ε，算法結(jié)束,X為方程(1)的近似解;

④ 通過(guò)式(8)更新Xk和Vk，重復(fù)步驟②.

3 例 證

舉例說(shuō)明廣義哈密頓算法的有效性. 將EHA算法的收斂速度與經(jīng)典的MSIM算法進(jìn)行比較. 關(guān)于MSIM算法的細(xì)節(jié)可以參考文獻(xiàn)[7]. 模擬中的誤差設(shè)為ρ(L(X)-X)<10-15.

例1考慮非線性矩陣方程

式中：

在EHA算法和MSIM算法中, 均取Q為初始值X0. 在MSIM算法中將3階單位矩陣I3作為X1. 在EHA算法中，關(guān)于η和μ的選取與算法迭代次數(shù)的關(guān)系通過(guò)圖1表示出來(lái). 可以看出EHA算法最少需要19次迭代即可實(shí)現(xiàn)事先給定的誤差, 這要求0.85<η<0.95,1.05<μ<1.15,并且這2個(gè)參數(shù)需要滿足式(7). 因此η和μ可以取值為圖2中的實(shí)驗(yàn)數(shù)據(jù). 對(duì)于事先給定的誤差, MSIM算法需要47次迭代, 所以EHA算法比MSIM算法具有更快的收斂速度(見(jiàn)圖3).

例2考慮非線性矩陣方程

式中：

Q=I3.

類似于例1的計(jì)算過(guò)程, 初始值取為X0=Q,X1=2Q,EHA算法和MSIM算法分別需要28次和50次迭代達(dá)到事先給定的誤差, 其中EHA算法中的參數(shù)分別取值為η=0.694,μ=1.328. 所以, EHA算法比MSIM算法具有更快的收斂速度.

4 結(jié) 論

基于正定矩陣流形的幾何結(jié)構(gòu), 提出求解一類非線性矩陣方程的廣義哈密頓算法. 這個(gè)算法具有近2階的收斂速度. 因?yàn)楦倪M(jìn)的Frobenius內(nèi)積所誘導(dǎo)的流形是測(cè)地完備的, 所以將流形上2個(gè)點(diǎn)之間的測(cè)地距離作為目標(biāo)函數(shù). 最后通過(guò)2個(gè)數(shù)值例子表明本文所提出的廣義哈密頓算法比經(jīng)典的多步定常迭代方法收斂速度更快.