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

張雪偉,江祝靈,段雪峰

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

張雪偉,江祝靈,段雪峰

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

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

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

記Rn×n為n×n實(shí)矩陣集合,SRn×n為n×n對(duì)稱半正定矩陣集合,rank(B)和tr(B)分別為矩陣B的秩和跡。對(duì)于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)。‖B‖F(xiàn)表示矩陣B的F范數(shù),則

研究對(duì)稱半正定矩陣的正則化低秩逼近問(wèn)題。

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

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

1 主要結(jié)論

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

由引理1,問(wèn)題1的可行集

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

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

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

其中,

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

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

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

又由泰勒定理可知,

比較上述兩式有

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

證畢。

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

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

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

3)確定搜索方向dk,

4)由Armijo線搜索方法確定搜索步長(zhǎng),即找一個(gè)最小的非負(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)過(guò)23次迭代,得到問(wèn)題2的最優(yōu)解

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

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

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

例2 給定矩陣

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

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

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

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

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

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

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

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

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

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

3 結(jié)束語(yǔ)

利用Gramian分解X=YYT,Y∈Rn×k,刻畫了問(wèn)題1的可行集,將搜索空間維數(shù)降低,并把問(wèn)題1轉(zhuǎn)化為無(wú)約束優(yōu)化問(wèn)題2,通過(guò)Armijo線性搜索的非線性共軛梯度方法對(duì)問(wèn)題2求解。在這個(gè)過(guò)程中,關(guān)鍵在于如何計(jì)算目標(biāo)函數(shù)的梯度,用2個(gè)數(shù)值例子說(shuō)明算法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] 李董輝,童小嬌,萬(wàn)中.數(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

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

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

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