張雅麗，樊明宇

(溫州大學(xué)數(shù)理與電子信息工程學(xué)院，浙江溫州 325035)

在如今的數(shù)字化時代，隨著網(wǎng)絡(luò)技術(shù)的巨大進步和發(fā)展，我們每天都面對著來自不同領(lǐng)域的海量數(shù)據(jù)，在機器學(xué)習(xí)、信號和圖像處理、計算機視覺、模式識別、生物信息等領(lǐng)域中，高維數(shù)據(jù)也是隨處可見的，例如，現(xiàn)代數(shù)碼相機采集到的圖像通常會包括幾千萬個像素，視頻可能由幾百萬幀圖像構(gòu)成，而文本和在線文檔可能與成百上千的特征相關(guān)．由于數(shù)據(jù)特征維度很高，構(gòu)成復(fù)雜，所以這些數(shù)據(jù)的存儲與關(guān)系分析就變得更為困難，要求也更高．近些年來，基于稀疏編碼的技術(shù)成為上述問題的有效解決方法之一，在計算機視覺、語音信號、圖像壓縮等領(lǐng)域引起了高度關(guān)注，得到了廣泛且成功的應(yīng)用．

所謂稀疏編碼就是在大量的數(shù)據(jù)集中，選取很小部分作為元素來重建新的數(shù)據(jù)，而在信號處理過程中的圖像去噪以及在視頻處理中更復(fù)雜的圖像分類和聚類中常常需要處理高維的特征空間，因此采用能夠降維的稀疏編碼表示就成為行之有效的想法．文獻(xiàn)[6]提出了圖形正則化非負(fù)矩陣分解（Graph Regularized Nonnegative Matrix Factorization，GNMF）方法，在此基礎(chǔ)上，文獻(xiàn)[7]提出了圖形正則化稀疏編碼（Graph Regularized Sparse Coding，Graph-SC），這種方法能夠考慮到數(shù)據(jù)的幾何結(jié)構(gòu)．與之前的非負(fù)矩陣分解和稀疏編碼相比，GNMF和Graph-SC在圖像分類和聚類上顯示出更大的進步．近年來基于稀疏表達(dá)的子空間聚類[1-4]獲得了極大關(guān)注，這些方法首先尋找數(shù)據(jù)的稀疏表示，然后通過稀疏系數(shù)矩陣建立一個相似圖從而將數(shù)據(jù)分割，這些方法并不需要知道子空間的個數(shù)和維度，尤其是稀疏子空間聚類（Sparse Subspace Clustering，SSC）[1-2]，不僅有很好的理論支持[5]，而且在很多數(shù)據(jù)集上效果顯著．一些關(guān)于稀疏子空間聚類的文章[1-5]中指出，每一個存在于聯(lián)合低維子空間中的數(shù)據(jù)點都可以由其它點的線性或者仿射組合得到，并且有多種組合方式．此方法關(guān)鍵在于找到少量與數(shù)據(jù)點來自相同子空間的點來重構(gòu)數(shù)據(jù)點，它克服了局部譜聚類對鄰域選擇較敏感的缺點，但是由于在解稀疏優(yōu)化問題時將l0-范數(shù)最小化凸松弛為l1-范數(shù)最小化，而所得解（即子空間重構(gòu)系數(shù)向量）不一定有意義，也就是說有可能得到稠密的重構(gòu)系數(shù)，而這違背了最初重構(gòu)系數(shù)稀疏性的假設(shè)．

為了避免上述問題，本文提出了k稀疏編碼（k-SC），不同于SCC有可能會得到未知稀疏程度的表示系數(shù)（甚至有可能是稠密解），k-SC通過改變k的值可以得到任意確定稀疏性的解，此外它的求解并未將l0-范數(shù)凸松弛為l1-范數(shù)，而是將原來的整數(shù)稀疏優(yōu)化問題轉(zhuǎn)化為一個與原問題完全等價的0-1整數(shù)規(guī)劃問題，并進一步采用lp-Box ADMM的方法來求解該等價問題．

