基于Lp-Box ADMM的k稀疏编码方法

张雅丽，樊明宇

(温州大学数理与电子信息工程学院，浙江温州 325035)

在如今的数字化时代，随着网络技术的巨大进步和发展，我们每天都面对着来自不同领域的海量数据，在机器学习、信号和图像处理、计算机视觉、模式识别、生物信息等领域中，高维数据也是随处可见的，例如，现代数码相机采集到的图像通常会包括几千万个像素，视频可能由几百万帧图像构成，而文本和在线文档可能与成百上千的特征相关．由于数据特征维度很高，构成复杂，所以这些数据的存储与关系分析就变得更为困难，要求也更高．近些年来，基于稀疏编码的技术成为上述问题的有效解决方法之一，在计算机视觉、语音信号、图像压缩等领域引起了高度关注，得到了广泛且成功的应用．

所谓稀疏编码就是在大量的数据集中，选取很小部分作为元素来重建新的数据，而在信号处理过程中的图像去噪以及在视频处理中更复杂的图像分类和聚类中常常需要处理高维的特征空间，因此采用能够降维的稀疏编码表示就成为行之有效的想法．文献[6]提出了图形正则化非负矩阵分解（Graph Regularized Nonnegative Matrix Factorization，GNMF）方法，在此基础上，文献[7]提出了图形正则化稀疏编码（Graph Regularized Sparse Coding，Graph-SC），这种方法能够考虑到数据的几何结构．与之前的非负矩阵分解和稀疏编码相比，GNMF和Graph-SC在图像分类和聚类上显示出更大的进步．近年来基于稀疏表达的子空间聚类[1-4]获得了极大关注，这些方法首先寻找数据的稀疏表示，然后通过稀疏系数矩阵建立一个相似图从而将数据分割，这些方法并不需要知道子空间的个数和维度，尤其是稀疏子空间聚类（Sparse Subspace Clustering，SSC）[1-2]，不仅有很好的理论支持[5]，而且在很多数据集上效果显著．一些关于稀疏子空间聚类的文章[1-5]中指出，每一个存在于联合低维子空间中的数据点都可以由其它点的线性或者仿射组合得到，并且有多种组合方式．此方法关键在于找到少量与数据点来自相同子空间的点来重构数据点，它克服了局部谱聚类对邻域选择较敏感的缺点，但是由于在解稀疏优化问题时将l0-范数最小化凸松弛为l1-范数最小化，而所得解（即子空间重构系数向量）不一定有意义，也就是说有可能得到稠密的重构系数，而这违背了最初重构系数稀疏性的假设．

为了避免上述问题，本文提出了k稀疏编码（k-SC），不同于SCC有可能会得到未知稀疏程度的表示系数（甚至有可能是稠密解），k-SC通过改变k的值可以得到任意确定稀疏性的解，此外它的求解并未将l0-范数凸松弛为l1-范数，而是将原来的整数稀疏优化问题转化为一个与原问题完全等价的0-1整数规划问题，并进一步采用lp-Box ADMM的方法来求解该等价问题．

本文安排如下：第2节介绍稀疏编码相关的一些研究工作；第3节介绍本文提出的k稀疏子空间线性编码模型，并且给出了模型的优化求解算法；第4节给出新算法在两个重要图像数据集上的实验结果，并与该领域比较成熟的其它算法相比较．

1 相关工作

在介绍相关工作前首先进行符号说明．本文中矩阵用大写字母表示，小写字母则表示列向量，数据点表示为为数据集，为数据点xi在字典下的表示系数，为系数矩阵，向量的l2-范数定义为，l1-范数定义为，l0-范数定义为xi中非零元素的个数．

给定数据矩阵X，字典其中每个di都是字典D的基向量，那么稀疏编码可以通过优化问题（1）来寻找xi的稀疏解

或者也可表示为问题（2）的形式：

然而问题（1）是NP难并且非凸，虽然可以通过组合搜索来求解，但计算复杂度较高．针对上述问题，研究者提出匹配追踪（MP）[8]、正交匹配追踪（OMP）[9]等替代算法，针对形式（2）则提出迭代硬阈值算法（ITH）[10]，此外实验表明在一定条件下l1-范数最小化等价于l0-范数最小化[11]，因此问题（1）等价于：

