matlab增广拉格朗日,[Opt] 拉格朗日乘子法 | ADMM | Basis pursuit

乘子法

本文先简要介绍三个乘子法,它们的收敛条件依次减弱(不做具体介绍),然后应用 ADMM 算法求解 Basis pursuit 问题最后读读 Boyd 给出的代码。

Lagrange Multiplier(拉格朗日乘子法)

Augmented Lagrangian Multiplier(增广拉格朗日乘子法,ALM)

Alternating Direction Method of Multipliers(变方向乘子法,ADMM)

乘子法部分参考文献[1]。

拉格朗日乘子法

考虑如下等式约束的优化问题

equation?tex=%5Cmin_%7Bx%7D+f%28x%29%5C%5C+%5Ctext%7Bs.t.%7D+Ax%3Db

用拉格朗日乘子法求解,先写出拉格朗日函数:

equation?tex=L%28x%2Cy%29%3Df%28x%29%2By%5ET%28Ax-b%29%5C%5C

其中

equation?tex=y 是拉格朗日乘子。它的对偶函数为

equation?tex=g%28y%29%3D%5Cinf_x+L%28x%2Cy%29%5C%5C

对偶问题为

equation?tex=%5Cmax_y+g%28y%29%5C%5C

我们需要的解是

equation?tex=x%5E%5Cast%3D%5Carg%5Cmin_x+L%28x%2Cy%5E%5Cast%29%5C%5C

这里的

equation?tex=y%5E%5Cast 是对偶问题的最优解,下面我们用梯度下降法求解它。迭代公式为

equation?tex=y%5E%7Bk%2B1%7D%3Dy%5E%7Bk%7D%2B%5Calpha%5Ek+%5Cnabla+g%28y%5Ek%29%5C%5C

其中

equation?tex=%5Cnabla+g%28y%5Ek%29%3DA%5Ctilde%7Bx%7D-b%2C%5Ctilde%7Bx%7D%3D%5Carg%5Cmin_xL%28x%2Cy%5Ek%29

①总结拉格朗日乘子法的更新方法如下:

equation?tex=%5Cbegin%7Bsplit%7D+x%5E%7Bk%2B1%7D%26%3A%3D%5Carg%5Cmin_xL%28x%2Cy%5Ek%29%5C%5C+y%5E%7Bk%2B1%7D%26%3A%3Dy%5Ek+%2B+%5Calpha%5Ek%28Ax%5E%7Bk%2B1%7D-b%29+%5Cend%7Bsplit%7D%5C%5C

增广拉格朗日乘子法

用增广拉格朗日乘子法求解,先写出增广拉格朗日函数:

equation?tex=L_%5Crho%28x%2Cy%29%3Df%28x%29%2By%5ET%28Ax-b%29%2B%5Cfrac%7B%5Crho%7D%7B2%7D%5C%7CAx-b%5C%7C%5E2_2%5C%5C

②类似拉格朗日乘子法,增广拉格朗日乘子法的更新方法如下:

equation?tex=%5Cbegin%7Bsplit%7D+x%5E%7Bk%2B1%7D%26%3A%3D%5Carg%5Cmin_xL_%5Crho%28x%2Cy%5Ek%29%5C%5C+y%5E%7Bk%2B1%7D%26%3A%3Dy%5Ek+%2B+%5Crho%28Ax%5E%7Bk%2B1%7D-b%29+%5Cend%7Bsplit%7D%5C%5C

ADMM

ADMM 问题形式(with

equation?tex=f%2Cg convex)

equation?tex=%5Cmin_%7Bx%2Cz%7D+f%28x%29%2Bg%28z%29%5C%5C+%5Ctext%7Bs.t.%7D+Ax%2BBz%3Dc

它的增广拉格朗日函数为

equation?tex=L_%5Crho%28x%2Cz%2Cy%29%3Df%28x%29%2Bg%28z%29%2By%5ET%28Ax%2BBz-c%29%2B%5Cfrac%7B%5Crho%7D%7B2%7D%5C%7CAx%2BBz-c%5C%7C%5E2_2%5C%5C

③ADMM 的更新方法如下:

