L1极小化问题的一种Gauss-Seidal算法

张梦兰， 李董辉

(华南师范大学数学科学学院，广州 510631)



L1极小化问题的一种Gauss-Seidal算法

张梦兰， 李董辉*

(华南师范大学数学科学学院，广州 510631)

采用罚函数法与Gauss-Seidal算法相结合的思想研究求解L1极小化问题的数值算法:把L1正则化问题视为对L1极小化问题的一种罚函数,由于该函数是非光滑函数,采用光滑化函数对其进行光滑逼近；在此基础上,对此无约束光滑极小化问题采用Gauss-Seidal迭代法求其某种形式的非精确解；再通过合理调整罚参数和光滑化参数, 使得算法产生点列收敛于L1极小化问题的解；最后,通过数值试验测试文中算法的效果, 并从数值计算角度与已有算法进行比较, 结果表明,文中算法具有很好的数值效果.

线性方程组稀疏解； L1极小化； 外点罚函数；Gauss-Seidal迭代

求线性方程组的稀疏解可以表述为L0极小化问题

min‖x‖0,

(1)

(2)

其中,υ>0是参数. 当υ→0时, 问题(2)的解逼近于问题(1)的解.

由于‖x‖0不连续,问题(1)、(2)都是NP-难问题[1], 考虑对其进行松弛，以下为常用的一种松弛形式

(3)

或

minF(x,υ)‖Ax-b‖‖x‖1,

(4)

在一定条件下, 问题(3)和问题(1)在很大概率上有着相同的最优解[2].

线性方程组稀疏解问题在很多领域中都有广泛的应用，如信号恢复[2]、图像去噪[3]、信号的稀疏表示[4]和图像压缩[5]等.关于问题(3)和问题(4), 已有一系列数值效果较好的算法[6].然而，已有的研究通常对这2个问题分别考虑. 本文把问题(4)作为问题(3)的罚函数，利用罚函数法结合Gauss-Seidal算法求解问题(3).

Gauss-Seidal算法是求解线性方程组最受欢迎的经典算法之一, 并已应用于求解非光滑优化问题；如文献[7]用Gauss-Seidal算法来求解逻辑分类问题；文献[8]利用Gauss-Seidal算法来加速近邻点法求L1/TV去噪模型.

本文结合求解最优化问题的罚函数法和解线性方程组的Gauss-Seidal算法提出求解问题(4)的1个算法，称之为求解L1极小化问题的Gauss-Seidal算法, 并证明其收敛性.数值结果表明,本文算法有很好的数值效果.

1　算法及其收敛性

首先,简单回顾求解线性方程组的Gauss-Seidal算法.

考察线性方程组

Ax=b,

(5)

求解该方程组的Gauss-Seidal 的迭代公式为:

xl+1=M-1Uxl+M-1b(l=0,1,…)，

(6)

其中, M和U分别是A的下三角部分和严格上三角部分.众所周知，该算法具有如下收敛性定理.

定理1[9]设A是具有正对角元的半正定矩阵,从任意初始点x0出发，{xl}是用Gauss-Seidal算法求解线性方程组(5)产生的点列, 则{xl}收敛于方程组(5)的一个解.

注意到问题(4)是1个严格凸最优化问题, 其稳定点方程与问题(4)等价.根据定理1, 采用类似于解线性方程组的Gaussi-Seidal算法求该问题的稳定点方程.

由于函数‖x‖1是不光滑的, 难以保证其收敛性.处理这种不光滑函数的1种有效方法是对其进行光滑逼近. 即用某个光滑函数近似‖x‖1, 通过求解相应的光滑问题得到原问题的1个近似解.函数‖x‖1的光滑化函数有很多[3,10-12].本文采用文献[12]的函数，该函数基于如下近似关系

(7)

(8)

其中g(x)是‖x‖1的广义梯度函数.

利用‖x‖1的光滑化函数fθ(x), 得到问题(4)的如下光滑化问题

(9)

其中υk是罚参数，θk是光滑化参数.

本文算法的基本框架是：在每次迭代，利用Gauss-Seidal算法非精确求解子问题(9), 使得点xk满足：

(10)

通过调节罚参数{υk}和光滑化参数{θk}, 使得子问题(9)的非精确解序列{xk}趋于问题(3)的解.

