% Generates figure 11.21 of Boyd & Vandenberghe, Convex Optimization
%
% Primal-dual method for a small LP.

function figure_11_21()

% Generate a problem with optimal value 1.

randn('state',0);
rand('state',0);
m = 100;
n = 50;
A = randn(m,n);
b = rand(m,1);
z0 = rand(m,1);
c = -A'*z0;
x = linprog(c, A, b);
s = (1-c'*x)/norm(c)^2;
b = b + s*A*c;
x0 = s*c; 

[x, z, iters, gaps, rs] = lp_pd(A, b, c, x0);

figure(1)
semilogy([0:length(gaps)-1], gaps,'-', [0:length(gaps)-1], gaps,'o');
xlabel('iteration number');  ylabel('surrogate gap');

figure(2)
semilogy([0:length(gaps)-1], rs, '-', [0:length(gaps)-1], rs, 'o');
axis([0 30 1e-16 1e5]);
xlabel('iteration number'); ylabel('residual');


function [x, z, iters, gaps, resdls] = lp_pd(A, b, c, x0)

% [x, z, iters, gaps, resdls] = lp_pd(A, b, c, x0)
%
% Solves 
%
%     minimize    c'*x                 maximize    -b'*z
%     subject to  A*x <= b             subject to  A'*z + c = 0
%                                                  z >= 0 
% 
% using the primal-dual method with simple line search on norm
% of residual.
%
% x0 must be primal feasible.


MAXITERS = 500;  
TOL = 1e-8;       
RESTOL = 1e-8;    
MU = 10;
ALPHA = 0.01;
BETA = 0.5;

[m,n] = size(A);
gaps = [];  resdls = [];
x = x0;  s = b-A*x;  z = 1./s;  
for iters = 1:MAXITERS
    gap = s'*z;  gaps = [gaps, gap];
    res = A'*z  + c;   resdls = [resdls, norm(res)];
    if ((gap < TOL) & (norm(res) < RESTOL)), return;  end;
    tinv = gap/(m*MU);
    sol = -[zeros(n,n), A'; A, diag(-s./z)] \ [A'*z+c; -s+tinv*(1./z)];
    dx = sol(1:n);  dz = sol(n+[1:m]);  ds = -A*dx;

    % backtracking line search
    r = [c+A'*z;  z.*s-tinv];
    step = min(1.0, 0.99/max(-dz./z));
    while (min(s+step*ds) <= 0), step = BETA*step; end;
    newz = z+step*dz;  newx = x+step*dx;  news = s+step*ds;
    newr = [c+A'*newz;  newz.*news-tinv];
    while (norm(newr) > (1-ALPHA*step)*norm(r))
        step = BETA*step;
        newz = z+step*dz;  newx = x+step*dx;  news = s+step*ds;
        newr = [c+A'*newz;  newz.*news-tinv];
    end;
    x = x+step*dx;  z = z +step*dz;  s = b-A*x;
end;