equation?tex=%5Cbegin%7Bsplit%7D+x%5E%7Bk%2B1%7D%26%3A%3D%5Carg%5Cmin_xL_%5Crho%28x%2Cz%5Ek%2Cy%5Ek%29%5C%5C+z%5E%7Bk%2B1%7D%26%3A%3D%5Carg%5Cmin_zL_%5Crho%28x%5E%7Bk%2B1%7D%2Cz%2Cy%5Ek%29%5C%5C+y%5E%7Bk%2B1%7D%26%3A%3Dy%5Ek+%2B+%5Crho%28Ax%5E%7Bk%2B1%7D%2BBz%5E%7Bk%2B1%7D-c%29+%5Cend%7Bsplit%7D%5C%5C

Basis pursuit

优化目标:

equation?tex=%5Cmin_x%5C%7Cx%5C%7C_1%5C%5C+%5Ctext%7Bs.t.%7D+Ax%3Db

写成 ADMM 的形式

equation?tex=%5Cmin_%7Bx%2Cz%7D+f%28x%29%2B%5C%7Cz%5C%7C_1%5C%5C+%5Ctext%7Bs.t.%7D+x%3Dz

其中

equation?tex=f%28x%29 表示指示函数,相关的内容我在 [Opt] 近端最小化算法 有写过。那么它的增广拉格朗日函数为

equation?tex=L_%5Crho%28x%2Cz%2Cy%29%3D+%5Cbegin%7Bcases%7D+%5C%7Cz%5C%7C_1%2By%5ET%28x-z%29%2B%5Cfrac%7B%5Crho%7D%7B2%7D%5C%7Cx-z%5C%7C%5E2_2+%26+%5Ctext%7Bif+%7D+Ax%3Db%2C%5C%5C++%5Cinfty+%26+%5Ctext%7Bif+%7D+Ax%5Cneq+b.+%5Cend%7Bcases%7D%5C%5C

第一步:更新

equation?tex=x :

equation?tex=%5Cbegin%7Bsplit%7D+x_%7Bk%2B1%7D%26%5Cin+%5Carg+%5Cmin_%7BAx%3Db%7D%5C%7B%28y%5Ek%29%5ET%28x-z%5Ek%29%2B%5Cfrac%7B%5Crho%7D%7B2%7D%5C%7Cx-z%5Ek%5C%7C%5E2_2%5C%7D%5C%5C+%26%5Cin%5Carg+%5Cmin_%7BAx%3Db%7D%5C%7B%5Cfrac%7B%5Crho%7D%7B2%7D%5C%7Cx-%28z%5Ek-%5Cfrac%7By%5Ek%7D%7B%5Crho%7D%29%5C%7C_2%5E2%5C%7D+%5Cend%7Bsplit%7D%5C%5C

我们记

equation?tex=u%5Ek%3D%5Cfrac%7By%5Ek%7D%7B%5Crho%7D 。上述问题为 

equation?tex=z%5Ek-u%5Ek 在仿射集上的投影问题,求解参考[4],具体原理等内容我在 [ML] 深入线性回归 提到过。

它的解为

equation?tex=P%28x%29+%3D+x+-+A%5ET%28AA%5ET%29%5E%7B-1%7D%28Ax-b%29%2Cx%3Dz%5Ek-u%5Ek 。

第二步:更新

equation?tex=z :

equation?tex=%5Cbegin%7Bsplit%7D+z%5E%7Bk%2B1%7D%26%5Cin%5Carg%5Cmin_%7Bz%5Cin%5Cmathbb%7BR%7D%5En%7D%5C%7Cz%5C%7C_1%2B%28y%5E%7Bk%7D%29%5ET%28x%5E%7Bk%2B1%7D-z%29%2B%5Cfrac%7B%5Crho%7D%7B2%7D%5C%7Cx%5E%7Bk%2B1%7D-z%5C%7C%5E2_2%5C%5C+%26%5Cin%5Carg%5Cmin_%7Bz%5Cin%5Cmathbb%7BR%7D%5En%7D%5C%7Cz%5C%7C_1%2B%5Cfrac%7B%5Crho%7D%7B2%7D%5C%7Cz-%28x%5E%7Bk%2B1%7D%2Bu%5Ek%29%5C%7C%5E2_2+%5Cend%7Bsplit%7D%5C%5C

equation?tex=z%5E%7Bt%2B1%7D%3DS_%7B%5Cfrac%7B1%7D%7B%5Crho%7D%7D%28x%5E%7Bk%2B1%7D%2Bu%5Ek%29%2CS_%7B%5Cfrac%7B1%7D%7B%5Crho%7D%7D%28x%29%3D%5Ctext%7Bsgn%7D%28x%29%5Cmax%5C%7B%7Cx%7C-%5Cfrac%7B1%7D%7B%5Crho%7D%5C%7D%5C%5C

