提示:本文大部分参考Boy-《Distributed_Optimization_and_Statistical_Learning》-2011
目录
前言
交替方向乘子法(Alternating Direction Method of Multipliers,ADMM)是一种求解具有可分离性的凸优化问题的计算框架, 由于其是对偶分解法和增广拉格朗日乘子法的结合,使该算法有分解性的同时保证了良好的收敛性,处理速度快。ADMM适用于求解分布式凸优化问题,主要应用在解空间规模很大的情况,可以进行分块求解,而且解的绝对精度要求不是太高。
ADMM以先分解再结合的形式求解问题,即先把原问题分解成若干个相对原问题较简单的子问题,再把子问题的解结合起来得到原问题的全局解。
注意:本文适合有一点点凸优化理论基础的读者阅读,一点点就ok.
一、预备知识
1. 对偶上升(Dual Ascent)
考虑等式约束凸优化问题
min
f
(
x
)
s.t.
A
x
=
b
(1)
\min \quad f(x)\\ \text{ s.t.} \quad Ax=b \tag{1}
minf(x) s.t.Ax=b(1) 其中变量
x
∈
R
n
x\in R^n
x∈Rn,
A
∈
R
m
×
n
A\in R^{m\times n}
A∈Rm×n, 并且函数
f
:
R
n
→
R
f:R^n \to R
f:Rn→R是凸的。
那么(1)式的拉格朗日函数( Lagrangian)可以写为:
L
(
x
,
y
)
=
f
(
x
)
+
y
T
(
A
x
−
b
)
(2)
L(x,y)=f(x)+y^T(Ax-b) \tag{2}
L(x,y)=f(x)+yT(Ax−b)(2)上式的对偶函数(Dual Function)为
h
(
y
)
=
inf
x
L
(
x
,
y
)
=
−
f
∗
(
−
A
T
y
)
−
b
T
y
(3)
h(y)=\inf_xL(x,y)=-f^*(-A^Ty)-b^Ty \tag{3}
h(y)=xinfL(x,y)=−f∗(−ATy)−bTy(3)其中
y
y
y叫做对偶变量(Dual variable)或者拉格朗日乘子(Lagrange multiplier),
y
∈
R
m
y\in R^m
y∈Rm,
f
∗
f^*
f∗是
f
f
f的共轭函数(Conjugate function)
注
1
^{注1}
注1。
注1:这里简单介绍一下共轭函数以及上式中对偶函数的推导,如想深入理解对偶函数的意义请参考Boyd 《Convex Optimization》中的3.3节。
\quad
函数
f
(
x
)
f(x)
f(x)的共轭函数为:
f
∗
(
z
)
=
sup
x
∈
dom
f
(
z
T
x
−
f
(
x
)
)
f^*(z)=\sup_{x\in \text{dom}f}(z^Tx-f(x))
f∗(z)=x∈domfsup(zTx−f(x))则公式(3)的推导如下
h
(
y
)
=
inf
x
L
(
x
,
y
)
=
−
sup
x
−
L
(
x
,
y
)
=
−
sup
(
−
y
T
A
x
−
f
(
x
)
+
y
T
b
)
=
−
sup
(
−
y
T
A
x
−
f
(
x
)
)
−
y
T
b
=
−
f
∗
(
−
A
T
y
)
−
y
T
b
\begin{aligned} h(y)&=\inf_{x}L(x,y) \\ &=-\sup_x-L(x,y) \\ &=-\sup(-y^TAx-f(x)+y^Tb) \\ &=-\sup(-y^TAx-f(x)) -y^Tb \\ &=-f^*(-A^Ty)-y^Tb \end{aligned}
h(y)=xinfL(x,y)=−xsup−L(x,y)=−sup(−yTAx−f(x)+yTb)=−sup(−yTAx−f(x))−yTb=−f∗(−ATy)−yTb
此时将(1)式的优化问题转化为其对偶问题:
max
h
(
y
)
(4)
\max h(y) \tag{4}
maxh(y)(4)需要注意的是,这里我们假设式(1)和(4)是满足强对偶条件的(如 slater条件,此处不展开介绍),即原问题(1)和对偶问题(4)的解是相等的,即我们可以通过对偶问题的解(最优点)
y
∗
y^*
y∗得到原问题的解(最优点):
x
∗
=
arg min
x
L
(
x
,
y
∗
)
(5)
x^*=\argmin_xL(x,y^*) \tag{5}
x∗=xargminL(x,y∗)(5)
由于对偶问题(4)是max,所以要用梯度上升而不是梯度下降来优化,这就是 “对偶上升” 名字的由来。假设对偶函数
h
h
h是可导的,即
∇
h
(
y
)
\nabla h(y)
∇h(y)存在,那么对偶上升法的迭代更新过程如下:
x
k
+
1
:
=
arg min
x
L
(
x
,
y
k
)
(6)
x^{k+1}:=\argmin_xL(x,y^k) \tag{6}
xk+1:=xargminL(x,yk)(6)
y
k
+
1
:
=
y
k
+
α
k
∇
h
(
y
)
=
y
k
+
α
k
(
A
x
k
+
1
−
b
)
(7)
y^{k+1}:=y^k+\alpha^k\nabla h(y)=y^k+\alpha^k(Ax^{k+1}-b) \tag{7}
yk+1:=yk+αk∇h(y)=yk+αk(Axk+1−b)(7)
其中
α
k
>
0
\alpha^k>0
αk>0表示步长。公式(6)称为
x
x
x-minimization step,公式(7)被称为 dual variable update。这里特别提前一下,这个更新过程能够正常工作的前提是
f
(
x
)
f(x)
f(x)需要满足一些强条件,比如
f
(
x
)
f(x)
f(x)必须要严格凸的或者有界的;再比如如果
f
(
x
)
f(x)
f(x)是线性的即
f
(
x
)
=
k
x
+
b
f(x)=kx+b
f(x)=kx+b类似的函数,那么(6)式是找不到极小值点的,就无法更新
x
x
x了。后面会介绍如何解决这个问题(增广拉格朗日)。
2. 对偶分解(Dual Decomposition)
现在我们假设目标函数
f
f
f是可分的,即
f
(
x
)
f(x)
f(x)可以写成
f
(
x
)
=
∑
i
=
1
N
f
i
(
x
i
)
(8)
f(x)=\sum_{i=1}^Nf_i(x_i) \tag{8}
f(x)=i=1∑Nfi(xi)(8)
其中
x
=
(
x
1
,
.
.
.
,
x
N
)
x=(x_1,...,x_N)
x=(x1,...,xN),
x
i
∈
R
n
i
x_i \in R^{n_i}
xi∈Rni是
x
x
x的子向量(subvector)。这样
f
i
f_i
fi之间相互独立,自变量
x
i
x_i
xi之间也没有关系,我们就可以并行的去计算每一个
f
i
(
x
i
)
f_i(x_i)
fi(xi)的最小值。然后,我们在将
A
A
A按列分块,即
A
x
=
[
A
1
,
.
.
.
,
A
N
]
x
=
A
1
x
i
+
,
.
.
.
,
A
N
x
N
=
∑
i
=
1
N
A
i
x
i
(9)
Ax=[A_1,...,A_N]x=A_1x_i+,...,A_Nx_N=\sum_{i=1}^N A_ix_i \tag{9}
Ax=[A1,...,AN]x=A1xi+,...,ANxN=i=1∑NAixi(9)所以,很自然地,拉格朗日函数也可进行分解,即
L
(
x
,
y
)
=
∑
i
=
1
N
L
i
(
x
i
,
y
)
=
∑
i
=
1
N
(
f
i
(
x
i
)
+
y
T
A
i
x
i
)
−
y
T
b
(10)
L(x,y)=\sum_{i=1}^N L_i(x_i,y)=\sum_{i=1}^N (f_i(x_i)+y^TA_ix_i)-y^Tb \tag{10}
L(x,y)=i=1∑NLi(xi,y)=i=1∑N(fi(xi)+yTAixi)−yTb(10)公式(10)意味着我们在做公式(6)x-minimazation的时候,可以并行来做。即
x
i
k
+
1
:
=
arg min
x
i
L
(
x
i
,
y
k
)
(11)
x_i^{k+1}:=\argmin_{x_i}L(x_i,y^k) \tag{11}
xik+1:=xiargminL(xi,yk)(11)
y
k
+
1
:
=
y
k
+
α
k
∇
h
(
y
)
=
y
k
+
α
k
(
A
x
k
+
1
−
b
)
(12)
y^{k+1}:=y^k+\alpha^k\nabla h(y)=y^k+\alpha^k(Ax^{k+1}-b) \tag{12}
yk+1:=yk+αk∇h(y)=yk+αk(Axk+1−b)(12)
在以上两个公式(11)和(12)的迭代过程中,我们需要“广播”(broadcast)和 “聚合”(gather)操作。通俗的讲,一方面,
A
i
x
i
k
+
1
A_ix_i^{k+1}
Aixik+1的值需要聚合在一起得到
(
A
x
k
+
1
−
b
)
(Ax^{k+1}-b)
(Axk+1−b)(这个也被称为残差residual);另一方面,当得到对偶变量
y
k
+
1
y^{k+1}
yk+1的值后,它得被广播到每一个
x
i
x_i
xi-minimization step中。
3. 增广拉格朗日乘子法(Augmented Lagrangians and the Method of Multipliers)
我们把公式(2)的拉格朗日函数加上一个二次正则项,便将其升级为增广拉格朗日乘子法
L
ρ
(
x
,
y
)
=
f
(
x
)
+
y
T
(
A
x
−
b
)
+
ρ
2
∣
∣
A
x
−
b
∣
∣
2
2
(13)
L_{\rho}(x,y)=f(x)+y^T(Ax-b) +\frac{\rho}{2}||Ax-b||_2^2 \tag{13}
Lρ(x,y)=f(x)+yT(Ax−b)+2ρ∣∣Ax−b∣∣22(13)其中
ρ
>
0
\rho>0
ρ>0称为惩罚系数(penalty parameter)。
现在我们来理解一下这个增广拉格朗日乘子法的作用:
- (1)有了后面这个二次项,我们不再需要假设 f ( x ) f(x) f(x)是严格凸的或者有界的了,比如此时即使 f ( x ) f(x) f(x)是线性的,也可以求导算极值了;
- (2)增加这个二次项后会使我们的求解过程更加“光滑”(smooth),使得迭代过程更加稳定(robust)更好收敛;
- (3)增加这一项不会改变原问题,因为我们受该问题的约束条件限制,最终得到的最优解是满足 A x − b = 0 Ax-b=0 Ax−b=0的,所以和原问题是等价的。
好了,到这里,我们重新书写一下我们的优化问题:
min
f
(
x
)
+
ρ
2
∣
∣
A
x
−
b
∣
∣
2
2
s.t.
A
x
=
b
(14)
\min \quad f(x)+\frac{\rho}{2}||Ax-b||_2^2\\ \text{ s.t.} \quad Ax=b \tag{14}
minf(x)+2ρ∣∣Ax−b∣∣22 s.t.Ax=b(14)
然后,对偶上升法也被重新书写为:
x
k
+
1
:
=
arg min
x
L
ρ
(
x
,
y
k
)
(15)
x^{k+1}:=\argmin_xL_\rho(x,y^k) \tag{15}
xk+1:=xargminLρ(x,yk)(15)
y
k
+
1
:
=
y
k
+
ρ
∇
h
(
y
)
=
y
k
+
ρ
(
A
x
k
+
1
−
b
)
(16)
y^{k+1}:=y^k+\rho\nabla h(y)=y^k+\rho(Ax^{k+1}-b) \tag{16}
yk+1:=yk+ρ∇h(y)=yk+ρ(Axk+1−b)(16)
此时这个步长要求是
ρ
\rho
ρ值
注
2
^{注2}
注2.
注2:这里解释一下为什么步长是
ρ
\rho
ρ。细节可以参考Boyd 《Convex Optimization》中的第5章。
\quad
假设我们得到了满足式(1)的最终优化结果
(
x
∗
,
y
∗
)
(x^*,y^*)
(x∗,y∗),那么将满足可行性(feasibility)条件如下:
A
x
∗
−
b
=
0
,
∇
f
(
x
∗
)
+
A
T
y
∗
=
0
(
这
个
就
是
对
式
(
2
)
求
导
)
Ax^*-b=0, \qquad \nabla f(x^*)+A^Ty^*=0 (这个就是对式(2)求导)
Ax∗−b=0,∇f(x∗)+ATy∗=0(这个就是对式(2)求导)我们令
x
k
+
1
x^{k+1}
xk+1为可以使
L
ρ
(
x
,
y
k
)
L_\rho (x,y^k)
Lρ(x,yk)最小化的值(公式(15)),则有
0
=
∇
x
L
ρ
(
x
k
+
1
,
y
k
)
=
∇
x
f
(
x
k
+
1
)
+
A
T
(
y
k
+
ρ
(
A
x
k
+
1
−
b
)
)
=
∇
x
f
(
x
k
+
1
)
+
A
T
y
k
\begin{aligned} 0&=\nabla_xL_\rho (x^{k+1},y^k) \\ &=\nabla_x f(x^{k+1})+A^T(y^k+\rho(Ax^{k+1}-b)) \\ &=\nabla_x f(x^{k+1})+A^Ty^k \end{aligned}
0=∇xLρ(xk+1,yk)=∇xf(xk+1)+AT(yk+ρ(Axk+1−b))=∇xf(xk+1)+ATyk注意观察后两个等号后面的第二个式子,有
y
k
+
1
=
y
k
+
ρ
(
A
x
k
+
1
−
b
)
y^{k+1}=y^k+\rho(Ax^{k+1}-b)
yk+1=yk+ρ(Axk+1−b),因此这里只有步长取
ρ
\rho
ρ才能使
(
x
k
+
1
,
y
k
+
1
)
(x^{k+1},y^{k+1})
(xk+1,yk+1)对偶可行性(dual feasibility)成立。
注意了,我们刚刚提到增广拉格朗日乘子法给优化过程中的收敛性带来了极大的改善,但这是有代价的。当函数 f f f是可分解的(separable)时,增广拉格朗日函数 L ρ L_\rho Lρ却不可分解了!原因是二次项展开后有会有一个交叉项。总之,此时我们无法将公式(15)像公式(11)那样分解了,无法将任务拆分成多个子任务并行计算了,导致在数据维度很高的情况下无法运行。为了解决这个问题,ADMM终于出场了。
- 需要指出的是,ADMM并不能将原问题分解成多个并行计算的子任务,而是退而求其次,将其转变为可以串行计算的子任务,即先更新 x 1 x_1 x1,完成后才能再更新 x 2 x_2 x2。这样便可以将一个复杂度过高无法求解的问题转化成多个复杂度低的子问题依次求解了。
二、交替方向乘子法(ADMM)
- ADMM将对偶上升法的可分解性与乘子法的收敛性相结合的算法!
ADMM基本算法
为了说明ADMM,我们将原问题即公式(1)的目标函数
f
(
x
)
f(x)
f(x)分解成两个
f
(
x
)
f(x)
f(x)和
g
(
z
)
g(z)
g(z),相应的原变量
x
x
x也分解成两个
(
x
,
z
)
(x,z)
(x,z),则我们得到优化问题
min
f
(
x
)
+
g
(
z
)
s.t.
A
x
+
B
z
=
c
(17)
\min \quad f(x)+g(z)\\ \text{s.t.} Ax+Bz=c \tag{17}
minf(x)+g(z)s.t.Ax+Bz=c(17)为了避免混乱,这里不再说明矩阵
A
A
A,
B
B
B以及向量
x
,
z
,
c
x,z,c
x,z,c的维度。同样,我们假设
f
f
f和
g
g
g是凸的(convex)。
现在我们直接写出公式(17)的增广拉格朗日乘子法的函数
L
ρ
(
x
,
z
,
y
)
=
f
(
x
)
+
g
(
z
)
+
y
T
(
A
x
+
B
z
−
c
)
+
ρ
2
∣
∣
A
x
+
B
z
−
c
∣
∣
2
2
(18)
L_\rho (x,z,y)=f(x)+g(z)+y^T(Ax+Bz-c)+\frac{\rho}{2}||Ax+Bz-c||_2^2 \tag{18}
Lρ(x,z,y)=f(x)+g(z)+yT(Ax+Bz−c)+2ρ∣∣Ax+Bz−c∣∣22(18)同样,我们直接给出ADMM的迭代优化过程
x
k
+
1
:
=
arg min
x
L
ρ
(
x
,
z
k
,
y
k
)
(19)
x^{k+1}:=\argmin_xL_\rho(x,z^k,y^k) \tag{19}
xk+1:=xargminLρ(x,zk,yk)(19)
z
k
+
1
:
=
arg min
z
L
ρ
(
x
k
+
1
,
z
,
y
k
)
(20)
z^{k+1}:=\argmin_zL_\rho(x^{k+1},z,y^k) \tag{20}
zk+1:=zargminLρ(xk+1,z,yk)(20)
y
k
+
1
:
=
y
k
+
ρ
(
A
x
k
+
1
+
B
z
k
+
1
−
c
)
(21)
y^{k+1}:=y^k+\rho(Ax^{k+1}+Bz^{k+1}-c) \tag{21}
yk+1:=yk+ρ(Axk+1+Bzk+1−c)(21)
其中
ρ
>
0
\rho>0
ρ>0。可以看到,ADMM和前面提到的对偶上升和乘子法非常相似,ADMM包含一个
x
x
x-minimization step (19)、一个
z
z
z-minimization step (20)和一个dual variable update (21)。需要注意的是,
x
x
x和
z
z
z的更新是交替的(alternating)。到这里,是不是觉得ADMM其实也没啥:)。
典型的例子
Soft Thresholding
Quatratic Objective
这是两个比较典型ADMM例子,有需要的可以自己去看看《Distributed_Optimization_and_Statistical_Learning》这本书的4.4.3和5.2。笔者这里懒得码字了hhh。
欢迎交流,谢谢指正!