韓如意， 王川龍

(1. 太原理工大學(xué) 數(shù)學(xué)學(xué)院， 山西 太原 030024； 2. 太原師范學(xué)院 數(shù)學(xué)系， 山西 太原 030619)

0 引 言

隨著互聯(lián)網(wǎng)公司 (如Google, 360, 百度, Netflix) 的快速發(fā)展, 搜索引擎也逐漸在普及, 使得數(shù)據(jù)的收集越來越方便, 數(shù)據(jù)量猛增, 大數(shù)據(jù)概念被越來越廣泛提及, 而這些大量的數(shù)據(jù)蘊含著非常重要的信息. 所以, 如果能充分運用好這些數(shù)據(jù)信息, 將會為自身的發(fā)展, 公司的利益和社會的進(jìn)步帶來巨大的推動作用. 既然大數(shù)據(jù)對社會和人們的實際生活有重大的價值, 那么就會出現(xiàn)因各種原因而導(dǎo)致的整體數(shù)據(jù)不完整的現(xiàn)象, 研究如何將這些缺失的數(shù)據(jù)補充完整也將變得越來越重要. 這就是矩陣填充問題.

矩陣填充問題近幾年發(fā)展迅速, 它主要研究的實際問題是如何在矩陣元素缺失的情況下通過已知的部分元素精確地填充這些缺失元素， 從而將原矩陣填充完整. 著名的Netflix推薦系統(tǒng)[1]就是矩陣填充問題的一個典型應(yīng)用. 在圖像恢復(fù)[2], 機器學(xué)習(xí)[3], 計算機視覺[4], 控制[5]等信息科學(xué)領(lǐng)域矩陣填充發(fā)揮著重要的作用, 并且其在系統(tǒng)識別[6], 頻譜感知[7], 多媒體編碼和通信[8]等方面也有應(yīng)用.

對于矩陣填充問題有如下模型,

minrank(Z)

s.t.PΩ(Z)=PΩ(Z0),

(1)

式中：Z∈Rn1×n2；Z0∈Rn1×n2且Z0是采樣矩陣；Ω?{1,…,n1}×{1,…,n2} 為已知的樣本元素下標(biāo)集；PΩ為集合Ω的正交投影. 該問題的意義在于將缺失的元素填充后使得矩陣的秩盡可能的低. 近年來, 在已知矩陣秩r的前提下, 人們提出許多算法, 如黎曼信賴域法(Riemannian trust-region method, 簡稱 RTRMC[9])， OptSpace 算法[10], ADMiRA算法[11].

(2)

為求解模型(2) , 人們提出一種交替最小算法PowerFactorization, 簡稱PF[12]. 而算法PF中的每次迭代實際上是求解一個標(biāo)準(zhǔn)的最小二乘問題, 這樣不僅帶來高的迭代計算成本, 還降低了計算效率. 為了改善這個缺陷， 人們提出另一種交替最小算法Low-rankMatrixFitting, 簡稱LMaFit[13]. 2015 年,TannerJ等[14]提出了交替最速下降法 (AlternatingSteepestDecent, 簡稱ASD).ASD比PF的迭代計算成本低很多, 精度更高; 比起LMaFit,ASD能更好地填充高秩矩陣, 并通過數(shù)值實驗也進(jìn)一步證明了ASD算法的可行性和有效性.

在文獻(xiàn)[14]中, 用ASD算法對一般低秩矩陣進(jìn)行了填充. 而在實際中, 采樣矩陣往往也會有許多特殊結(jié)構(gòu)的矩陣, 其中對稱矩陣就是特殊結(jié)構(gòu)矩陣中的一種, 并且對稱矩陣在通信工程和電力系統(tǒng)等諸多工程應(yīng)用領(lǐng)域中也發(fā)揮著重要作用. 因此研究對稱矩陣填充問題很有意義.

本文主要研究對稱矩陣的填充問題, 并結(jié)合最優(yōu)化中的非精確線性搜索方法來搜索最優(yōu)步長, 不僅從理論上證明該算法的可行性， 而且通過不同采樣率的矩陣填充數(shù)值試驗展示該算法的可行性和有效性.