第三步,更新乘子

equation?tex=y

equation?tex=y%5E%7Bk%2B1%7D%3Dy%5Ek+%2B+%5Crho%28x%5E%7Bk%2B1%7D-z%5E%7Bk%2B1%7D%29%5C%5C

因此我们只要按如下等价地更新

equation?tex=u 即可

equation?tex=u%5E%7Bk%2B1%7D%3Du%5Ek%2Bx%5E%7Bk%2B1%7D-z%5E%7Bk%2B1%7D%5C%5C

④ 总结 Basis pursuit 更新方法如下:

equation?tex=%5Cbegin%7Bsplit%7D++%5Cend%7Bsplit%7D

equation?tex=%5Cbegin%7Bsplit%7D+x%5E%7Bk%2B1%7D%26%3DP%28z%5Ek-u%5Ek%29%2C%5C%5C+z%5E%7Bk%2B1%7D%26%3DS_%7B%5Cfrac%7B1%7D%7B%5Crho%7D%7D%28x%5E%7Bk%2B1%7D%2Bu%5Ek%29%2C%5C%5C+u%5E%7Bk%2B1%7D%26%3Du%5Ek%2Bx%5E%7Bk%2B1%7D-z%5E%7Bk%2B1%7D.++%5Cend%7Bsplit%7D%5C%5C

MATLAB 代码

运行结果

左实线:目标函数值迭代更新变化曲线

右上实线:primal residual

equation?tex=%5C%7Cx-z%5C%7C_2

右上虚线:primal feasibility tolerance

右下实线:dual residual

equation?tex=-%5Crho%5C%7Cz%5E%7Bk%2B1%7D-z%5E%7Bk%7D%5C%7C_2

右下虚线:dual feasibility tolerance

317f52d653b5ebc059544a3825e94c3b.png

MATLAB 中在命令行窗口输入如下指令获取 basis_pursuit 函数

grabcode("http://web.stanford.edu/~boyd/papers/admm/basis_pursuit/basis_pursuit.html")

输入以下指令获取 Basis pursuit example 测试代码

grabcode("http://web.stanford.edu/~boyd/papers/admm/basis_pursuit/basis_pursuit_example.html")

basis_pursuit.m

function[z, history] =basis_pursuit(A, b, rho, alpha)% basis_pursuit Solve basis pursuit via ADMM

%

% [x, history] = basis_pursuit(A, b, rho, alpha)

%

% Solves the following problem via ADMM:

%

% minimize ||x||_1

% subject to Ax = b

%

% The solution is returned in the vector x.

%

% history is a structure that contains the objective value, the primal and

% dual residual norms, and the tolerances for the primal and dual residual

% norms at each iteration.

%

% rho is the augmented Lagrangian parameter.

%

% alpha is the over-relaxation parameter (typical values for alpha are

% between 1.0 and 1.8).

%

% More information can be found in the paper linked at:

% http://www.stanford.edu/~boyd/papers/distr_opt_stat_learning_admm.html

%

t_start = tic;

%% Global constants and defaults

QUIET = 0;

MAX_ITER = 1000;

ABSTOL = 1e-4; % 绝对误差容限

RELTOL = 1e-2; % 相对误差容限

%% Data preprocessing

[~, n] = size(A);

%% ADMM solver

x = zeros(n,1);

z = zeros(n,1);

u = zeros(n,1);

if ~QUIET

fprintf(‘%3s\t%10s\t%10s\t%10s\t%10s\t%10s\n‘, ‘iter‘, ...

‘r norm‘, ‘eps pri‘, ‘s norm‘, ‘eps dual‘, ‘objective‘);

end

% precompute static variables for x-update (projection on to Ax=b)

AAt = A*A‘;

P = eye(n) - A‘ * (AAt \ A);

q = A‘ * (AAt \ b);

for k = 1:MAX_ITER

% x-update

x = P*(z - u) + q;

% z-update with relaxation

zold = z;

x_hat = alpha*x + (1 - alpha)*zold; % 带松弛的更新方法

z = shrinkage(x_hat + u, 1/rho);

u = u + (x_hat - z);

% diagnostics, reporting, termination checks

history.objval(k) = objective(A, b, x);

history.r_norm(k) = norm(x - z); % primal residual

history.s_norm(k) = norm(-

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
©️2022 CSDN 皮肤主题:1024 设计师:我叫白小胖 返回首页
评论
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值