{0,1} Integer Programming; dual ascent

ip-dualascent.m

%
% solve   min c'x s.t. Ax >= b, x in {0,1}^n
% 
% Lagrangian L(x,u) := cx + ub - uAx
%
% original problem:
%   min[x] max[u] L(x,u) 
% Lagrangian dual
%   max[u] min[x] L(x,u)
%
% alternating optimization:
%   Step P. x = argmin L(x,u), u is fixed
%   Step D. u = argmax L(x,u), x is fixed
%

RandStream.setGlobalStream( RandStream('mt19937ar','Seed',1) );

m = 500;
n = 200;
A = randi(2, m, n)-1; % in Map(U -> V)
b = ones(m, 1);       % in V
c = randi(100, n, 1); % in U*


x = randi(2, n,1)-1;  % in U;  x in {0,1}^n
u = ones(m, 1);       % in V*; u >= 0

optx = ones(n,1);
optv = c'*optx;
optvs = [];
for iter = 1 : 1000
  % u = argmax Lt(x,u),
  % L(x,u) = cx + ub - uAx 
  %        = u bhat + const.
  bhat = b - A * x;      % reduced dual cost
  u = u + bhat/iter;     % projected subgradient method 
  u(u < 0) = 0;          %   with O(1/iter) convergence

  % x = argmin Lt(x,u),
  % L(x,u) = cx + ub - uAx
  %        = chat x + const.
  chat = c - A' * u;     % reduced primal cost
  x = double(chat <= 0); % closed form solution

  % save optimal solution
  if A*x >= b
    if c'*x < optv
      optv = c'*x;
      optx = x;
    end
  end
  optvs = [optvs, optv];
end
disp(optv);
plot(optvs);
%disp(optx);