我们主要考虑如下形式的无约束非光滑凸优化问题
许多现代统计和信号应用领域的问题 (如:信号噪声处理,压缩感知等) 都可以归结为具有问题 (1.1) 形式的非光滑凸优化问题. 众所周知,一个常用来估计稀疏变量 x 的技术手段是LASSO[16],利用它可以同时选择和估计变量. LASSO 可以被看作 ℓ1 -范数正则化问题,也就是问题 (1.1) 中的函数形式如
Bakin[1] 在 1999 年提出了适用于变量分块的 LASSO,同时也给出了问题的算法. 他当时提出的问题可以看成是问题 (1.1) 中函数形式为
这个形式的问题相当于用在变量为分块形式时的 LASSO. 更进一步, 当取 KJ=I,这个问题就是块 LASSO 手段[19]. 大家可以参考文献 [2, 10]去了解块 LASSO 技术的应用, 也可以阅读文献[11] 去进一步了解在统计方面块 LASSO 的应用.
但是,由于块 LASSO 手段并不能使变量的分块稀疏, 文献 [7]中提出一个可以使变量分块稀疏的手段. 这个技术手段可以称为稀疏块 LASSO,同时也可以看成问题 (1.1) 中函数形式为
本文主要结果如下.
在本文中,我们提出了一类修正的邻近梯度法来求解相同问题. 在保证算法是线性收敛的基础上, 与邻近梯度法相比,主要是以下方面做出了修正.
(1) 我们提出的修正邻近梯度法中的步长 αk 是自适应的, 而邻近梯度法中的步长需要 3 次或更多次线搜索来得到. 从这个意义上来讲, 本文中的修正邻近梯度法的计算代价要小.
(2) 邻近梯度法的线性收敛性需要步长有一个理论上的下界 α_>0 和其他的前提条件,但是方法本身并不能保证步长不会太小. 而我们提出的修正邻近梯度法中的步长 αk 有一个明确的下界. 这也解决了邻近梯度法的步长可能要很小的缺点.
本文内容安排如下.
修正邻近梯度法的算法框架是在第3节给出, 它的线性收敛性是在第4节给出. 第5节也是最后一节是我们的数值实验.
我们首先介绍几个基本的定义.
定义 2.1 设集合 S⊂Rn,如果对任意的 x1,x2∈S 以及任意的 α∈[0,1],有 αx1+(1−α)x2∈S, 那么称 S 是凸集.
定义 2.2 设集合 S⊂Rn 为凸集,α∈[0,1],f 是定义在 S 上的函数. 如果对任意的 x1,x2∈S,有 f(αx1+(1−α)x2)≤αf(x1)+(1−α)f(x2), 称 f 为 S 上的凸函数.
函数 f 在集合 S 上是μ -强凸函数当且仅当函数 f−(1/2)μ‖⋅‖2 在集合S上是凸函数.
定义 2.3 函数 f 称为 coercive 的,如果 lim‖x‖→+∞f(x)=+∞.
定义 2.4 函数h被称为 L -Lipschitz 连续的,如果存在正数 L,使得对任意 x,y 都有下式成立 ‖h(x)−h(y)‖≤L‖x−y‖. 这里也称 L 是函数 h 的 Lipschitz 常数. 同理,如果称函数 h 是梯度 L -Lipschitz 连续的,就表示其梯度满足 ‖∇h(x)−∇h(y)‖≤L‖x−y‖.
接下来介绍一个求解非光滑凸优化问题 (1.1) 的常用方法: 邻近梯度法 (PGM)[3, 4].
先介绍邻近算子的定义. 对于凸函数 φ(x) (可能非光滑),Moreu-Yoshida 邻近算子[14] proxφ:Rn↦Rn 定义如下
使用邻近算子,我们可以把问题问题 (1.1) 的最优性条件写成一个不动点方程
即使被经常采用,邻近梯度法的收敛性分析还是非常少的. 例如,通过 文献[5,定理3.4],我们只知道如果步长 αk 满足
那么由邻近梯度法 (2.3) 得到的迭代序列 {xk}k≥0 收敛到问题 (1.1) 的一个解. 收敛速率是次线性 O(1/k) (参见文献[12]). 除了一些特殊情况外, 它的线性收敛速率还是没有得到证明. Luo[9],Tseng[13], Luo[20] 和 Zhang[21] 已经证明了邻近梯度法 (PGM) 在 f(x) 满足某些条件时的线性收敛速率,并且这不需要 h(Ax) 是强凸函数的前提. 特别地,在一些应用中,不论是 LASSO,块 LASSO,或者稀疏块 LASSO, 测量值的个数一般是远小于未知量的个数,也就是说 m≪n,所以矩阵 A 不可能列满秩,这意味着 h(Ax) 不可能是强凸函数.
考虑求解问题 (1.1),取函数形式为 (1.3) 或者 (1.4)式. 在这一节,我们将给出解此问题的修正邻近梯度法 (MPGM) 的框架以及分析.
对问题(1.1)的MPGM算法
步 1 任意初始点 x0 和一个比较小的正数 ε,取 k=0.
步 2 计算 yk=proxf(xk−AT∇h(Axk)), 如果 ‖xk−yk‖≤ε,停止;
其它,计算 zk=xk−yk+AT(∇h(Axk)−∇h(Ayk)).
取αk=‖xk−yk‖2/‖zk‖2,计算 xk+1=xk−αkzk.
步 3 如果 ‖xk+1−xk‖≤ε,停止; 否则, 取 k=k+1,转步 2.
引理 3.1 考虑问题 (1.1),ˉX 表示该问题的最优解集合. 假设 ˉX≠∅,那么
(1) ˉX 是一个闭凸集.
(2) x→Ax 在集合ˉX内是个常值函数,即存在 ˉs∈dom(h) 使得对任意 x∈ˉX,有
证 (1) 因为问题 (1.1) 的目标函数 F(x) 是个下半连续凸函数,根据文献 [8,性质1.1.6 和 性质 1.2.2],水平集 {x|F(x)≤minF} 是个闭凸集.
(2) 相同的证明大家可以参考文献 [21,引理 2] 或者文献 [20,性质 2],这里我们就不再赘述.
引理 3.2 假设问题 (1.1) 的最优解集合 ˉX 非空. 对任意 x∈dom(f+h∘A), 令 y(x)=proxf(x−AT∇h(Ax)), ˉX1={x|x=y(x)} 和 ˉX2={x|x=y(x)−AT(∇h(Ax)−∇h(Ay(x)))}. 那么 ˉX 是一个闭凸集,并且 ˉX1=ˉX2=ˉX.
证 对任意 x∈ˉX1,y(x)=x,∇h(Ax)=∇h(Ay(x)),我们有 x∈ˉX2. 因此,ˉX1⊆ˉX2. 假设存在 x∈ˉX2 但是 x 不属于 ˉX1,即‖x−y(x)‖≠0. 那么 ‖x−y(x)+AT(∇h(Ax)−∇h(Ay(x)))‖2=‖x−y(x)‖2+2⟨x−y(x),AT(∇h(Ax)−∇h(Ay(x)))⟩+‖AT(∇h(Ax)−∇h(Ay(x)))‖2≥‖x−y(x)‖2+2⟨Ax−Ay(x),∇h(Ax)−∇h(Ay(x))⟩≥‖x−y(x)‖2+2μ‖Ax−Ay(x)‖2>0. 这与 x∈ˉX2 矛盾. 也就是说,ˉX2⊆ˉX1. 从而, ˉX1=ˉX2.
对任意 x∈ˉX,由问题 (1.1) 的最优性条件和性质, 当且仅当 0∈∂f(x)+AT∇h(Ax)⇔x−AT∇h(Ax)∈x+∂f(x)⇔x=proxf(x−AT∇h(Ax))⇔x∈ˉX1. 因此,ˉX1=ˉX2=ˉX.
注 3.1 引理 3.2 说明如果 MPGM 算法产生的序列收敛,那么它肯定收敛于问题 (1.1) 的某个最优点,这与邻近梯度法的结果一致.
引理 3.3 MPGM 算法中的步长 αk 满足
证 1αk=‖xk−yk+AT(∇h(Axk)−∇h(Ayk))‖2‖xk−yk‖2≤1+2⟨xk−yk,AT(∇h(Axk)−∇h(Ayk))⟩‖xk−yk‖2+‖AT(∇h(Axk)−∇h(Ayk))‖2‖xk−yk‖2≤1+2L‖A‖2+(L‖A‖2)2=(1+L‖A‖2)2. 证毕.
注 3.2 引理 3.3 说明 MPGM 算法的步长有一个明确的正数下界, 因此不可能无限地小. 而邻近梯度法并不能满足这一点.
另外,在每次迭代中,MPGM 算法需要另外计算一次函数梯度, 但是邻近梯度法需要3次或更多次线搜索. 所以 PGM 算法的计算代价要高于 MPGM 算法.
在这一节,我们将证明在上一节提出的修正邻近梯度法 (MPGM) 的线性收敛性. 在此之前,我们回顾并写出我们解决问题的前提条件.
前提条件
(A1) f,h 是下半连续凸函数; f+h∘A 是 coercive 的;
(A2) h 是 μ -强凸函数并且是梯度 L -Lipschitz 连续的,其中 μ,L>0;
(A3) H(s):=∇2h(s) 是 L′-Lipschitz 连续函数,L′>0.
(A4) {误差界条件:} 令 ˉX 表示问题 (1.1) 的最优解集, distˉX(x):=minx∗∈ˉX‖x−x∗‖.
那么对任意的 ξ≥minF(x),存在 κ,ε>0 使得
对任意满足 ‖proxf(x−AT∇h(Ax))−x‖≤ε 和 F(x)≤ξ 的 x 成立.
我们给出下面的引理来推导出 MPGM 算法的线性收敛性证明.
引理 4.1 对任意 x∗∈ˉX,由引理 3.1,ˉs=Ax∗, 存在邻域 B(ˉs,δ),使得当 Ax∈B(ˉs,δ), y=proxf(x−AT∇h(Ax)),η∈[Ax,ˉs],ξ∈[Ax,Ay],有下式成立 ⟨A(x−y),(H(η)−H(ξ))A(x−x∗)⟩−μ‖A(x−x∗)‖2≤0.
证 注意到 ˉX 是闭凸集. 如果 A(x−x∗)=0,那么 x∈ˉX, 并且由引理 3.2,y∈ˉX. 假设 A(x−x∗)≠0. 当 x 充分靠近 x∗ 且 Ax≠Ax∗ 时, 由于临近算子是非扩张的,再由文献 [20,引理 2]有 ‖x−ˉx‖≤κ′‖Ax−ˉs‖, 这里 ˉx:=argminw∈ˉX{‖x−w‖}, 我们有 ‖Ay−Ax∗‖=‖Ay−Aˉx‖≤‖A‖‖y−ˉx‖≤‖A‖‖x−AT∇h(Ax)−ˉx+AT∇h(Aˉx)‖≤‖A‖(‖x−ˉx‖+‖AT(∇h(Ax)−∇h(Aˉx))‖)≤‖A‖(κ′‖Ax−ˉs‖+‖AT(∇h(Ax)−∇h(ˉs))‖)≤(κ′‖A‖+L‖A‖2)‖Ax−ˉs‖=M‖Ax−ˉs‖=M‖A(x−x∗)‖, 其中 M:=κ′‖A‖+L‖A‖2,再由条件 (A3) 和 η∈[Ax,ˉs],ξ∈[Ax,Ay], ⟨A(x−y),(H(η)−H(ξ))A(x−x∗)⟩≤‖A(x−y)‖L′‖η−ξ‖‖A(x−x∗)‖≤‖A(x−y)‖L′{‖A(x−x∗)‖+‖A(x−y)‖}‖A(x−x∗)‖=L′‖A(x−y)‖2‖A(x−x∗)‖+L′‖A(x−y)‖‖A(x−x∗)‖2=L′‖A(x−y)‖{‖A(x−y)‖‖A(x−x∗)‖+‖A(x−x∗)‖2}≤L′‖A(x−y)‖{(‖A(x−x∗)‖+‖A(x∗−y)‖)‖A(x−x∗)‖+‖A(x−x∗)‖2}≤L′‖A(x−y)‖{‖A(x∗−y)‖‖A(x−x∗)‖+2‖A(x−x∗)‖2}≤L′‖A(x−y)‖{M‖A(x−x∗)‖2+2‖A(x−x∗)‖2}≤L′(M+2)‖A(x−y)‖‖A(x−x∗)‖2. 当 x 充分靠近 x∗ 时,‖A(x−y)‖ 无限接近 0, 所以我们有 L′(M+2)‖A(x−y)‖≤μ. 结论得以证明.
定理 4.1 考虑求解问题 (1.1),取函数形式为 (1.3) 或者 (1.4). 如果 ˉX≠∅,那么
对任意 ζ∈R,{x|F(x)≤ζ} 是有界的. 并且误差界条件 (4.1) 成立.
证 证明见文献 [21,定理 2] 和文献 [20,定理1].
现在,我们可以证明我们提出的修正邻近梯度法的线性收敛性了.
定理 4.2 假设函数 f 和 h 满足条件 (A1)--(A4),并且问题 (1.1) 的最优解集 ˉX 非空. 那么由上面的 MPGM 算法得到的迭代序列是收敛的,并且当初始点 x0 离最优解集 ˉX 充分近时,该迭代序列以线性的速率收敛到最优解集 ˉX 中的一个点.
证 在 MPGM 算法中, yk=argminx{f(x)+12‖x−xk+AT∇h(Axk)‖2}, k=0,1,2,⋯. 所以, xk−yk−AT∇h(Axk)∈∂f(yk). 因此对任意的 x∗∈ˉX,由问题 (1.1) 的最优性条件, −AT∇h(Ax∗)∈∂f(x∗). 由 ∂f 的单调性[15],我们有
综上,我们得到序列 {xk} 是有界的并且随着 k→+∞, xk−yk→0. 因此,它存在收敛的子序列 {xk(j)},并且它的极限点 x∞ 是同时满足 {xk(j)}→x∞ 和 yk(j)→x∞ 的. 取任意 x∈Rn 和 w∈∂f(x). 由 (4.2) 式和 ∂f 的单调性 ⟨xk−yk−AT(∇h(Axk))−w, yk−x⟩≥0. 那么沿 k(j) 取极限得到 ⟨−AT(∇h(Ax∞))−w, x∞−x⟩≥0. 因为 ∂f 是极大单调算子[6],我们可以得到 −AT(∇h(Ax∞))∈∂f(x∞). 因此,x∞∈ˉX. 由 (4.5) 式知 ‖xk+1−x∞‖2≤‖xk−x∞‖2−ρ‖xk−yk‖2, 那么,序列 {‖xk−x∞‖} 收敛. 又由于它的子序列 {‖xk(j)−x∞‖} 收敛到 0,所以 {‖xk−x∞‖} 收敛到 0. 从而有,xk→x∞. 记 ˉxk:=argminx∈ˉX‖x−xk‖, k=1,2,⋯, 由于 ˉX 是闭集,从而 ˉxk∈ˉX. 由 (4.5) 式和定理 4.1推得 ‖xk+1−ˉxk+1‖2≤‖xk+1−ˉxk‖2≤‖xk−ˉxk‖2−ρ‖xk−yk‖2≤(1−ρκ2)‖xk−ˉxk‖2, 这意味着 {xk} 以线性速率收敛到最优解集 ˉX 中的一点.
在本节中,我们同时用我们提出的修正邻近梯度法 (MPGM) 和原始的邻近梯度法 (PGM) 来求解取函数形式为 (1.4)式的问题(1.1). 我们选取了几类不同的数据来测试两个算法的效果,并把它们做了比较与分析. 数值实验是在一台具有 Intel Core CPU 2.20 GHz 和2 GB RAM的电脑上进行的. 软件平台是 MATLAB 8.0.0 (R2012b). 终止误差精度 (ε) 取 10−6.
在两个算法中,矩阵 A,向量d 和 初始向量 x0 是随机选取的. 至于函数 (1.4) 中的参数与分块,我们选取 λ=wJ=1,并把集合{1,2,⋯,n} 分为 10 块. 不失一般性,对称正定矩阵 KJ 我们设为单位矩阵. 这是由于对称正定矩阵与单位矩阵合同,选取单位矩阵与其情况等价.
我们测试了两个算法的计算时间和最优目标函数值来显示它们在计算效果方面的差异,数值实验结果见表 1.
在表 1 中,第一列中是矩阵 A 的维数. 邻近梯度法 (PGM) 的测试结果 (计算时间和最优函数值) 分别放在第二和第三列. 对应的,修正邻近梯度法 (MPGM)的结果放在第四和第五列. 表中的数据显示两个算法的最优值是几乎相同的,但是 MPGM 算法是在较少的计算时间内完成的. 这很好地对应了我们在引言的最后部分所 说明的我们所提出的修正邻近梯度法 (MPGM) 相对于邻近梯度法 (PGM) 的改进. 主要原因是我们提出的修正邻近梯度法中的步长 αk 是自适应的,而邻近梯度法中的步长需要3次或更多次线搜索来得到. 同时还有步长 αk 在我们提出的 MPGM 算法中有个明确的正数下界,从而每一步迭代都有明确的行进. 而 PGM 算法可能会出现步长很小的情况,这样可能会导致迭代的效果相对改变不明显.