对于问题（3）的优化可以采用经典的内点法[12]，但计算复杂度较高，因此后来提出一些其它方法，比如增广拉格朗日法（ALM）[13]、原始增广拉格朗日法（PALM）[14]和快速迭代收缩阈值算法（FISTA）[15]．理论上l1-范数是l0-范数的最优凸近似且比l0-范数优化问题更易于求解，但l1-范数优化问题倾向于选择数目较少的一些非常大的值或者数目较多的小值，因此有可能得到稠密解．

2 基于0-1整数规划方法的k稀疏编码方法

我们知道大部分的稀疏编码主要是用l1范数来凸近似l0范数，从而得到优化问题如SSC中的优化目标为：

但是基于凸近似实现的稀疏编码在很多问题中可能是次优解，人们仍期待去直接求解最稀疏的编码问题．本文提出一种k稀疏编码即k-SC，该方法从字典中选取k且仅k个基向量来重构新数据，具体模型如下：

特别地，我们可以将数据矩阵X除去x那一列的部分作为字典D，稀疏系数α的l0范数为k也就意味着α中有且仅有k个非零元，也就是只用字典中对应的k列来重构数据x．直接求解离散非凸优化问题（5）是比较困难的，因此先将问题（5）转化为以下等价问题：

其中对于优化向量v中的元素为0或1．相较于连续优化问题，优化问题（6）对于变量v是一个0-1整数规划问题，很难直接求解，其难点在于离散的整数约束条件．为了解决该问题，可以引入l2-Box将离散的二进制约束条件等价替代为一系列连续约束条件．

命题1 0-1整数集合等价于两个连续集合的交集[16]，

其中p∈ （ 0, ∞ ），1n为n维向量，每个元素均为一，

这里的Sp可看做一个定义在lp空间内，以为中心为半径的（n-1）维的球（也就是说在向量空间n R中，两点间的距离通过lp范数来计算），于是p取2时问题（7）就可以表示为：

其中，那么上述问题便可以采用ADMM 算法来进行优化，从而求得其局部最优解，具体来说，就是算法依次优化各个变量来实现目标函数的优化．

1）固定其它变量，通过解最小二成回归问题优化α；2）固定其它变量，通过解二次优化问题优化v；3）将v分别投影到Sb和Sp上得到v1和v2；4）更新拉格朗日乘子．

针对优化问题（8），得到其增广拉格朗日函数为：

这里是拉格朗日乘子，ADMM每次只更新一个变量而保持其它变量固定，本文用表示第t次迭代时的变量值，则表示第t次迭代时的拉格朗日乘子．

第一步．优化α固定其它变量，关于α的优化可以看做是最小二乘回归问题：

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

第二步．优化v固定其它变量，从（9）知有交叉项和平方项，那么等价于以下二次优化问题：

上式对v求导并让其为零，并令，则有：

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

对任意的x，将其映射在Sb是一个分段函数：

对于在Sp上的映射，首先给出两个候选点x1,x2：

那么x在Sp上的映射为：

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

其中的μ＞1为给定参数．

重复上述步骤直到满足收敛条件和或者达到最大迭代次数．关于更新变量的具体细节在算法1中进行总结．

算法1．

输入：数据矩阵X，正则化参数λ，初始化

输出：稀疏系数α和v

3 实 验

将k稀疏编码算法（k-SC）运用在两个真实数据集（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.）上，实验结果与Hard-10[14]、LLE[15]、logBarrier[16]、PrimalDoul[16]的进行比较，其中Hard-10为l0-稀疏表示问题求解算法，后3种为l1-稀疏表示问题求解算法．此外选取一定的损失函数来评估k稀疏编码算法的收敛性．

3.1 数据集描述

本文使用的YelaB数据集包含38类共2 414张人脸图像，每类图像均包括光照条件和姿态的变化，为减少计算量和内存的限制，将每张图像压缩为32 × 32的分辨率并转化为1 024维的向量．