下面的定理表明，若参数选取适当，则子问题(9)的非精确解收敛于问题(3)的1个解.

其中g(x*)由式(9)定义.

▽Qk(xk)=AT(Axk-b)+υk▽fθk(xk) ≤tk*υk→0.

因为A行满秩, 所以x*是问题(3)的可行解, 即

Ax*=b.

(11)

注意到

(12)

从式(11)、(12)可知x*是问题(3)的KKT点.证毕.

下面给出非精确求解子问题(9)的Gauss-Seidal迭代法：即求满足式(10)的点xk. 为了方便, 在本节中省略下标k.

子问题(9)的一阶必要条件为

▽Q(x)=AT(Ax-b)+υ▽fθ(x)=0.

(13)

考虑到▽fθ(x)是非线性的，采用如下 Gauss-Seidal 公式得其迭代格式

Mxl+1=Uxl+ATb-υ▽fθ(xl),

(14)

其中ATA=D-L-U 是ATA矩阵的三角分裂, M=D-L.由于ATA对称半正定，且其对角元是正的，因此方向

dl=xl+1-xl=-M-1▽Q(xl)

(15)

是函数Q(x)在xl的下降方向.

为了保证{Q(xl)}单调递减, 设定步长小于1/2, 且υ/θ< dmin,其中dmin是ATA的对角元中的最小值 .

定理3对子问题(9), 假设步长小于1/2, 且υ/θ < dmin,{xl}由式(14)生成. 那么{Q(xl)}是单调减数列.

证明采用反证法. 假设{Q(xl)}不是单调减数列, 则存在某个l, 使得Q(xl)

dlT(ATA+D)dl=2dlTMdl=-2▽Q(xl)Tdl.

另一方面,

-▽Q(xl)Tdl=[▽Q(xl*)-▽Q(xl)]Tdl=

[ATA(xl*-xl)+υ(▽fθk(xl*)-▽fθ(xl))]Tdl=

δ*dlTATAdl+υ(▽fθ(xl*)-▽fθ(xl))Tdl.

所以

dlT(ATA+D)dl=2δ*dlTATAdl+

2υ(▽fθ(xl*)-▽fθ(xl))Tdl.dlTDdl+(1-2δ*)dlTATAdl=υ(▽fθ(xl*)-▽fθ(xl))Tdl.因为δ*<1/2,则

dlTDdl≤2υ(▽fθ(xl*)-▽fθ(xl))Tdl.

所以

即

(16)

式(16)与υ/θ

以下定理证明{xl}收敛于子问题(9)的解.

定理4设{xl}由Gauss-Seidal迭代法产生，其中步长小于1/2, 且υ/θ

证明由定理3知,Q(xl+1)

2　数值实验

对本文算法(记为SPGSA算法)进行数值试验, 检验其数值效果,并与已有的求解L1极小化问题的常用算法(Homotopy、PALM和FISTA算法)进行数值比较，运行环境为PC电脑、Intel Core i3、Quad CPU 2.40GHz和2G RAM memory，采用Matlab(R2012a)编程，Homotopy、PALM和FISTA算法的代码来源于文献[13]提供的网址.

A=randn(m,n)，xg=zeros(n,1)， r=randperm(n)，

x0(r(1∶p))=10*(randn(p,1))， bg=A*xg,

其中n是xg的维数,m是方程的个数,p是xg中非零元的个数. 令xEst为输出的迭代值,bEstAxEst，采用以下误差: xError‖xEst-xg‖∞,以下实验中xError≤1e-2. 2.1无噪声情况下的算法比较

