正定矩陣流形上的系統(tǒng)控制算法

納文, 孫福鵬, 曾澍楠, 孫華飛

(北京理工大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，北京 100081)

對(duì)于非線性問題如果強(qiáng)行利用線性方法處理，可能帶來很大誤差，達(dá)不到必要的精度. 幾何方法在處理非線性問題時(shí)發(fā)揮重要的作用. 把所研究的非線性問題納入幾何框架，在流形上每一點(diǎn)處的切空間上定義內(nèi)積-黎曼度量，使之成為黎曼流形，就可以給出幾何刻畫，利用黎曼梯度算法求解定義在黎曼流形上的目標(biāo)函數(shù)的最優(yōu)值. 本文利用黎曼幾何研究正定矩陣系統(tǒng)的控制問題. 在文獻(xiàn)[1-8]中，作者們?cè)O(shè)計(jì)了隨控制系統(tǒng)，使得系統(tǒng)的輸出所服從的概率密度函數(shù)與事先指定的目標(biāo)越接近越好. 其中利用了信息幾何的方法，特別是利用了Fisher信息矩陣充當(dāng)黎曼度量，把Kullback-Leibler散度作為所謂的距離函數(shù)，并給出了算法. 對(duì)于上述問題，由于經(jīng)典信息幾何理論框架的制約，上述算法無法給出更加深入的研究. 另一方面，矩陣信息幾何的主要載體是一般線性群的子流形或李子群，黎曼幾何、拓?fù)鋵W(xué)、纖維叢、代數(shù)拓?fù)?、李群與李代數(shù)等深刻的數(shù)學(xué)都可以發(fā)揮作用. 目前，信息幾何用來研究非隨機(jī)的非線性問題，并在圖像處理、信號(hào)處理等方面獲得了成功的應(yīng)用. 作為啟發(fā)，在文獻(xiàn)[9-14]中，作者利用矩陣信息幾何的思想，設(shè)計(jì)了輸出是正定矩陣的控制系統(tǒng)，使得輸出與事先指定的目標(biāo)越接近越好. 其中，充分利用了矩陣信息幾何優(yōu)越性，即把正定矩陣全體看成黎曼流形，并引入了仿射的黎曼度量來定義測(cè)地距離，在該測(cè)地距離意義下給出問題的自然梯度算法，模擬仿真結(jié)果說明了算法的效果. 一個(gè)問題是：能不能給出更好的距離函數(shù)來求解上述問題. 本文利用對(duì)數(shù)歐氏梯度方法給出改進(jìn)的新算法，關(guān)鍵是在正定矩陣流形以及每一點(diǎn)的切空間建立等距同構(gòu)，使得正定矩陣流形成為一個(gè)歐氏空間，進(jìn)而定義對(duì)數(shù)歐氏距離，使得算法更加簡(jiǎn)單.

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

本文將正定矩陣全體看成一般線性群GL(n,R)={A∈Rn×n|det(A)≠0}的子流形. 由于正定矩陣具有良好的性質(zhì)，特別地，其特征值均為正的性質(zhì)使得矩陣分解具有很好的結(jié)果，在圖像處理、信號(hào)處理等領(lǐng)域獲得成功地應(yīng)用. 用SPD(n)表示n×n正定矩陣所構(gòu)成的微分流形，對(duì)于同一個(gè)流形可以定義不同的黎曼度量使之成為不同的黎曼流形. 定義特殊的黎曼度量-歐氏度量

gA(X,Y)=tr(XY)

(1)

式中X,Y∈TASPD(n),A∈SPD(n).

在文獻(xiàn)[10-11]，作者們提出了仿射黎曼度量

gA(X,Y)=tr(A-1XA-1Y)

(2)

經(jīng)過計(jì)算，可以獲得連接SPD(n)上任意兩點(diǎn)的測(cè)地距離

(3)

在文獻(xiàn)[11]中，作者們利用指數(shù)映射在SPD(n)上定義乘法運(yùn)算

AB=exp(logA+logB)

(4)

在SPD(n)上定義黎曼度量

gA(X,Y)=tr(((dlog)AX)T(dlog)AY))

(5)

式中d表示對(duì)數(shù)映射的微分.

