凸差算法(DC Algorithm)

一、基本概念

1.1凸差问题

为了解决非凸规划问题，1975年引入，一般形式如下：
$inf\{f(x)=g(x)-h(x),x\in \mathbb R^n\}$inf{f(x)=g(x)−h(x),x∈Rn}

其对偶形式为：
$inf\{f^*(y)=h^*(y)-g^*(y):y \in \mathbb R^n\}$inf{f∗(y)=h∗(y)−g∗(y):y∈Rn}

1.2共轭函数(conjugate function)

这个概念和复数的共轭没有一点关系，事实上从牛津词典上conjugate一词只是表示(动词)的变化形式，数学维基百科里也有四种共轭函数的形式，第四种才是凸分析中的共轭函数。函数$f:X\rightarrow R$f:X→R的共轭函数为
$g^*(y)=\mathop{sup}\limits_{x\in X}\{<x,y>-f(x)\}$g∗(y)=x∈Xsup​{<x,y>−f(x)}
共轭函数的自变量变成了$y$y，$x$x变成了它的一个参数。它有一个很好的性质：一定是凸函数，甚至是多面凸函数（下面将讲到什么是多面凸函数）。因为它是一系列关于$y$y的仿射函数(线性函数)中取的最大值，所以斜率肯定是分段增大的，可以从下面这张图很清晰地看出这一点。另一方面也可以直接对$y$y求导。

这个共轭函数通过对$x$x求导并令其等于0可以求出确切解，想进一步了解的可以阅读参考文献[1] 中的3.3节，或者它的译著参考文献[3] 。

1.3次微分(subdifferntial)和次梯度(subgradient)

函数$\varphi$φ在点$x_0$x0​处的次梯度定义如下：
$\partial \varphi(x_0):=\{y_0 \in \mathbb{R}^n:\varphi(x)\geq \varphi(x_0)+<x-x_0,y_0>,\forall x \in \mathbb{R}^n\}$∂φ(x0​):={y0​∈Rn:φ(x)≥φ(x0​)+<x−x0​,y0​>,∀x∈Rn}
一般来说，次梯度是闭凸集，是导数概念在凸函数上的扩展。次微分中的一个元素$y_0 \in \partial \varphi(x_0)$y0​∈∂φ(x0​)称作函数$\varphi$φ在点$x_0$x0​处的次梯度。当且仅当函数$\varphi$φ在$x_0$x0​处可微，那么闭集$\partial \varphi(x_0)$∂φ(x0​)变为导数$\{\nabla \varphi(x_0)\}${∇φ(x0​)}。

比如，绝对值函数$\varphi(x)=|x|$φ(x)=∣x∣：$x<0$x<0时，$\partial \varphi(x)=-1$∂φ(x)=−1;$x>0$x>0时，$\partial \varphi(x)=1$∂φ(x)=1；当$x=0$x=0时，$\partial \varphi=[-1,1]$∂φ=[−1,1]这个闭区间。

1.4$\rho -convex$ρ−convex函数

如果定义在凸集$C$C上的函数$\varphi$φ对于任意$x,y \in C,\lambda\in[0,1]$x,y∈C,λ∈[0,1]，对于某个$\rho \geq0$ρ≥0满足以下不等式
$\varphi(\lambda x+(1-\lambda)y)\leq\lambda \varphi(x)+(1-\lambda)\varphi(y)-\frac{\rho}{2}\lambda(1-\lambda)||x-y||^2$φ(λx+(1−λ)y)≤λφ(x)+(1−λ)φ(y)−2ρ​λ(1−λ)∣∣x−y∣∣2
那么称$\varphi$φ为凸集$C$C上的$\rho -convex$ρ−convex函数。

1.5局部最优条件

原DC规划的充要局部最优条件为：（对偶DC规划问题类似，就不提了）
$\partial h(x^*) \bigcap \partial g(x^*) \neq \empty$∂h(x∗)⋂∂g(x∗)​=∅
这个点就称作临界点，或者说上式是原规划问题的一般KKT条件。并且
$\empty \neq \partial h(x^*)\subset\partial g(x^*)$∅​=∂h(x∗)⊂∂g(x∗)
在很多重要的DC规划问题中也是局部最优的充分条件。

1.6多面凸函数(polyhedral function)

如果一个函数是有限个仿射函数集合中的最大值组成的，那么这个函数一定是多面凸的。也就是说，共轭函数就是多面凸函数。

多面凸差规划是指凸差问题中的函数$g$g或$h$h至少有一个是多面凸的凸差问题。多面凸规划问题在非凸优化和全局优化上发挥了核心作用，是凸差规划和凸差算法的基础，从理论和算法的角度看在局部最优条件和DCA算法的收敛限度上都有非常让人感兴趣的特性，有如下两点：