為方便, Rn1×n2表示n1×n2實矩陣全體. ‖Z‖*表示矩陣Z的核范數(shù), ‖Z‖F(xiàn)表示矩陣Z的Frobenius范數(shù)，ZT表示矩陣Z的轉(zhuǎn)置. 矩陣X和Y的內(nèi)積用〈X,Y〉=trace(XTY) 表示. 集合Ω?{1,…,n1}×{1,…,n2}表示采樣矩陣已知元素的下標(biāo)集.PΩ表示集合Ω上的正交投影.

本文第1節(jié)詳細(xì)介紹了非精確線性搜索方法和提出的對稱矩陣填充算法; 第2節(jié)分析了算法的收斂性; 第3節(jié)通過取不同的采樣率進(jìn)行數(shù)值試驗， 進(jìn)一步驗證了算法的可行性和有效性; 第4節(jié)總結(jié)全文.

1 算 法

本節(jié)介紹基于非精確線性搜索的交替最速下降算法去填充一個秩為r的對稱矩陣.

引理 1[15]在數(shù)域P上, 任意一個對稱矩陣都合同于一對角矩陣. 于是, 對稱矩陣填充問題可以描述為如下模型

(3)

式中：Z0∈Rn× n且為對稱矩陣,X∈Rn× r,D∈Rr× r且為對角矩陣.

為了提高算法PF中最小二乘問題的計算效率, 本文算法采用了沿梯度下降方向簡單的線性搜索方法. 設(shè)fD(X)表示f(X,D)對X的梯度;fX(D) 表示f(X,D) 對D的梯度. 于是梯度上升方向為

fD(X)=-2(PΩ(Z0)-PΩ(XDXT))XD,

(4)

沿梯度下降方向的最優(yōu)步長也能被準(zhǔn)確地計算.

引理 2 設(shè)tx,td是下降方向-fD(X) 和-fX(D)對應(yīng)的最優(yōu)步長, 則

(5)

而tx不能被精確求出, 于是用最優(yōu)化中的非精確線性搜索方法求解.

g′(td)=〈PΩ(Z0)-PΩ(X(D-tdfX(D))XT),XfX(D)XT〉=〈PΩ(Z0)-PΩ(XDXT)+

tdPΩ(XfX(D)XT),XfX(D)XT〉=〈PΩ(Z0)-PΩ(XDXT),XfX(D)XT〉+

td〈PΩ(XfX(D)XT),XfX(D)XT〉=〈(PΩ(Z0)-PΩ(XDXT))X,XfX(D)〉+

td〈PΩ(XfX(D)XT),XfX(D)XT〉=〈XT(PΩ(Z0)-PΩ(XDXT))X,fX(D)〉+

td〈PΩ(XfX(D)XT),PΩ(XfX(D)XT)〉=〈-fX(D),fX(D)〉+td〈PΩ(XfX(D)XT),

PΩ(XfX(D)XT)〉=-‖

(6)

顯然h(t)是關(guān)于t的四次函數(shù), 我們用最優(yōu)化中的非精確線性搜索方法求解.

在最優(yōu)化中求步長的非精確線性搜索方法很多, 經(jīng)過多次試驗， 最后選擇帶保險措施的簡單準(zhǔn)則來搜索. 在給出該搜索方法之前, 做如下說明:

設(shè)φ(tx)=h(tx), 則

tx2PΩ(fD(X)D

tx(PΩ(XDfD(X)T)+PΩ(fD(X)DXT))-tx2PΩ(fD(X)DfD(X)T),

(7)

于是

(8)

φ′(tx)=〈PΩ(Z0)-PΩ(XDXT)+txPΩ(XDfD(X)T+fD(X)DXT)-

tx2PΩ(fD(X)DfD(X)T),PΩ(XDfD(X)T+fD(X)DXT)-

2txPΩ(fD(X)DfD(X)T)〉=trace((PΩ(Z0)-PΩ(XDXT)+txPΩ(XDfD(X)T+

2txPΩ(fD(X)DfD(X)T))),

(9)

則

φ′(0)=trace((PΩ(Z0)-PΩ(XDXT))TPΩ(XDfD(X)T+fD(X)DXT)).

(10)

這樣非精確線性搜索中帶保險措施的簡單準(zhǔn)則如下:

1) 選取初始數(shù)據(jù).tx=1,ρ=0.1,a=0.5, minstep=10-5,k∶=1.

2) 檢驗準(zhǔn)則. 計算φ(tx),φ(0),φ′(0).

如果φ(tx)≤φ(0) +ρtxφ′(0), 停止， 輸出tx； 否則， 轉(zhuǎn)第 3) 步.