首先比较SPGSA算法与Homotopy、PALM和FISTA算法在无噪声情况下的数值结果.从图1可以看出:(1)Homotopy算法在非零元很少(少于10%)的时候消耗时间很少,但是随着非零元增多,其计算消耗时间越来越长,增长得很快；(2)PALM、SPGSA和FISTA算法在解不定方程组(n>m)问题中差异不大,但在解超定方程组(n

图1　无噪声情况下各算法的耗时比较

2.2有噪声情况下的算法比较

下面测试算法对于有噪声问题的处理能力.问题模型为

其中根据实际问题, ‖·‖可以是L1范数, 也可以是l∞范数、l2范数.这个问题模型也可以用以下问题模型来近似

图3比较了不同算法在解超定方程组(n

蓝点表示噪声, 红点表示bEst-bg.

图2SPGSA算法中bEst和bg的差异(n=4 000,m=5 000)

Figure2ThedifferencebetweenbEstandbginthealgorithms(n=4 000,m=5 000)

图3　有噪声情况下各算法耗时比较(n=4 000,m=5 000)

Figure3ComparisonoftheCPUtimeofthealgorithmswithnoise(n=4 000,m=5 000)

[1]GIBBONSLE,HEARNDW,PARDALOSPM,etal.Continuouscharacterizationsofthemaximalcliqueproblem[J].MathematicsofOperationsResearch, 1997,22(3):754-768.

[2]DONOHODL.Compressedsensing[J].IEEETransactionsonInformationTheory, 2006,52(4):1289-1306.

[3]CANDESEJ,ROMBERGJK,TAOT.Robustuncertaintyprinciples:exactsignalreconstructionfromhighlyincompletefrequencyinformation[J].IEEETransactionsonInformationTheory, 2006,52(2):489-509.

[4]FUCHSJJ.Onsparserepresentationsinarbitraryredundantbases[J].IEEETransactionsonInformationTheory, 2004,50(6):1341-1344.

[5]ELADM,BRUCKSTEINAM.Ageneralizeduncertaintyprincipleandsparserepresentationsinpairsofbases[J].IEEETransactionsonInformationTheory, 2002,48(9): 2558-2567.

[6]文再文, 印卧涛, 刘歆,等. 压缩感知和稀疏优化简介[J].运筹学学报,2012,16(3):49-64.

WENZW,YINWT,LIUX,etal.Introductiontocompressivesensingandsparseoptimization[J].OperationsResearchTransactions,2012,16(3):49-64.

[7]SHEVADESK,KEERTHISS.Asimpleandeffcientalgorithmforgeneselectionusingsparselogisticregression[J].Bioinformatics, 2003,19(17):2246-2253.

[8]LIQ,MICCHELLICA,SHENLX,etal.AproximityalgorithmacceleratedbyGauss-SeideliterationsforL1/TVdenoisingmodels[J].InverseProblems, 2012,28(9):101-120.

[9]SERRED.Matrices:theoryandapplication[M].Berlin:Springer,2002:152-155.

[10]BECKA,TEBOULLEM.Afastiterativeshrinkage-thresholdingalgorithmforlinearinverseproblems[J].SIAMJournalonImagingScience, 2009,2(1):183-202.

[11]DONOHODL,MALEKIA,MONTANARIA.Messagepassingalgorithmsforcompressedsensing[J].ProceedingsoftheNationalAcademyofSciencesoftheUnitedStatesofAmerica, 2009,106(45):18914-18919.

[12]AYBATNS,IYENGARYG.Afirst-ordersmoothedpenaltymethodforcompressedsensing[J].SIAMJournalonOptimization, 2011,21(1):287-313.

[13]YANGAY,SASTRYSS,GANESHA,etal.Fastl1-minimizationalgorithmsandanapplicationinrobustfacerecognition:areview[C]∥IEEEInternationalConferenceonImageProcessing.HongKong:IEEE,2010:1849-1852.

【中文责编：庄晓琼英文责编：肖菁】

Gauss-Seidal Algorithm to L1Minimization

ZHANG Menglan, LI Donghui*

(School of Mathematics Science,South China Normal University,Guangzhou 510631, China)

The numerical method for solving theL1minimization problem is studied.The penalty method and the Gauss-Seidal iteration technique are adopted to develop the method. TheL1regularization problem is regarded as a penalizedL1minimizationproblem. Taking into account that theL1norm function is nonsmooth, a smoothing function is adopted as an approximation to it. Then a smooth unconstrained optimization problem is solved inexactly via Gauss-Seidal iteration. At last,the penalty parameter is adjusted in an appropriate way so that the generated sequence of iterates converges to a solution of theL1minimization problem. The proposed method is tested by numerical experiments, and its performance is compared with some existing methods. The results show that the proposed method is practically effective.

sparse solution of the system of linear equations;L1norm minimization problem; penalty function; Gauss-Seidal iteration

2016-04-08《华南师范大学学报(自然科学版)》网址：http://journal.scnu.edu.cn/n

国家自然科学基金项目(11371154)

李董辉，教授，Email:20091126@m.scnu.edu.cn.

O22;O

A