張雪偉,江祝靈,段雪峰

基于Gramian分解的對稱半正定矩陣的正則化低秩逼近

張雪偉,江祝靈,段雪峰

(桂林電子科技大學(xué)數(shù)學(xué)與計算科學(xué)學(xué)院,廣西桂林 541004)

針對對稱半正定矩陣的正則化低秩逼近問題,基于對稱半正定矩陣的Gramian分解,將對稱半正定矩陣的正則化低秩逼近問題轉(zhuǎn)化為等價的無約束優(yōu)化問題,并構(gòu)造非線性共軛梯度方法求解轉(zhuǎn)化后的無約束優(yōu)化問題。數(shù)值實(shí)驗(yàn)驗(yàn)證了新方法的可行性。

對稱半正定矩陣;正則化低秩逼近;非線性共軛梯度法;Gramian分解

記Rn×n為n×n實(shí)矩陣集合,SRn×n為n×n對稱半正定矩陣集合,rank(B)和tr(B)分別為矩陣B的秩和跡。對于n×n階矩陣B,設(shè)bi為矩陣B的第i列,則B可表示為B=(b1,b2,…,bn)。設(shè)[Bij]為矩陣B的第i行第j列元素,即[Bij]=bij。用vec(B)表示矩陣B的按列拉直向量,即vec(B)=(b1T,b2T,…,bnT)?！珺‖F(xiàn)表示矩陣B的F范數(shù),則

研究對稱半正定矩陣的正則化低秩逼近問題。

問題1 給定矩陣A,B,C∈Rn×n,正整數(shù)k≤n,參數(shù)0＜α＜1,求一個秩小于等于k的對稱半正定矩陣ˉX,使得

研究結(jié)構(gòu)矩陣逼近問題采用的方法有低階近似的截斷奇異值分解法[1-2]、Lanczos雙對角化過程[3]、Monte Carlo算法[4]?；诮Y(jié)構(gòu)化總體最小二乘, Park等[5]提出Hankel低階近似的數(shù)值方法,此方法后被用于Sylvester的低秩近似計算,并且在計算Sylvester逼近過程中給出了一元多項(xiàng)式的近似最大公約數(shù)[6-8]。Higham[9]和段雪峰[10]分別采用分半算法和共軛梯度法求解對稱半正定矩陣的近似逼近; Suliman[11]采用交替投影算法和擬牛頓法求解Toeplitz矩陣逼近問題;唐鳴[12]采用牛頓法求解,將問題轉(zhuǎn)化為無約束優(yōu)化問題,但計算量較大,而且牛頓法對于大多數(shù)問題不是整體收斂,而改用擬牛頓法對Hankel矩陣逼近問題進(jìn)行求解得到了較好的收斂效果。對稱半正定矩陣的正則化低秩逼近問題幾乎無人研究,原因在于難以刻畫該問題的可行集。鑒于此,利用Gramian分解X=YYT,Y∈Rn×k,刻畫該問題的可行集,將原問題轉(zhuǎn)化為一個等價的無約束優(yōu)化問題,應(yīng)用帶Armijo線性搜索的非線性共軛梯度法對等價問題求解。

1 主要結(jié)論

引理1 n×n階對稱半正定矩陣X的秩小于等于k,當(dāng)且僅當(dāng)存在矩陣Y∈Rn×k,使得X=YYT。

由引理1,問題1的可行集

Ω={X|X∈SRn×n+,ran k(X)≤k}

可刻畫成X=YYT,Y∈Rn×k,因此,問題1可等價轉(zhuǎn)化為問題2。

問題2 給定矩陣A,B,C∈Rn×n,參數(shù)α∈(0, 1),求一個n×k階矩陣ˉY,使得

其中,

易知,若ˉY是問題2的解,則X=YY是問題1的解。問題2的搜索空間為nk,當(dāng)k較小時,比問題1的搜索空間n2小很多。

定理1 目標(biāo)函數(shù)f(Y)的梯度為:

證明 由增量矩陣E∈Rn×k,有

又由泰勒定理可知,

比較上述兩式有

設(shè)E=eieTj,其中ei為n×n單位矩陣的第i列元素,ej為k×k單位矩陣的第j列元素,可得

證畢。

利用Armijo線搜索下的非線性共軛梯度方法求解問題2。

算法1(非線性共軛梯度算法) 1)初始化ρ∈(0,1),δ∈(0,0.5),σ∈(δ,0.5)。選取誤差0≤t?1,取初始值迭代矩陣Y0∈Rn×k,置k∶=0。

2)計算gk=▽f(Yk),若‖gk‖F(xiàn)＜t,停機(jī)并輸出≈Yk。

3)確定搜索方向dk,

4)由Armijo線搜索方法確定搜索步長,即找一個最小的非負(fù)整數(shù)mk,使得

置αk=,ζk+1=ζk+αkdk。

5)令k∶=k+1,并轉(zhuǎn)步驟2)。由文獻(xiàn)[13]的定理,可建立算法1的收斂性定理。定理2 假設(shè)f連續(xù)可微且有界,若其梯度▽f是Lipschitz連續(xù)的,那么存在常數(shù)L＞0,使得

‖▽f(X)-▽f(Y)‖F(xiàn)≤L‖X-Y‖F(xiàn),?X,Y∈Rn×k。那么由算法1生成的序列{Yk}滿足

2 數(shù)值實(shí)驗(yàn)

實(shí)驗(yàn)在Matlab 2009a環(huán)境下運(yùn)行,設(shè)▽fk=▽f(Yk),其中Yk是算法1中第k次迭代值,停機(jī)標(biāo)準(zhǔn)為‖▽fk‖＜1×10-4。