3) 檢驗保險措施.

如果 ‖tx(-fD(X))‖2

4)tx=atx,k∶=k+1, 轉(zhuǎn)第 2) 步.

以上引理2和非精確線性搜索中帶保險措施的簡單準(zhǔn)則給出沿負(fù)梯度方向的最優(yōu)步長, 接下來介紹對稱矩陣填充算法的具體實施過程.

算法 1 基于非精確線性搜索的交替最速下降算法.

1) 給定矩陣大小n×n, 秩r, 采樣密度δ, 采樣個數(shù)p, 樣本矩陣PΩ(Z0), 初始矩陣X0,D0,ε>0,k∶=1.

2) 計算梯度方向和相應(yīng)的步長, 更新X0,D0.

fDk(Xk)=-2(PΩ(Z0)-

PΩ(XkDkXkT))XkDk,

tx用非精確線性搜索方法求解,

Xk+1=Xk-txfDk(Xk)，

Dk+1=Dk-tdfXk+1(Dk).

5)k∶=k+1, 轉(zhuǎn)第 2) 步.

2 收斂性分析

本節(jié)在一定條件下證明了算法1的收斂性. 為方便,dx和dd分別對應(yīng)X,D的方向, 即負(fù)梯度方向;tx,td分別表示沿負(fù)梯度方向搜索到的最優(yōu)步長.

由于問題(3)的特殊結(jié)構(gòu), 結(jié)合算法中求步長的規(guī)則(5), 能直接證明算法的極值點就是問題(3)的穩(wěn)定點. 也就是說, 極值點(X,D)滿足

fD(X)=0,fX(D)=0.

(11)

定理 1 設(shè)(Xi,Di)是算法1生成的序?qū)? 且在每次迭代中都是非奇異的. 設(shè) (Xik,Dik)是(Xi,Di)的子序?qū)η覍τ谝粋€非奇異序?qū)?(X*,D*)都有

則(X*,D*)是穩(wěn)定點, 即滿足式(11).

為了證明定理1, 由梯度的連續(xù)性, 只須證明

接下來, 利用以下引理3和引理4來證明.

引理 3 收斂于(X*,D*) 的定理1中的子序?qū)?Xik,Dik)有界, 且滿足

(12)

f(Xi,Di-1)=f(Xi-1,Di-1)+〈fDi-1(Xi-1),Xi-Xi-1〉=

f(Xi-1,Di-1)+〈fDi-1(Xi-1),txi-1dxi-1〉=

f(Xi-1,Di-1)+txi-1〈fDi-1(Xi-1),dxi-1〉≥

(13)

(14)

取式(13)下界并從i=0到i=l對式(13)～(14)求和得

而對所有的l, f(Xl,Dl)非負(fù), 于是

則

(15)

將dxi和ddi-1代入式(15)得

(16)

即

(17)

特別地, 將i換成ik, 式(17)也是成立的, 即

(18)

因為Xik,Dik有界, 于是

‖PΩ(

(19)

‖PΩ(Xik

‖Xik‖F(xiàn)‖fXik(Dik-1)‖F(xiàn),

(20)

式中：C>0. 結(jié)合式(18)～(20)有

(21)

引理 4 定理1中的子序?qū)?Xik,Dik)滿足

(22)

所以

‖fXik(Dik)-

(23)

3 數(shù)值實驗

本節(jié)通過數(shù)值實驗給出非精確線性交替最速下降算法對低秩的對稱矩陣填充的實驗結(jié)果, 所有的數(shù)值實驗均是在相同的工作環(huán)境下進(jìn)行的.

算法結(jié)果見表 1～3.p/(n×n)表示采樣密度, 其中p是已知元素的個數(shù).Zout表示得到的填充矩陣,Z0表示采樣矩陣. 在本文算法中, 參數(shù)δ∈{0.3,0.4,0.5},ε=10-3. 此外, 用s表示時間單位,iter表示迭代次數(shù).

表 1 取樣為0.3的填充結(jié)果

表 2 取樣為0.4的填充結(jié)果

表 3 取樣為0.5的填充結(jié)果

表 1～表 3 展示了不同采樣密度的低秩對稱矩陣的填充數(shù)值結(jié)果， 結(jié)果表明該算法是可行的. 隨著矩陣規(guī)模的增大, 計算時間將增大. 但是由于實驗中隨機產(chǎn)生的樣本矩陣和采樣密度, 造成矩陣完成的復(fù)雜度有差異. 為了展示同一樣本矩陣在不同采樣密度下的完成結(jié)果, 進(jìn)行了表 4～表 8 的測試.

