% Generates figures 11.4 and 11.5 of Boyd & Vandenberghe, Convex % Optimization % % Barrier method for a small LP. function figure_11_4() % Generate a strictly primal and dual feasible LP. % % (primal) minimize c'*x % subject to A*x <= b % % (dual) maximize -b'*z % subject to A'*z + c = 0, z >= 0 randn('state',0); rand('state',0); m = 100; n = 50; A = randn(m,n); b = rand(m,1); c = A'*rand(m,1); % Solve using linprog. x = linprog(c, A, b); optval = c'*x; % Make a change of variables x := x + s*c, so that optimal value is 1. s = (1-optval)/(c'*c); b = b + s*A*c; x0 = s*c; % xc is the point on the central path with barrier parameter t=1. t = 1; x = x0; for k=1:100 d = b-A*x; val = t*c'*x - sum(log(d)); g = t*c + A'*(1./d); H = A'*diag(1./(d.^2))*A; v = -H\g; fprime = g'*v; s = 1; while (min(b-A*(x+s*v)) < 0), s = s/2; end; while (t*c'*(x+s*v) - sum(log(b-A*(x+s*v))) > val + 0.01*s*fprime), s = s/2; end; x = x+s*v; if ((-fprime/2) < 1e-6), break; end; end; xc = x; % Figure 11.4. Barrier method starting at xc, for three values of mu. muvals = [2 50 150]; [x,iters,gaps] = lp(A,b,c,xc,muvals(1),1e-6); l = length(gaps); iters1 = []; gaps1 = []; for i=1:l-1 iters1 = [iters1 iters(i)-1 iters(i+1)-1]; gaps1 = [gaps1 gaps(i) gaps(i)]; end; iters1 = [iters1 iters(l)-1]; gaps1 = [gaps1 gaps(l)]; [x,iters,gaps] = lp(A,b,c,xc,muvals(2),1e-6); l = length(gaps); iters2 = []; gaps2 = []; for i=1:l-1 iters2 = [iters2 iters(i)-1 iters(i+1)-1]; gaps2 = [gaps2 gaps(i) gaps(i)]; end; iters2 = [iters2 iters(l)-1]; gaps2 = [gaps2 gaps(l)]; [x,iters,gaps] = lp(A,b,c,xc,muvals(3),1e-6); l = length(gaps); iters3 = []; gaps3 = []; for i=1:l-1 iters3 = [iters3 iters(i)-1 iters(i+1)-1]; gaps3 = [gaps3 gaps(i) gaps(i)]; end; iters3 = [iters3 iters(l)-1]; gaps3 = [gaps3 gaps(l)]; figure(1) semilogy(iters1,gaps1, iters2,gaps2, '-', iters3,gaps3, '--'); axis([0 90 0.5e-7 1e3]); text(iters1(length(iters1)), gaps1(length(iters1)),'mu=2'); text(iters2(length(iters2)), gaps2(length(iters2)),'mu=50'); text(iters3(length(iters3)), gaps3(length(iters3)),'mu=100'); axis([0 90 0.5e-7 1e3]); xlabel('Newton iterations'); ylabel('duality gap'); % Figure 11.5. Number of Newton iterations vs. mu. muvals = [1.5 2 4 6 8 10:10:200]; noiters=zeros(1,length(muvals)); for i=1:length(muvals) [x,iters,gaps] = lp(A,b,c,xc,muvals(i),1e-3); noiters(i) = iters(length(iters)); end; figure(2) plot(muvals,noiters,'o', muvals,noiters,'-'); axis([ 0 200 0 150]); xlabel('mu'); ylabel('Newton iterations'); function [x, inniters, gaps] = lp(A, b, c, x0, mu, tol) % [x, inniters, gaps] = lp(A, b, c, x0, mu, tol) % % minimize c'*x % subject to A*x <= b % % maximize -b'*z % subject to A'*z + c = 0 % z >= 0 % % Barrier method for solving LP within absolute accuracy tol, starting % with initial t = 1, at a strictly feasible x0. We assume the % problem is strictly dual feasible. % % inniters: array with number of Newton iters per outer iteration % gaps: array with duality gaps at the end of each outer iteration MAXITERS = 500; ALPHA = 0.01; BETA = 0.5; NTTOL = 1e-5; % stop inner iteration if lambda^2/2 < NTTOL [m,n] = size(A); t = 1; x = x0; gaps = []; inniters = []; for k=1:MAXITERS d = b-A*x; val = t*c'*x - sum(log(d)); g = t*c + A'*(1./d); H = A'*diag(1./(d.^2))*A; v = -H\g; fprime = g'*v; s = 1; while (min(b-A*(x+s*v)) < 0), s = BETA*s; end; while (t*c'*(x+s*v) - sum(log(b-A*(x+s*v))) > val + ALPHA*s*fprime), s = BETA*s; end; x = x+s*v; if ((-fprime/2) < NTTOL) inniters = [inniters, k]; z = (1./d) .* (1 + (A*v)./d) / t; gap = (b-A*x)'*z; gaps = [gaps,gap]; if (gap < tol), return; end; t = min(t*mu, (m+1)/tol); end; end; inniters = [inniters, k]; gaps = [gaps, gap]; disp(['Maxiters (', int2str(MAXITERS), ') exceeded.']);