經(jīng)過計(jì)算可獲得連接SPD(n)上任意兩點(diǎn)A,B的測(cè)地距離

d2(A,B)=tr(logA-logB)2

(6)

在文獻(xiàn)[12]中，根據(jù)建立纖維叢上的水平子空間與底流形SPD(n)的切空間之間的等距映射，給出了新的黎曼度量以及測(cè)地距離. 基本想法如下：建立GL(n,R)與SPD(n)之間的黎曼淹沒

(7)

可以獲得一個(gè)直和分解

(8)

由于等距保持測(cè)地線和測(cè)地距離，SPD(n)上的局部測(cè)地距離可由GL(n,R)上的測(cè)地距離來表示. 由于GL(n,R)上任意兩點(diǎn)的測(cè)地距離可以由歐氏距離定義為

(9)

則可得到SPD(n)上任意兩點(diǎn)A,B的測(cè)地距離

(10)

2 控制系統(tǒng)及其算法

本文假設(shè)系統(tǒng)輸出僅由控制輸入決定. 設(shè)x=[x1x2…xm]∈Rm為輸入變量，則系統(tǒng)輸出滿足P(x)∈SPD(n). 本文的目標(biāo)是設(shè)計(jì)輸入x，使輸出y接近給定的目標(biāo)Q(正定矩陣). 此外，希望獲得從初始矩陣P(x0)到最終矩陣P(x*)的控制軌跡.

圖1 控制系統(tǒng)Fig.1 Control system

為了描述輸出矩陣P(x)和目標(biāo)矩陣Q之間的距離，將測(cè)地距離作為它們之間的度量. 因此，目標(biāo)函數(shù)得到如下

J(x)=J(P(x))=d2(P(x),Q)

(11)

正定矩陣系統(tǒng)的最佳輸出x*滿足

(12)

設(shè)輸出流形S={P(x)|x=[x1x2…xm]∈Rm}是SPD(n)的子流形,P(x)是控制輸出矩陣. 下面給出算法.

對(duì)于定義在歐氏空間上的目標(biāo)函數(shù)，則可以利用梯度下降法求解其最小值. 但對(duì)于定義在黎曼流形上的目標(biāo)函數(shù)，此時(shí)的梯度方向已經(jīng)不是歐氏空間的梯度下降方向，黎曼梯度才是最速下降方向. 因此在研究黎曼流形上函數(shù)的最優(yōu)值時(shí)要考慮黎曼梯度. 具體表達(dá)如下.

命題1設(shè)f:M→R是一個(gè)光滑函數(shù),g是黎曼流形M上的黎曼度量. 則有下面的求解f的最小值的迭代公式

θt+1=θt-λgradf(θt)

(13)

注1從命題1可以發(fā)現(xiàn)，自然梯度中要含有黎曼度量g，而且出現(xiàn)了求逆運(yùn)算. 但是在矩陣流形上定義適當(dāng)?shù)睦杪攘浚梢院?jiǎn)單地表達(dá)黎曼度量.

注2無論利用歐氏梯度法還是利用黎曼梯度法求解目標(biāo)函數(shù)的最小值時(shí)，都要陷入所求的的最小值不是整體最小值，而是局部最小值的問題，除非目標(biāo)函數(shù)是一個(gè)凸函數(shù).

既然在文獻(xiàn)[9-10]中，作者們已經(jīng)利用仿射黎曼度量給出了問題的求解過程，在本文中，分別利用定義在SPD(n)上的對(duì)數(shù)歐氏度量以及基于纖維叢黎曼度量給出問題的求解過程. 對(duì)于定義在SPD(n)上的黎曼度量，事實(shí)上，由于對(duì)稱矩陣全體所形成的流形與SPD(n)的等距關(guān)系，SPD(n)成為平坦的流形. 利用歐氏梯度來計(jì)算定義其上面的目標(biāo)函數(shù)的最小值.

命題2設(shè)目標(biāo)Q∈S，則有梯度下降算法

(14)

式中：η為迭代步長,J(xt)為式(11)中的距離函數(shù)

J(xt)=d(P(xt),Q)=tr(logP(xt)-logQ)2

(15)