表 4 矩陣大小為100×100的填充結(jié)果

表 5 矩陣大小為300×300的填充結(jié)果

表 6 矩陣大小為500×500的填充結(jié)果

表 7 矩陣大小為800×800的填充結(jié)果

表 8 矩陣大小為1 000×1 000的填充結(jié)果

從表 4～表 8 可以看到相同環(huán)境下同一矩陣在不同采樣密度下的具體填充實驗結(jié)果. 發(fā)現(xiàn)隨著采樣密度的增加, 填充時間和迭代次數(shù)也隨著增大, 這在理論分析上也是合理的.

4 結(jié) 論

本文基于一般矩陣填充的ASD算法對具有特殊結(jié)構(gòu)的對稱矩陣填充進(jìn)行了研究. 將對稱矩陣進(jìn)行簡單因式分解, 用ASD算法中的求導(dǎo)思想找到最優(yōu)下降方向, 然后求出相對應(yīng)的步長. 對于不能精確求解的步長, 提出用最優(yōu)化中的非精確線性搜索方法來搜索最優(yōu)步長, 進(jìn)而迭代更新形成新矩陣. 理論上給出了算法的收斂性分析, 表 1～8 的數(shù)值實驗結(jié)果也驗證了算法的可行性和有效性. 同時還可以看到, 填充的矩陣規(guī)模相對較小, 希望可以繼續(xù)探討更大規(guī)模矩陣的填充算法, 這將是今后需要研究的問題.

[1]Bennett J, Lanning S. The netflix prize[C]. Proceedings of KDD Cup and Workshop, 2007: 35.

[2]Bertalmio M, Sapiro G, Caselles V, et al. Image inpainting[C]. Proceedings of the 27th Annual Conference on Computer Graphics and Interactive Techniques, 2000: 417-424.

[3]Sch?lkopf B, Platt J, Hofmann T. Multi-task feature learing[C]. Conference on Advances in Neural Information Processing Systems, 2006: 41-48.

[4]Tomasi C, Kanade T. Shape and motion from image streams under orthography: a factorization method[J]. International Journal of Computer Vision, 1992, 9(2): 137-154.

[5]Mesbahi M, Papavassilopoulos G P. On the rank minimization problem over a positive semidefinite linear matrix inequality[J]. IEEE Transactions on Automatic Control, 1997, 42(2): 239-243.

[6]Moon Y S, Park P, Kwon W H, et al. Delay-dependent robust stabilization of uncertain state-delayed systems[J]. Internation Journal of Control, 2004, 74(14): 1447-1455.

[7]Wang Y, Pandharipande A, Polo Y. Distributed compressive wide-band spectrum sensing[C]. Information Theory and Application Workshop (ITA), San Diego, CA, 2009: 2337-2340.

[8]Pudlewski S, Melodia T. On the performance of compressive video streaming for wireless multimedia sensor networks[C]. IEEE International Conference on Communication, Cape Town, 2010: 1-5.

[9]Boumal N, Absil P A. RTRMC: A Riemannian trust-region method for low-rank matrix completion[C]. International Conference on Neural Information Processing Systems, 2011: 406-414.

[10]Keshavan R H, Oh S, Montanari A. Matrix completion from a few entries[J]. IEEE International Symposium on Information Theory, 2009, 56(6): 324-328.

[11]Lee K, Bresler Y. ADMiRA: Atomic decomposition for minimum rank approximation[M]. Wiley: IEEE Press, 2010: 4402-4416.

[12]Haldar J P, Hernando D. Rank-constrained solutions to linear matrix equations using power factorization[J]. IEEE Signal Processing Letters, 2009, 16(7): 584-587.

[13]Wen Z W, Yin W, Zhang Y. Solving a low-rank factorization model for matrix completion by a non-linear successive over-relaxation algorithm[J]. Mathematical Programming Computation, 2012, 4(4): 333-361.

[14]Tanner J, Wei K. Low rank matrix completion by alternating steepest descent methods[J]. Applied and Computational Harmonic Analysis, 2016, 40(2)： 417-429.

[15]王萼芳, 石生明. 高等代數(shù)[M]. 北京: 高等教育出版社, 2003.