章建軍*① 曹 杰② 王源源①

?

Stiefel流形上的梯度算法及其在特征提取中的應用

章建軍曹 杰王源源

(南京航空航天大學電子信息工程學院 南京 210016)(南京航空航天大學無人機研究院 南京 210016)

為了提高系統(tǒng)特征提取算法的計算效率、減少占用的存儲空間和簡化程序設計，該文基于Riemann流形上優(yōu)化算法的幾何框架，提出了改進的Stiefel流形上的梯度下降算法。根據不同要求采用不同的測地線計算公式，并使用多項式逼近測地線方程，同時采用了秦九韶-Horner多項式算法及線搜索、變步長的方法。以主分量分析問題為例，詳細討論了Stiefel流形上的梯度算法在其中的應用。理論分析和實驗結果均表明，此方法可以在確保迭代矩陣列向量單位正交性的同時獲得更好的計算效率和收斂速度，并且更容易實現(xiàn)。

Stiefel流形；特征提?。惶荻人惴?；測地線；主分量分析

1 引言

系統(tǒng)特征信息分析方法包括主分量分析(Principal Component Analysis, PCA)、次分量分析(Minor Component Analysis, MCA)和獨立分量分析(Independent Component Analysis, ICA)等，其廣泛應用于數(shù)據壓縮、波達方向估計、信號濾波等信號處理領域。近十年來，眾多專家學者提出了各種算法。然而，現(xiàn)行的各種算法大都是基于歐氏空間的，這些學習算法并沒有充分利用約束集為Stiefel流形這一事實，因而也就不能充分利用Stiefel流形的幾何性質。另外，這類算法往往需要進行特征分解，而特征分解運算量大，算法復雜，工程難以實現(xiàn)。

本文首先詳細討論了Riemann流形上梯度算法，重點討論了Stiefel流形上的梯度算法。其次，對測地線公式作了深入分析與比較，為了能夠在工程實際中應用，作了較大改進，提出了改進的梯度下降算法M-GDM(Modified Gradient Descent Method)。最后，以主分量分析問題為例，詳細討論如何在Stiefel流形上應用梯度下降算法。

2 Riemann流形上的梯度算法

作為求解非線性最優(yōu)化問題的新方法，基于流形的幾何優(yōu)化算法核心思想是把約束集視為高維空間中的曲面，也就是(微分)流形。當在流形上指定了一個Riemann度量(也就是在流形的每一點的切空間上指定一個對稱、正定的2階協(xié)變張量)后，成為Riemann流形。

測地線是Riemann幾何中非常重要的概念，定義如下：

一方面，由于Riemann流形本質上是一種彎曲的空間，所以迭代估計應該沿著流形的測地線進行；另一方面，維Euclidean空間是一種曲率處處恒為0的特殊的Riemann流形，所以無約束最優(yōu)化和約束最優(yōu)化都可以統(tǒng)一成Riemann流形上的最優(yōu)化。然而由于Euclidean空間中的直線就是測地線，所以無約束優(yōu)化的搜索方向可以沿著空間中任何方向，而約束最優(yōu)化的搜索方向則應沿著Riemann流形上測地線所對應的切方向。據此，可以解決梯度算法中的一個關鍵問題搜索方向，并可以得出Riemann流形上無約束優(yōu)化和約束優(yōu)化的梯度類算法的統(tǒng)一框架。

(2) 計算負梯度方向對應的測地線。因為Euclidean空間中測地線就是直線，所以對于無約束最優(yōu)化，就是。然而對于約束最優(yōu)化，則應先計算，然后把投影到流形在點處的切空間上。

在一定條件下，可以證明上面的梯度算法是收斂的。證明過程較長，可以參閱文獻[8]。

3 Stiefel流形及其測地線方程

值得指出的是，Stiefel流形是緊致的，即參數(shù)的取值范圍是一個有限閉集。關于這一點，可以簡單的說明一下，設有，則有

(3)

Riemann流形上梯度下降算法的另一個核心問題測地線的計算。在微分流形與Riemann幾何中測地線方程的推導往往采用內蘊的方式，測地線對應一組微分方程組的解，直接求解并應用到實際中有很大困難。然而Stiefel流形上的測地線卻有相對簡單的表達式，文獻[10]給出了Stiefel流形上測地線的一種計算公式，以定理表述如下：

(5)

(7)

文獻[11]同時也給出了相應的測地線逼近公式，但是需要矩陣求逆和矩陣乘方。

測地線計算公式(4)需要QR分解及矩陣指數(shù)運算，測地線計算公式(6)需要矩陣求逆及矩陣乘方運算，這兩種方法不僅計算量大，而且占用較大的存儲空間、編程復雜。為了減小計算量，同時簡化程序設計，必須簡化測地線的計算方法(當然必須盡可能地滿足單位列正交性的約束條件，這是在流形上進行優(yōu)化的根本目的)。