本文安排如下：第2節(jié)介紹稀疏編碼相關(guān)的一些研究工作；第3節(jié)介紹本文提出的k稀疏子空間線性編碼模型，并且給出了模型的優(yōu)化求解算法；第4節(jié)給出新算法在兩個重要圖像數(shù)據(jù)集上的實驗結(jié)果，并與該領(lǐng)域比較成熟的其它算法相比較．

1 相關(guān)工作

在介紹相關(guān)工作前首先進行符號說明．本文中矩陣用大寫字母表示，小寫字母則表示列向量，數(shù)據(jù)點表示為為數(shù)據(jù)集，為數(shù)據(jù)點xi在字典下的表示系數(shù)，為系數(shù)矩陣，向量的l2-范數(shù)定義為，l1-范數(shù)定義為，l0-范數(shù)定義為xi中非零元素的個數(shù)．

給定數(shù)據(jù)矩陣X，字典其中每個di都是字典D的基向量，那么稀疏編碼可以通過優(yōu)化問題（1）來尋找xi的稀疏解

或者也可表示為問題（2）的形式：

然而問題（1）是NP難并且非凸，雖然可以通過組合搜索來求解，但計算復(fù)雜度較高．針對上述問題，研究者提出匹配追蹤（MP）[8]、正交匹配追蹤（OMP）[9]等替代算法，針對形式（2）則提出迭代硬閾值算法（ITH）[10]，此外實驗表明在一定條件下l1-范數(shù)最小化等價于l0-范數(shù)最小化[11]，因此問題（1）等價于：

對于問題（3）的優(yōu)化可以采用經(jīng)典的內(nèi)點法[12]，但計算復(fù)雜度較高，因此后來提出一些其它方法，比如增廣拉格朗日法（ALM）[13]、原始增廣拉格朗日法（PALM）[14]和快速迭代收縮閾值算法（FISTA）[15]．理論上l1-范數(shù)是l0-范數(shù)的最優(yōu)凸近似且比l0-范數(shù)優(yōu)化問題更易于求解，但l1-范數(shù)優(yōu)化問題傾向于選擇數(shù)目較少的一些非常大的值或者數(shù)目較多的小值，因此有可能得到稠密解．

2 基于0-1整數(shù)規(guī)劃方法的k稀疏編碼方法

我們知道大部分的稀疏編碼主要是用l1范數(shù)來凸近似l0范數(shù)，從而得到優(yōu)化問題如SSC中的優(yōu)化目標(biāo)為：

但是基于凸近似實現(xiàn)的稀疏編碼在很多問題中可能是次優(yōu)解，人們?nèi)云诖ブ苯忧蠼庾钕∈璧木幋a問題．本文提出一種k稀疏編碼即k-SC，該方法從字典中選取k且僅k個基向量來重構(gòu)新數(shù)據(jù)，具體模型如下：

特別地，我們可以將數(shù)據(jù)矩陣X除去x那一列的部分作為字典D，稀疏系數(shù)α的l0范數(shù)為k也就意味著α中有且僅有k個非零元，也就是只用字典中對應(yīng)的k列來重構(gòu)數(shù)據(jù)x．直接求解離散非凸優(yōu)化問題（5）是比較困難的，因此先將問題（5）轉(zhuǎn)化為以下等價問題：

其中對于優(yōu)化向量v中的元素為0或1．相較于連續(xù)優(yōu)化問題，優(yōu)化問題（6）對于變量v是一個0-1整數(shù)規(guī)劃問題，很難直接求解，其難點在于離散的整數(shù)約束條件．為了解決該問題，可以引入l2-Box將離散的二進制約束條件等價替代為一系列連續(xù)約束條件．

命題1 0-1整數(shù)集合等價于兩個連續(xù)集合的交集[16]，

其中p∈ （ 0, ∞ ），1n為n維向量，每個元素均為一，

這里的Sp可看做一個定義在lp空間內(nèi)，以為中心為半徑的（n-1）維的球（也就是說在向量空間n R中，兩點間的距離通過lp范數(shù)來計算），于是p取2時問題（7）就可以表示為：

