原问题
min x c T x s . t . A x = b x ≥ 0 (1) \min_x \;c^Tx\\s.t. \;Ax=b\\x\geq 0 \tag{1} xmincTxs.t.Ax=bx≥0(1)
对偶问题
max
y
b
T
y
s
.
t
.
A
T
y
+
s
=
c
s
≥
0
(2)
\max_y\;b^Ty\\s.t.\;A^Ty+s=c\\s\geq 0\tag{2}
ymaxbTys.t.ATy+s=cs≥0(2)
A
∈
R
m
×
n
,
x
∈
R
n
,
s
∈
R
n
,
y
∈
R
m
A \in \R^{m\times n}, x \in \R^{n}, s \in \R^{n}, y \in \R^{m}
A∈Rm×n,x∈Rn,s∈Rn,y∈Rm
本文用交替方向乘子法(Alternating Direction Method of Multipilers)实现线性规划问题的求解器。
首先,写出对偶问题的增广拉格朗日函数,
L
t
(
x
,
y
,
s
)
=
−
b
T
y
+
x
T
(
A
T
y
+
s
−
c
)
+
t
2
∣
∣
A
T
y
+
s
−
c
∣
∣
2
2
s
.
t
.
s
≥
0
L_t(x,y,s) = -b^Ty+x^T(A^Ty+s-c) + \frac{t}{2}||A^Ty+s-c||_2^2 \\ s.t. \quad s \geq 0
Lt(x,y,s)=−bTy+xT(ATy+s−c)+2t∣∣ATy+s−c∣∣22s.t.s≥0
基于ADMM的迭代规则:
y
k
+
1
=
a
r
g
min
y
L
t
(
y
;
x
k
,
s
k
)
s
k
+
1
=
a
r
g
min
s
L
t
(
s
;
x
k
,
y
k
+
1
)
,
s
≥
0
x
k
+
1
=
x
k
+
t
∗
(
A
T
y
k
+
1
+
s
k
+
1
−
c
)
y^{k+1} = arg \min_y L_t(y;x^k,s^k) \\ s^{k+1} = arg \min_s L_t(s;x^k,y^{k+1}), \quad s\geq 0 \\ x^{k+1} = x^{k} +t*(A^Ty^{k+1}+s^{k+1}-c)
yk+1=argyminLt(y;xk,sk)sk+1=argsminLt(s;xk,yk+1),s≥0xk+1=xk+t∗(ATyk+1+sk+1−c)直到结果收敛。针对该问题,其增广拉格朗日函数非常简单,可以写出迭代步的显式形式:
y
k
+
1
=
a
r
g
min
y
L
t
(
y
;
x
k
,
s
k
)
=
−
(
A
A
T
)
−
1
(
(
A
x
k
−
b
)
/
t
+
A
(
s
k
−
c
)
)
y^{k+1} = arg \min_y L_t(y;x^k,s^k) \\ =-(AA^T)^{-1}\bigg((Ax^k-b)/t+A(s^k-c)\bigg)
yk+1=argyminLt(y;xk,sk)=−(AAT)−1((Axk−b)/t+A(sk−c))
s
k
+
1
=
a
r
g
min
s
L
t
(
s
;
x
k
,
y
k
+
1
)
,
s
≥
0
=
m
a
x
(
c
−
A
T
y
k
+
1
−
x
k
t
,
0
)
s^{k+1} = arg \min_s L_t(s;x^k,y^{k+1}), \quad s\geq 0 \\ = max(c-A^Ty^{k+1}-\frac{x^k}{t},0)
sk+1=argsminLt(s;xk,yk+1),s≥0=max(c−ATyk+1−txk,0)
x
k
+
1
=
x
k
+
t
∗
(
A
T
y
k
+
1
+
s
k
+
1
−
c
)
x^{k+1} = x^{k} +t*(A^Ty^{k+1}+s^{k+1}-c)
xk+1=xk+t∗(ATyk+1+sk+1−c)
代码如下:
function [x] = LP_ADMM_dual(c, A, b, opts, x0)
m = size(A,1);
n = size(A,2);
% 随机初始化
y = randn(m,1);
x = x0;
s = randn(n,1);
t = 0.1; % ALM rate
S = A*A';
U = chol(S);
L = U'; %cholesky decomposition: S = L*U = U'*U
err = 1;
x_old = x;
while(err > 1e-6)
y = -U\(L\((A*x-b)/t+A*(s-c))); % 固定 x,s, 更新 y
s = max(c-A'*y-x/t,0); % 固定 y,x, 更新 s
x = x + (A'*y+s-c)*t; % 固定 y,s, 更新 x
err = norm(x-x_old);
x_old = x;
end
end
测试:
% 生成数据
n = 100;
m = 20;
A = rand(m,n);
xs = full(abs(sprandn(n,1,m/n)));
b = A*xs;
y = randn(m,1);
s = rand(n,1).*(xs==0);
c = A'*y + s;
% 计算误差
errfun = @(x1, x2) norm(x1-x2)/(1+norm(x1));
% 标准答案
figure(1);
subplot(2,1,1);
stem(xs,'fill','k-.')
title('exact solu');
% ADMM 求解
opts = [];
tic;
x1 = LP_ADMM_dual(c, A, b, opts, x0);
t1 = toc;
subplot(2,1,2);
stem(x1,'fill','k-.');
title('lp_admm_dual');
fprintf('lp-alm-dual: cpu: %5.2f, err-to-exact: %3.2e\n', t1, errfun(x1, xs));
lp-admm-dual: cpu: 0.08, err-to-exact: 1.17e-04
又快又准!!!