• 在$h$h是多面凸函数的多面凸问题中，如果$h$h在临界点$x^*$x∗处可导，那么$x^*$x∗实际上就是原DC问题的局部最小值。实际上多面凸函数一般处处可导，所以找到了临界点一般就等价于找到了原DC问题的局部极小值。
• 函数$f$f在临界点$x^*$x∗处是局部凸的，所以就可以用解决凸优化问题的方法来解决原DC问题。这一点的原因是因为$f=g-h$f=g−h在$h$h可导的地方都是局部凸的。

二、凸差算法(DC Algorithm)

2.1迭代格式

基于凸差规划中的局部最优条件和对偶性，凸差算法(DCA)由分别构建原问题和对偶问题的尝试解的两个序列${x_k}$xk​和$y_k$yk​组成，这两个序列满足$g(x^k)-h(x^k)$g(xk)−h(xk)和$h^*(y^k)-g^*(y^k)$h∗(yk)−g∗(yk)是递减的，并且都收敛至局部最优条件，以及
$x^* \in \partial g^*(y^*),y^* \in \partial h(x^*)$x∗∈∂g∗(y∗),y∗∈∂h(x∗)
这两个序列是分别以$x^{k+1}$xk+1、$y^{k+1}$yk+1是凸规化问题$P_k$Pk​和$D_{k+1}$Dk+1​的解的形式计算的，初始点$x^0 \in dom \partial h$x0∈dom∂h，$y^0 \in \partial h(x_0)$y0∈∂h(x0​)。
$(P_k) \quad inf\{g(x)-[h(x^k)+<x-x^k,y^k>]:x \in \mathbb{R}^n\}$(Pk​)inf{g(x)−[h(xk)+<x−xk,yk>]:x∈Rn}

$(D_{k+1})\quad inf\{h^*(y)-[g^*(y^k)+<y-y^k,x^{k+1}>]:y \in \mathbb{R}^n\}$(Dk+1​)inf{h∗(y)−[g∗(yk)+<y−yk,xk+1>]:y∈Rn}

DCA算法有一个非常简单的理解：在第k次迭代中，我们把初始凸差问题中的第二个部分$h$h替换成它的仿射较小值(affine minorization)$h(x^k)+<x-x^k,y^k>$h(xk)+<x−xk,yk>。这个仿射较小值是通过函数h在$x^k$xk点处的次梯度$y^k$yk来定义的，目的是为了生成第$k$k个原始凸规化$(P_k)$(Pk​)，这个规划问题的解恰恰就是$\partial g^*(y^*)$∂g∗(y∗)。对偶地，$(P_k)$(Pk​)的一个解$x^{k+1}$xk+1进而就被用来生成凸规化问题$(D_{k+1})$(Dk+1​)——通过把对偶规划问题$(D_{dc})$(Ddc​)中的第二项$g^*$g∗换成其仿射最小值$(g^*)^{(k)}(y):=g^*(y^k)+<y-y^k,x^{k+1}>$(g∗)(k)(y):=g∗(yk)+<y−yk,xk+1>。这个对偶规划问题的解集也恰恰就是$h$h函数在点$x^{k+1}$xk+1处的微分$\partial h(x^{k+1})$∂h(xk+1)。这个过程一直重复直至收敛。DCA借助函数$h$h和$g*$g∗的次梯度执行了一次双重线性化，从而生成下一个迭代凸规化问题：
$y^k \in \partial h(x^k);\quad x^{k+1} \in \partial g^*(y^k),\forall k\geq0.\\ starting \quad from\quad x^0\in dom \partial h$yk∈∂h(xk);xk+1∈∂g∗(yk),∀k≥0.startingfromx0∈dom∂h

这就是标准的DCA算法迭代格式，中心思想就是把第二项通过线性化替换成局部极小值，然后根据多面凸规化的特性，原问题（对偶问题）就变成了标准的凸优化问题。

2.2收敛特性

DCA是一个不带线性搜索的下降算法，但是具有全局收敛性，享有下面的特点：

(C和D是两个分别包含序列$\{x^k\}${xk}$\{y^k\}${yk}在$\{\mathbb R^n\}${Rn}的凸集)