其中，那么上述問題便可以采用ADMM 算法來進行優(yōu)化，從而求得其局部最優(yōu)解，具體來說，就是算法依次優(yōu)化各個變量來實現(xiàn)目標(biāo)函數(shù)的優(yōu)化．

1）固定其它變量，通過解最小二成回歸問題優(yōu)化α；2）固定其它變量，通過解二次優(yōu)化問題優(yōu)化v；3）將v分別投影到Sb和Sp上得到v1和v2；4）更新拉格朗日乘子．

針對優(yōu)化問題（8），得到其增廣拉格朗日函數(shù)為：

這里是拉格朗日乘子，ADMM每次只更新一個變量而保持其它變量固定，本文用表示第t次迭代時的變量值，則表示第t次迭代時的拉格朗日乘子．

第一步．優(yōu)化α固定其它變量，關(guān)于α的優(yōu)化可以看做是最小二乘回歸問題：

如果令那么第t+1次迭代時有：

第二步．優(yōu)化v固定其它變量，從（9）知有交叉項和平方項，那么等價于以下二次優(yōu)化問題：

上式對v求導(dǎo)并讓其為零，并令，則有：

第三步．將v分別投影到Sb和Sp得到v1和v2，

對任意的x，將其映射在Sb是一個分段函數(shù)：

對于在Sp上的映射，首先給出兩個候選點x1,x2：

那么x在Sp上的映射為：

第四步．更新拉格朗日乘子和步長ρ如下：

其中的μ＞1為給定參數(shù)．

重復(fù)上述步驟直到滿足收斂條件和或者達(dá)到最大迭代次數(shù)．關(guān)于更新變量的具體細(xì)節(jié)在算法1中進行總結(jié)．

算法1．

輸入：數(shù)據(jù)矩陣X，正則化參數(shù)λ，初始化

輸出：稀疏系數(shù)α和v

3 實 驗

將k稀疏編碼算法（k-SC）運用在兩個真實數(shù)據(jù)集（YaleB[17-18]和Coil20①Nene S A, Nayar S K, Murase H. Columbia object image library (coil-20). Dept Comput Sci, Columbia Univ, New York, NY, USA,Tech Rep,1996, CUCS-005-96.）上，實驗結(jié)果與Hard-10[14]、LLE[15]、logBarrier[16]、PrimalDoul[16]的進行比較，其中Hard-10為l0-稀疏表示問題求解算法，后3種為l1-稀疏表示問題求解算法．此外選取一定的損失函數(shù)來評估k稀疏編碼算法的收斂性．

3.1 數(shù)據(jù)集描述

本文使用的YelaB數(shù)據(jù)集包含38類共2 414張人臉圖像，每類圖像均包括光照條件和姿態(tài)的變化，為減少計算量和內(nèi)存的限制，將每張圖像壓縮為32 × 32的分辨率并轉(zhuǎn)化為1 024維的向量．

Coil20數(shù)據(jù)集是灰度圖片集合，包含20個物體從不同角度的拍攝，每隔5度拍攝一副圖像，每個物體72張圖像，總共1 440張圖像，每張圖像為128 × 128的分辨率，同樣地，為了方便將其處理為32 × 32大小的圖像并轉(zhuǎn)化為1 024維的向量．

3.2 參數(shù)選擇和實驗設(shè)計

實驗中參數(shù)的設(shè)置直接影響到算法的效果．在這里將所有算法的稀疏系數(shù)k取為20，也就是只從數(shù)據(jù)集中選取20個元素來對某個樣本進行線性表示，對于k稀疏編碼算法，我們將懲罰項參數(shù)λ設(shè)為0.01，迭代步長μ設(shè)為1.05，而迭代次數(shù)則設(shè)為300．

實驗中選取每個數(shù)據(jù)集中的第十個數(shù)據(jù)為重構(gòu)樣本，依次采用上述幾種編碼方法求解其對應(yīng)的稀疏系數(shù)．實驗結(jié)果見圖1和圖2．

