ADMM 交替方向乘子法

ADMM交替方向乘子法，适合接触这一类算法的新手来使用。

function [z history] = quadprog（P q r lb ub rho alpha）
% quadprog  Solve standard form box-constrained QP via ADMM
%
% [x history] = quadprog（P q r lb ub rho alpha）
%
% Solves the following problem via ADMM:
%
%   minimize     （1/2）*x‘*P*x + q‘*x + r
%   subject to   lb <= x <= ub
%
% 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（P1）;
%% ADMM solver

x = zeros（n1）;
z = zeros（n1）;
u = zeros（n1）;

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

for k = 1:MAX_ITER

    if k > 1
        x = R \ （R‘ \ （rho*（z - u） - q））;
    else
        R = chol（P + rho*eye（n））;   %% 分解为一个上三角矩阵
        x = R \ （R‘ \ （rho*（z - u） - q））;
    end

    % z-update with relaxation
    zold = z;
    x_hat = alpha*x +（1-alpha）*zold;
    z = min（ub max（lb x_hat + u））;

    % u-update
    u = u + （x_hat - z）;

    % diagnostics reporting termination checks
    history.objval（k）  = objective（P q r x）;

    history.r_norm（k）  = norm（x - z）;
    history.s_norm（k）  = norm（-rho*（z - zold））;

    history.eps_pri（k） = sqrt（n）*ABSTOL + RELTOL*max（norm（x） norm（-z））;
    history.eps_dual（k）= sqrt（n）*ABSTOL + RELTOL*norm（rho*u）;

    if ~QUIET
        fprintf（‘%3d\t%10.4f\t%10.4f\t%10.4f\t%10.4f\t%10.2f\n‘ k ...
            history.r_norm（k） history.eps_pri（k） ...
            history.s_norm（k） history.eps_dual（k） history.objval（k））;
    end

    if （history.r_norm（k） < history.eps_pri（k） && ...
       history.s_norm（k） < history.eps_dual（k））
         break;
    end
end

if ~QUIET
    toc（t_start）;
end
end

function obj = objective（P q r x）
    obj = 0.5*x‘*P*x + q‘*x + r;
end

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       2395  2016-08-11 11:40  ADMM_quadprog\quadprog.m

     文件        930  2016-09-26 21:54  ADMM_quadprog\quadprog_main.m

     目录          0  2016-09-26 21:55  ADMM_quadprog

----------- ---------  ---------- -----  ----

                 3325                    3