具體計(jì)算上述距離函數(shù)的梯度，可獲得梯度J(xt)的分量

(16)

把式(16)帶入式(14)可以獲得計(jì)算上述目標(biāo)函數(shù)的最小值.

命題3設(shè)目標(biāo)Q∈S，則有梯度下降算法

xt+1=xt-λgradJ(xt)

(17)

式中：λ為迭代步長；gradJ(xt)為函數(shù)J(xt)的自然梯度；J(xt)為式(11)中的距離函數(shù)

J(xt)=d(P(xt),Q)=tr(P(xt))+tr(Q)-

(18)

從文獻(xiàn)[14]可知，對(duì)于定義在在SPD(n)上的光滑函數(shù)f，在這種由纖維叢誘導(dǎo)的度量下其自然梯度可以表示為

(19)

利用式(19)可以得到梯度grad(P(xt)),從而得到算法的具體迭代公式.

總結(jié)上述過程，則有下面的算法.

命題4對(duì)于控制輸入u=[x1x2…xm]有如下的步驟：

④ 增加t的值，再轉(zhuǎn)到步驟②.

3 模擬仿真

在下面的例子中，將模擬命題2給出2種算法，允許的誤差為10-5. 在第1個(gè)例子中，輸入矩陣和目標(biāo)矩陣為給定的3階正定矩陣，在第2個(gè)例子中，輸入矩陣和目標(biāo)矩陣為隨機(jī)生成的10階正定矩陣.

例1假設(shè)目標(biāo)B是輸出子流形的一點(diǎn). 選取三維正定矩陣

(20)

給定輸出矩陣

(21)

根據(jù)命題2，可以得到從A到B距離的軌跡，如圖2. 在圖2中可以清晰看到基于纖維叢的黎曼梯度在時(shí)間和收斂速度上都明顯優(yōu)于其他兩者，對(duì)數(shù)歐式梯度比黎曼梯度收斂速度略慢，但是所用的時(shí)間卻大大減少，這是由于對(duì)數(shù)歐式梯度簡(jiǎn)化了黎

圖2 輸入矩陣到輸出矩陣在3種梯度下的軌跡(3階)Fig.2 The trajectory from the input to the target under three gradients(3rd order)

曼梯度中復(fù)雜的求逆運(yùn)算，使得運(yùn)算效率大大提高. 模擬得到的結(jié)果與預(yù)期相同.

本文將學(xué)習(xí)速度分別設(shè)定為0.05，0.20，0.40，得到的對(duì)數(shù)歐式梯度和基于纖維叢的黎曼梯度的收斂速度如圖3.

圖3 對(duì)數(shù)歐式梯度和基于纖維叢的黎曼梯度的收斂性Fig.3 Convergence of logarithmic euclidean gradient and Riemannian gradient based on fiber bundle

例1將隨機(jī)生成的2個(gè)10階正定矩陣分別作為輸入和輸出. 根據(jù)命題2得到從輸入到輸出距離的軌跡，如圖4. 從圖4中可以看到基于纖維叢的黎曼梯度的收斂速度依然最優(yōu)，對(duì)數(shù)歐式梯度的收斂速度快于黎曼梯度，這說明隨著輸入和輸出復(fù)雜度的改變，對(duì)數(shù)歐式梯度的優(yōu)勢(shì)已經(jīng)開始顯現(xiàn).同時(shí)本例說明命題2給出的算法不依賴于初值和終值的選取，具有很好的普遍性.

圖4 輸入矩陣到輸出矩陣在3種梯度下的軌跡(10階)Fig.4 The trajectory from the input to the target under three gradients(10th order)

4 結(jié) 論

根據(jù)定義在正定矩陣流形上的對(duì)數(shù)歐氏度量以及基于纖維叢的黎曼度量，本文給出了正定矩陣流形系統(tǒng)的控制問題的求解過程. 與仿射黎曼度量相比，這兩個(gè)度量下測(cè)地距離的形式簡(jiǎn)單，求解速度比仿射黎曼度量的方法快. 問題的解決過程涉及深刻的李群與李代數(shù)、纖維叢的巧妙應(yīng)用，為解決類似的問題提供了參考.