圖1圖2中圖的橫坐標(biāo)均為數(shù)據(jù)集中樣本的索引，縱坐標(biāo)為重構(gòu)數(shù)據(jù)（第十個數(shù)據(jù)），以該數(shù)據(jù)集為字典通過稀疏編碼得到的稀疏編碼系．若某橫坐標(biāo)對應(yīng)的縱坐標(biāo)為零，則表示在重構(gòu)樣本時并未選取該數(shù)據(jù)作為重構(gòu)原子，因此圖像中非零縱坐標(biāo)越少表示求解的系數(shù)解越稀疏．此外由于第十個數(shù)據(jù)屬于數(shù)據(jù)集中的第一類，而我們希望盡可能用同類數(shù)據(jù)的線性表示來重構(gòu)樣本，所以在所得實驗結(jié)果中非零縱坐標(biāo)是否集中在第一類也是重要的評估標(biāo)準(zhǔn)．

以算法在Coil20數(shù)據(jù)集上的實現(xiàn)為例，圖1中的a - e對應(yīng)的算法依次為log-barrier、Hard-l0、LLE、PrimalDoul和本文的算法k-SC．對比圖a與圖b - e，可發(fā)現(xiàn)通過log-barrier算法得到的系數(shù)中非零項數(shù)量遠(yuǎn)超其它算法并且數(shù)據(jù)分散于整個坐標(biāo)，也就是該算法得到的系數(shù)既不具備稀疏性也并未集中于第一類，這主要是因為該方法用1范數(shù)替代0范數(shù)，并允許較大誤差存在，因此不能實現(xiàn)真正的稀疏編碼；另外幾個方法均較好地實現(xiàn)了由本類字典重構(gòu)測試樣本．對于Coil-20數(shù)據(jù)集來說，測試樣本對應(yīng)的稀疏表示應(yīng)該由前71個字典數(shù)據(jù)實現(xiàn)，對于YaleB數(shù)據(jù)集來說，測試樣本應(yīng)該由前 60個字典數(shù)據(jù)來稀疏表示，根據(jù)稀疏系數(shù)圖像來看，基本上所有的方法都能夠用合理的字典數(shù)據(jù)來重構(gòu)表示測試樣本，說明本文所提方法與當(dāng)前著名方法的有效性相當(dāng)．

圖1 各類稀疏編碼在數(shù)據(jù)集Coil20上的實驗結(jié)果Fig 1 Experimental Results on Coil20 Datasets of All Sparse Coding

圖2 各類稀疏編碼在數(shù)據(jù)集YaleB上的實驗結(jié)果Fig 2 Experimental Results on Yale B Datasets of All Sparse Coding

圖1中的圖f是k-SC算法選取三個損失函數(shù)在Coil-20數(shù)據(jù)集得到的收斂曲線，cost1、cost2和 cost3分別表示隨迭代次數(shù)增加和函數(shù)值的變化曲線，這里的可由上文中的（12）和（13）式得到．從cost1曲線可以看出，在迭代約15次時，函數(shù)值由剛開始的5快速減小并接近0，雖然之后有小幅度的波動，但是在迭代約70次后基本穩(wěn)定在0附近；相比較cost1曲線，cost2曲線也有類似的變化，在快速減小到一定程度后有較小浮動，但在更少的迭代次數(shù)約50次時逐漸穩(wěn)定在0；在三種曲線中收斂最快的為函數(shù)值所代表的cost3曲線，不但未出現(xiàn)波動，而且僅在迭代不到10次就已經(jīng)收斂．因此無論是從哪方面看本文的算法都能夠快速收斂，這說明了本文提出的lp-Box ADMM優(yōu)化算法在求解k-稀疏問題方面具有有效性．

4 結(jié) 論