鑒于此，本文對上面兩種算法中的矩陣指數(shù)運算作近似?？紤]矩陣指數(shù)函數(shù)的Taylor展開式：

由于流形上指數(shù)映射往往限于局部，取值往往比較小(比如取0.01)，又由于Stiefel流形為緊致的，所以最終切向量中的元素取值較小，和中各個元素相乘后，再進行矩陣乘法運算，最終的矩陣元素取值就更小，略去這些項，誤差會非常小，對正交性約束影響也自然非常小。

以上定性分析了誤差較小的原因，下面進行定量分析，給出誤差一個較好的上界。設截取項，表示矩陣范數(shù)，并令，則有

(10)

為了進一步減少計算量，對多項式運算采用秦九韶-Horner算法，即：

梯度算法的第3個關鍵問題是迭代步長的選取。迭代步長的選取本質上是一種“線搜索”或“1維搜索”技術，為了使算法具有更好的性能，應區(qū)分批處理和在線處理這兩種不同的情形，并采用不用的策略，而以往大多數(shù)算法都忽略了這一重要問題。

對于批處理情形，此時相關矩陣已經比較好地估計出來(用時間平均取代統(tǒng)計平均)，一定程度上，此時已經是確定性的最優(yōu)化問題，所以采用不精確線搜索技術Armijo準則。歐氏空間中，Armijo準則是指對于給定，令步長因子，其中為滿足下列不等式的最小非負整數(shù)：

在一些信號處理領域，為了能夠實時處理，必須采用在線處理的方式，由于此時相關矩陣的估計僅僅采用當前時刻的采樣值，所以相關矩陣是非常不精確的和隨機的。為了減小穩(wěn)態(tài)誤差、提高性能，應采用變步長的方法。變步長的方法較多，本文采用比較簡單的變步長公式，其它的變步長方法可以參考自適應濾波相關的文獻等。

4 Stiefel流形上的特征提取

定理2設Stiefel流形上1階連續(xù)可微函數(shù)為,為Stiefel流形上的點。當把Stiefel流形嵌入到Euclidean空間后，在Euclidean空間中的梯度為，投影到流形上點處的切空間得，則和有關系：。

據此，由Riemann流形上梯度算法的一般框架，可以得出Stiefel流形上梯度算法的一般步驟：

必須指出的是，測地線計算公式的選擇很重要，當待求的主分量數(shù)目較少時，應選用式(4)，否則應選用式(6)。

下面以PCA問題為例，討論如何在Stiefel流形上應用梯度下降算法。PCA問題可以等價為Stiefel流形上非線性約束最優(yōu)化問題：

(16)

至于MCA和ICA等問題，只要知道代價函數(shù)，很容易依次類推。

5 數(shù)值仿真與結論

記文獻[11,16]提出的算法為E-GDM(Exponential Gradient Descent Method)。

實驗1 比較正交性。求5維隨機向量的第1，第2主分量，分別設為和。理論上，和應相互正交，即。但是，采用多項式近似后，必然有誤差，即，直接取絕對值作為正交性誤差的度量，即。實驗時隨機選取的5 維隨機向量的協(xié)方差矩陣為：

多項式取7項，步長參數(shù)取0.01，迭代100次，結果如圖1。

從實驗可以看出，用多項式近似和使用指數(shù)映射(exp)其性能是非常接近的，對正交性約束的影響極小，可以認為正交性得以保持。

實驗2 比較運行時間。由于matalb中QR分解算法是用內嵌匯編實現(xiàn)的，為了便于比較，故使用C語言編程實現(xiàn)兩種算法。多項式取7項，向量維數(shù)從40取到100，每次取10000個采樣數(shù)據，E-GDM算法與M-GDM算法占用的時間比如圖2。

可以看出時間比基本穩(wěn)定在15左右，然而本文的算法簡單、占用存儲空間也小。

從實驗可以看出，由于改進的算法采用了“線搜索”技術，收斂速度更快，不僅迭代步數(shù)變小，而且誤差也更小。

6 結束語

如文中所指出，如果在高維數(shù)據(設為)中僅僅提取少量主分量(設為)，則應選用測地線計算公式(4)；當提取的主分量數(shù)接近或者當時，則應選用第2個測地線計算公式(6)。在獨立分量分析問題中，幾乎總是使用計算公式(6)?？偠灾?，應具體問題具體對待。

在Riemann流形(包括Stiefel流形)上發(fā)展其它類型的優(yōu)化算法，并將其應用到信號處理中，以提高算法的效率，不僅是重要的思路，也是一個值得研究的方向。

圖1 正交性誤差比較

圖2 E-GDM和M-GDM占用的時間比