例1 給定矩陣

若k=1,α=0.3,初始值

用算法1求解。經(jīng)過23次迭代,得到問題2的最優(yōu)解

從而問題1的最優(yōu)解為

其中rank(X*)=1,X*為對稱半正定矩陣,其目標(biāo)函數(shù)的梯度范數(shù)曲線如圖1所示。

圖1 n=3,k=1時的梯度范數(shù)曲線Fig.1 Gradient norm curve when n=3,k=1

例2 給定矩陣

1)若k=2,α=0.3,初始值

用算法1求解。經(jīng)過280次迭代,得到問題2的最優(yōu)解

從而問題1的最優(yōu)解為

其中rank(X*)=2,X*為對稱半正定矩陣,其目標(biāo)函數(shù)的梯度范數(shù)曲線如圖2所示。

圖2 n=5,k=2時的梯度范數(shù)曲線Fig.2 Gradient norm curve when n=5,k=2

2)若k=1,α=0.3,初始值

用算法1求解。經(jīng)過75次迭代,得到問題2的最優(yōu)解

從而問題1的最優(yōu)解為

其中rank(X*)=1,X*為對稱半正定矩陣,其目標(biāo)函數(shù)的梯度范數(shù)曲線如圖3所示。

圖3 n=5,k=1時的梯度范數(shù)曲線Fig.3 Gradient norm curve when n=5,k=1

3 結(jié)束語

利用Gramian分解X=YYT,Y∈Rn×k,刻畫了問題1的可行集,將搜索空間維數(shù)降低,并把問題1轉(zhuǎn)化為無約束優(yōu)化問題2,通過Armijo線性搜索的非線性共軛梯度方法對問題2求解。在這個過程中,關(guān)鍵在于如何計算目標(biāo)函數(shù)的梯度,用2個數(shù)值例子說明算法1的可行性。

[1] Golub G H,Van Loan C F.Matrix Computations[M]. Baltimore:Johns Hopkins University Press,1996:76-81.

[2] Hansen P C.The truncated SVD as a method for regularization[J].BIT Numerical Mathematics,1987,27: 534-553.

[3] Siman H D,Zha Hongyuan.Low rank matrix approximation using the Lanczos Bidiagonalization process withapplications[J].SIAM Journal on Scientific Computing, 2000,21:2257-2274.

[4] Drineas P,Kannan R,Mahoney M W.Fast Monte Carlo algorithms for matrices II:computing a low rank approximation to a matrix[J].SIAM Journal on Computing,2006,36:158-183.

[5] Park H,Zhang Lei,Rosen J B.Low rank approximation of a Hankel matrix by structured total norm[J].BIT Numerical Mathematics,1999,39:757-779.

[6] Li Bingyu,Yang Zhengfeng,Zhi Lihong.Fast low rank approximation of a Sylvester matrix by structured total least norm[J].Japan Society for Symbolic and Algebraic Computation,2005,11:165-174.

[7] Winkler J,Allan J.Structured low rank approximations of the Sylvester resultant matrix for approximate gcds of bernstein basis polynomials[J].Electronic Transactions on Numerical Analysis,2008,31:141-155.

[8] Winkler J,Allan J.Structured total least norm and approximate gcds of inexact polynomials[J].Journal of Computational and Applied Mathematics,2008,215:1-13.

[9] Higham N J.Computing a nearest symmetric positive semi-definite matrix[J].Linear Algebra and its Applications,1988,103:103-118.

[10] Duan Xuefeng,Li Jiaofeng,Wang Qingwen,et al.Low rank approximation of the symmetric positive semidefinite matrix[J].Journal of Computational and Applied Mathematics,2014,260:236-243.

[11] Al-Homidan S.Toeplitz matrix approximation[J]. Mathematical Sciences Research Journal,2002,2:104-111.

[12] 唐鳴.低秩Hankel矩陣逼近及其加權(quán)逼近的算法[D].南京:南京航空航天大學(xué),2006:7-33.

[13] 李董輝,童小嬌,萬中.數(shù)值最優(yōu)化算法與理論[M].北京:科學(xué)出版社,2005:75-81.

編輯:翁史振

Regularized low rank approximation of symmetric positive semi-definite matrix based on Gramian decomposition

Zhang Xuewei,Jiang Zhuling,Duan Xuefeng
(School of Mathematics and Computational Science,Guilin University of Electronic Technology,Guilin 541004,China)

The regularized low rank approximation of symmetric positive semi-definite matrix is studied.Based on Gramian decomposition of a symmetric positive semi-definite matrix,the regularized low rank approximation of the symmetric positive semi-definite matrix problem is transformed into an equivalent unconstrained optimization problem,and the nonlinear conjugate gradient method is constructed to solve the equivalent unconstrained optimization problem.The numerical experiments verify that the new method is feasible.

symmetric positive semi-definite matrix;regularized low rank approximation;nonlinear conjugate gradient method;Gramian decomposition

O241.6

A

1673-808X(2015)05-0419-05

2015-04-22

國家自然科學(xué)基金(11101100);廣西自然科學(xué)基金(2012GXNSFBA053006)

段雪峰(1982-),男,湖南祁陽人,教授,博士,研究方向?yàn)閼?yīng)用矩陣分析和信息安全。E-mail:zhangxuewwei@163.com

張雪偉,江祝靈,段雪峰.基于Gramian分解的對稱半正定矩陣的正則化低秩逼近[J].桂林電子科技大學(xué)學(xué)報,2015,35(6):419-423.