1. 序列$\{g(x^k)-h(x^k)\}${g(xk)−h(xk)}和$\{h^*(y^k)-g^*(y^k)\}${h∗(yk)−g∗(yk)}是递减的并且

   • $g(x^{k+1})-h(x^{k+1})=g(x^k)-h(x^k)$g(xk+1)−h(xk+1)=g(xk)−h(xk)当且仅当$y^k \in \partial g(x^k)\bigcap \partial h(x^k)$yk∈∂g(xk)⋂∂h(xk),$y^k \in \partial g(x^{k+1})\bigcap \partial h(x^{k+1})$yk∈∂g(xk+1)⋂∂h(xk+1)，并且$[\rho(g,C)+\rho(h,C)]||x^{k+1}-x^k||=0$[ρ(g,C)+ρ(h,C)]∣∣xk+1−xk∣∣=0.而且如果$g$g和$h$h在$C$C上是严格凸的，那么$x^k=x^{k+1}$xk=xk+1.

     此时DCA算法就在第$k$k次迭代终止了(DCA算法的有限收敛)

   • 对偶地有，$h^*(y^{k+1})-g^*(y^{k+1})=h^*(y^k)-g^*(y^k)$h∗(yk+1)−g∗(yk+1)=h∗(yk)−g∗(yk)当且仅当$x^{k+1} \in \partial g^*(x^k)\bigcap \partial h^*(y^k)$xk+1∈∂g∗(xk)⋂∂h∗(yk),$x^{k+1} \in \partial g^*(y^{k+1})\bigcap \partial h^*(y^{k+1})$xk+1∈∂g∗(yk+1)⋂∂h∗(yk+1)，并且$[\rho(g^*,D)+\rho(h^*,D)]||y^{k+1}-y^k||=0$[ρ(g∗,D)+ρ(h∗,D)]∣∣yk+1−yk∣∣=0.而且如果$g^*$g∗和$h^*$h∗在$D$D上是严格凸的，那么$y^{k+1}=y^k$yk+1=yk。

     此时DCA算法就在第$k$k次迭代终止了(DCA算法的有限收敛)

2. 如果$\rho(g,C)+\rho(h,C)>0$ρ(g,C)+ρ(h,C)>0(对应地，$\rho(g^*,D)+\rho(h^*,D)>0$ρ(g∗,D)+ρ(h∗,D)>0)),那么序列${||x^{k+1}-x^k||^2}$∣∣xk+1−xk∣∣2（对应地,$||y^{k+1}-y^k||$∣∣yk+1−yk∣∣)收敛。

3. 如果原DC问题的最优值是有限个的并且无穷序列$x^k$xk和$y^k$yk是有界的，那么每个无穷序列$x^k$xk和$y^k$yk的极限点$\tilde x$x~(对应地,$\tilde y$y~​)都是$g-h$g−h的临界点(对应地,$h^*-g^*$h∗−g∗)。

4. 对于一般地凸差规划，DCA线性收敛。

5. 对于多面凸差问题，DCA有限收敛。

2.3相关问题

凸差分解，DCA中凸规化的求解，DCA初始点的选取，DCA多启动，凸差问题的最近分解技术等都是与凸差规划相关的问题，想进一步研究的可以去阅读参考文献[2].

2.4算法流程

算法：DCA

输入:
T:最大迭代次数
$\epsilon$ϵ:收敛误差
算法流程：
1.令t=0,初始化$x_t$xt​
2.while (t<T) do
3. 计算$y_t \in \partial h(x_t)$yt​∈∂h(xt​);
4. 求解凸优化问题
$argmin \{ g(x)-[h(x_t)+<x-x_t,y_t>]\}$argmin{g(x)−[h(xt​)+<x−xt​,yt​>]}
得到$x_{t+1}$xt+1​
5. if $|x_{t+1}-x_t|\leq\epsilon$∣xt+1​−xt​∣≤ϵ then
6. 算法收敛，跳出循环
7. end if
8. 令t=t+1
9.end while
10.return $x_t$xt​

有一个小细节值得一提的是，步骤4中的凸优化问题中的一项$h(x_t)$h(xt​)一项是一个常数，实际算的时候可以省去。

三、小案例

考虑以下非凸问题：
$inf\{f(x)=x^4-3x^2-x,x\in \mathbb R\}$inf{f(x)=x4−3x2−x,x∈R}
将其DC分解为两个凸函数：
$f(x)=g(x)-h(x)\\ where\quad g(x)=x^4,h(x)=3x^2+x$f(x)=g(x)−h(x)whereg(x)=x4,h(x)=3x2+x

x0=input('请输入初始点位置\n');

ah0=6*x0+1; %h的导数*

Iter=10;it=1;

x=zeros(Iter,1);

y=zeros(Iter,1);

x(1)=x0;y(1)=ah0;

while it<=Iter

​ %求取g的导数

​ [x(it+1),~]=fminunc(@(a)a^4-(a-x(it))*y(it),0);*

​ y(it+1)=6*x(it+1)+1;

​ it=it+1;

end

disp([x,y])

当选取初始点$x=0$x=0时，序列$x_k$xk​经过9次迭代收敛至1.3008：x=[ 0 0.6300 1.0612 1.2258 1.2783 1.2941 1.2989 1.3003 1.3007 1.3008 1.3008]。在$x^*=1.3008$x∗=1.3008这个临界点有$\partial h(x*)=6x^*+1=\partial g(x^*)=4x^{*3}=8.8048$∂h(x∗)=6x∗+1=∂g(x∗)=4x∗3=8.8048,也就是满足了凸差规划问题的局部最优条件。

参考文献

[1]Convex Analysis[2004],Stephen Boyd,et al

[2]Recent Advances in DC Programming and DCA [2014],Tao Pham Dinh

[3]凸优化[2013]，王书宁等译

[4]走进中神通Fenchel