圖3 收斂速度性能曲線

[1] Kong X Y, Hu C H, and Han C Z. A dual purpose principal and minor subspace gradient flow[J]., 2012, 60(1): 197-201.

[2] Eldar Y C and Oppenheim A V. MMSE whitening and subspace whitening[J]., 2003, 49(7): 1846-1851.

[3] Chang D X, Feng D Z, and Zheng W X. A fast recursive total least squares algorithm for adaptive IIR filtering[J]., 2005, 53(3): 957-965.

[4] Ouyang S and Bao Z. Fast principal component extraction by a weighted information criterion[J]., 2002, 50(8): 1994-2002.

[5] Moller R and Konies A. Coupled principal component analysis[J].2004, 15(1): 214-222.

[6] Rapcsak T. On minimization on Stiefel manifold[J]., 2002, 143(2): 365-376.

[7] Jurgen J. Riemannian Geometry and Geometric Analysis[M]. Berlin Heidelberg, NY, US, Springer, 2005: 13-20.

[8] Ferrira O P and Oliveira P R. Sub gradient algorithm on Riemannian manifolds[J]., 1998, 97(1): 93-104.

[9] Yang Y. Global convergent optimization algorithm on Riemannian manifold: uniform framework for unconstrained and constrained optimization[J]., 2007, 132(2): 245-265.

[10] Edelman A, Arias T, and Smith S. The geometry of algorithm with orthogonality constraints[J]., 1998, 20(2): 302-353.

[11] Nishimor Y and Akaho S. Learning algorithm utilizing quasi- geodesic flows on the Stiefel manifold[J]., 2005, 67: 106-135.

[12] 徐樹方, 錢江. 矩陣計算六講[M]. 北京: 高等教育出版社, 2011: 70-71.

Xu Shu-fang and Qian Jiang. Six Lectures on Matrix Computation[M]. Beijing: Higher Education Press, 2011: 70-71.

[13] 黃建國, 孫連山, 葉中行. 黎曼流形上帶Armijo準則步長準則優(yōu)化算法[J]. 上海交通大學學報, 2002, 36(2): 267-271.

Huang J G, Sun L S, and Ye Z X. Optimization algorithm with Armijo rule on Riemann manifold[J]., 2002, 36(2): 267-271.

[14] 張賢達. 現(xiàn)代信號處理[M]. 北京: 清華大學出版社, 2012: 199-200.

Zhang X D. Modern Signal Processing[M]. Beijing: Tsinghua University Press, 2012: 199-200.

[15] Yang H H. Series updating rule for blind separation derived from scoring[J]., 1999, 47: 2279-2285.

[16] 段玲, 黃建國. 主成分分析的一個黎曼幾何隨機算法[J]. 上海交通大學學報, 2004, 38(1): 71-74.

Duan L and Huang J G. A Riemannian geometry underlying stochastic algorithm for adaptive principal component analysis[J]., 2004, 38(1): 71-74.

Gradient Algorithm on Stiefel Manifold and Application in Feature Extraction

Zhang Jian-junCao JieWang Yuan-yuan

(College of Electronic and Information Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing210016, China)(Unmanned Aerial Vehicle Research Department, Nanjing University of Aeronautics and Astronautics, Nanjing210016, China)

To improve the computational efficiency of system feature extraction, reduce the occupied memory space, and simplify the program design, a modified gradient descent method on Stiefel manifold is proposed based on the optimization algorithm of geometry frame on the Riemann manifold. Different geodesic calculation formulas are used for different scenarios. A polynomial is also used to lie close to the geodesic equations. JiuZhaoQin-Horner polynomial algorithm and the strategies of line-searching technique and change of the step size of iteration are also adopted. The gradient descent algorithm on Stiefel manifold applied in Principal Component Analysis (PCA) is discussed in detail as an example of system feature extraction. Theoretical analysis and simulation experiments show that the new method can achieve superior performance in both the convergence rate and calculation efficiency while ensuring the unitary column orthogonality. In addition, it is easier to implement by software or hardware.

Stiefel manifold; Feature extraction; Gradient algorithm; Geodesic; Principal component analysis

TP391

A

2095-283X(2013)03-0309-05

10.3724/SP.J.1300.2013.13048

2013-05-22收到，2013-08-30改回；2013-09-03網絡優(yōu)先出版

國家自然科學基金(61106018)資助課題

章建軍 hblyup@163.com

章建軍(1988-)，男，江蘇南京，碩士研究生，主要研究方向為優(yōu)化算法與盲信號處理。

E-mail: hblyup@163.com

曹 杰(1963-)，男，研究員，研究方向為信號處理、數(shù)字圖像處理。

王源源(1988-)，女，江蘇鎮(zhèn)江，碩士研究生，主要研究方向為數(shù)據壓縮。