Coil20数据集是灰度图片集合，包含20个物体从不同角度的拍摄，每隔5度拍摄一副图像，每个物体72张图像，总共1 440张图像，每张图像为128 × 128的分辨率，同样地，为了方便将其处理为32 × 32大小的图像并转化为1 024维的向量．

3.2 参数选择和实验设计

实验中参数的设置直接影响到算法的效果．在这里将所有算法的稀疏系数k取为20，也就是只从数据集中选取20个元素来对某个样本进行线性表示，对于k稀疏编码算法，我们将惩罚项参数λ设为0.01，迭代步长μ设为1.05，而迭代次数则设为300．

实验中选取每个数据集中的第十个数据为重构样本，依次采用上述几种编码方法求解其对应的稀疏系数．实验结果见图1和图2．

图1图2中图的横坐标均为数据集中样本的索引，纵坐标为重构数据（第十个数据），以该数据集为字典通过稀疏编码得到的稀疏编码系．若某横坐标对应的纵坐标为零，则表示在重构样本时并未选取该数据作为重构原子，因此图像中非零纵坐标越少表示求解的系数解越稀疏．此外由于第十个数据属于数据集中的第一类，而我们希望尽可能用同类数据的线性表示来重构样本，所以在所得实验结果中非零纵坐标是否集中在第一类也是重要的评估标准．

以算法在Coil20数据集上的实现为例，图1中的a - e对应的算法依次为log-barrier、Hard-l0、LLE、PrimalDoul和本文的算法k-SC．对比图a与图b - e，可发现通过log-barrier算法得到的系数中非零项数量远超其它算法并且数据分散于整个坐标，也就是该算法得到的系数既不具备稀疏性也并未集中于第一类，这主要是因为该方法用1范数替代0范数，并允许较大误差存在，因此不能实现真正的稀疏编码；另外几个方法均较好地实现了由本类字典重构测试样本．对于Coil-20数据集来说，测试样本对应的稀疏表示应该由前71个字典数据实现，对于YaleB数据集来说，测试样本应该由前 60个字典数据来稀疏表示，根据稀疏系数图像来看，基本上所有的方法都能够用合理的字典数据来重构表示测试样本，说明本文所提方法与当前著名方法的有效性相当．

图1 各类稀疏编码在数据集Coil20上的实验结果Fig 1 Experimental Results on Coil20 Datasets of All Sparse Coding

图2 各类稀疏编码在数据集YaleB上的实验结果Fig 2 Experimental Results on Yale B Datasets of All Sparse Coding

图1中的图f是k-SC算法选取三个损失函数在Coil-20数据集得到的收敛曲线，cost1、cost2和 cost3分别表示随迭代次数增加和函数值的变化曲线，这里的可由上文中的（12）和（13）式得到．从cost1曲线可以看出，在迭代约15次时，函数值由刚开始的5快速减小并接近0，虽然之后有小幅度的波动，但是在迭代约70次后基本稳定在0附近；相比较cost1曲线，cost2曲线也有类似的变化，在快速减小到一定程度后有较小浮动，但在更少的迭代次数约50次时逐渐稳定在0；在三种曲线中收敛最快的为函数值所代表的cost3曲线，不但未出现波动，而且仅在迭代不到10次就已经收敛．因此无论是从哪方面看本文的算法都能够快速收敛，这说明了本文提出的lp-Box ADMM优化算法在求解k-稀疏问题方面具有有效性．

4 结 论

本文提出一种k稀疏编码方法，这种方法可以通过lp-Box ADMM算法得到具有特定k稀疏性的编码系数，也就是用数据集中k个样本来线性重构样本数据，这样不但提高了稀疏编码的可解释性，而且可以得到任意稀疏性的编码系数，方法的有效性与当前著名方法的相当，自身的收敛性也很好．本文提出的这种全新的算法可以为稀疏编码扩展新的解决思路．此外通过设定k的值可以实现数据的降维，从而使数据分析更加容易，在图像聚类等领域有广泛的应用前景．

参考文献

[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.