本文提出一種k稀疏編碼方法，這種方法可以通過lp-Box ADMM算法得到具有特定k稀疏性的編碼系數(shù)，也就是用數(shù)據(jù)集中k個樣本來線性重構(gòu)樣本數(shù)據(jù)，這樣不但提高了稀疏編碼的可解釋性，而且可以得到任意稀疏性的編碼系數(shù)，方法的有效性與當(dāng)前著名方法的相當(dāng)，自身的收斂性也很好．本文提出的這種全新的算法可以為稀疏編碼擴展新的解決思路．此外通過設(shè)定k的值可以實現(xiàn)數(shù)據(jù)的降維，從而使數(shù)據(jù)分析更加容易，在圖像聚類等領(lǐng)域有廣泛的應(yīng)用前景．

參考文獻(xiàn)

[1]Elhamifar E, Vidal R. Sparse subspace clustering [EB/OL]. [2017-08-18]. http://www.vision.jhu.edu/assets/ElhamifarCVPR09.pdf.

[2]Elhamifar E, Vidal R. Sparse subspace clustering: algorithm, theory, and applications [J]. IEEE Trans Pattern Anal Mach Intell, 2013, 35(11): 2765-2781.

[3]Favaro P, Vidal R, Ravichandran A. A closed form solution to robust subspace estimation and clustering [EB/OL].[2017-08-28]. https://www.researchgate.net/publication/221363914_A_Closed_Form_Solution_to_Robust_Subspace_Estimation_and_Clustering.

[4]Qiu Q, Sapiro G. Learning transformations for clustering and classification [J]. J Mach Learn Res, 2015, 16(1):187-225.

[5]Soltanolkotabi M, Candès E J. A geometric analysis of subspace clustering with outliers [J]. Ann Appl Stat, 2012, DOI:10.1214/12-AOS1034.

[6]Cai D, He X, Han J, et al. Graph regularized non-negative matrix factorization for data representation [J]. IEEE Trans Pattern Anal Mach Intell, 2011, 33(8): 1548-1560.

[7]Zheng M, Bu J, Chen C, et al. Graph regularized sparse coding for image representation [J]. IEEE Trans Image Process, 2011, 21(5): 1327-1336.

[8]Mallat S G, Zhang Z. Matching pursuits with time-frequency dictionaries [J]. IEEE Trans Signal Process, 1993,41(12): 3397-3415.

[9]Pati Y C, Rezaiifar R, Krishnaprasad P S, et al. Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition [EB/OL]. [2017-08-25]. http://www.doc88.com/p-2753097534946.html.

[10]Blumensath T, Davies M. Iterative hard thresholding for compressed sensing [J]. Appl Comput Harmon A, 2009,27(3): 265-274.

[11]Donoho D L. For most large underdetermined systems of linear equations the minimal 1-norm solution is also the sparsest solution [J]. Commun Pur Appl Math, 2006, 59(6): 797-829.

[12]Monteiro R D, Adler I. Interior path following primal-dual algorithms: part I, linear programming [J]. Math Program,1989, 44(1): 27-41.

[13]Yang J, Zhang Y. Alternating direction algorithms forl1-Problems in Compressive Sensing [J]. Siam J Sci Comput,2011, 33(1): 250-278.

[14]Yang A Y, Sastry S, Ganesh A, et al. Fastl1-minimization algorithms and an application in robust face recognition: a review [EB/OL]. [2017-07-18]. http://digitalassets.lib. berkeley. edu/techreports/ucb/text/EECS-2010-13.pdf.

[15]Beck A, Teboulle M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems [J]. Siam J Imaging Sci, 2009, 2(1): 183-202.

[16]Wu B, Ghanem B. lp-Box ADMM: a versatile framework for integer programming [EB/OL]. [2017-08-30].https://arxiv.org/pdf/1604.07666.pdf.

[17]Georghiades A S, Belhumeur P N, Kriegman D J. From few to many: illumination cone models for face recognition under variable lighting and pose [J]. IEEE Trans Pattern Anal Mach Intell, 2012, 3(6): 643-660.

[18]Lee K, Ho J, Kriegman D J. Acquiring linear subspaces for face recognition under variable lighting [J]. IEEE Trans Pattern Anal Mach Intell, 2005, 27(5): 684-698.