% Generates figure 11.22 of Boyd & Vandenberghe, Convex Optimization % % Primal-dual method for a small GP. function figure_11_22() % Generate a problem with optimal value 1 randn('state',0); rand('state',0); m = 100; n = 50; szs = 5*ones(m+1,1); N = sum(szs); m0 = szs(1); % Random postive nu0 with sum(nu0) = 1. nu = rand(N,1); nu(1:szs(1)) = nu(1:szs(1)) / sum(nu(1:szs(1))); % Random A that satisfies A'*nu = 0. A = randn(N,n); A = A - nu*(nu'*A)/(nu'*nu); % Random x0. x0 = randn(n,1); % Random b0. Other bi's are constructed so that f(Ai*x0+bi) + zi = 0 % with random positive initial slacks z. b = zeros(N,1); b(1:m0) = randn(m0,1); z = rand(m,1); ind = m0; for i = 1:m mi = szs(i+1); b(ind+(1:mi)) = -A(ind+(1:mi),:)*x0 - z(i) - log(mi); ind = ind + mi; end % Compute optimal value and shift b0 so that optimal value is 1. [x, z, iters, gaps, resdls] = gp_pd(A, b, szs, x0); optval = log(sum(exp(A(1:m0,:)*x + b(1:m0)))); b(1:m0) = b(1:m0) + 1 - optval; [x, z, iters, gaps, rs] = gp_pd(A, b, szs, 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(resdls)-1], rs, '-', [0:length(rs)-1], rs, 'o'); axis([0 length(resdls) 1e-15 1e5]); xlabel('iteration number'); ylabel('residual'); function [x, z, iters, gaps, resdls] = gp_pd(A, b, szs, x0) % [x, z, iters, gaps, resdls] = gp_pd(A, b, szs, x0) % % Solves the geometric program % % minimize f(y0) % subject to f(yi) <= 0, i=1,...,m % Ai*x+bi = yi, i=0,...,m % % and its dual % % maximize b0'*nu0 + ... + bm'*num + g(nu0, 1) + ... % g(nu1,lambda1) + ,.. + g(num, lambdam) % subject to sum(nu0) = 1 % sum(nui) = lambdai, i=1,...,m % A0'*nu0 + ... + Am'*num = 0 % nui >= 0, i=0,...,m % lambdai >= 0, i=1,...,m % % where f(y) = log(sum(exp(y), g(y,z) = -sum(y.*log(y/z)) % using a primal-dual method with a simplified line search. % Input arguments: % - A: (sum_i n_i)xn matrix; A = [A0; A1; ...; Am] % - b: (sum_i n_i) vector; b = [b0; b1; ...; bm] % - szs: array with row dimensions of Ai and bi % - x0: n-vector; must be strictly feasible MAXITERS = 500; TOL = 1e-8; RESTOL = 1e-8; MU = 10; ALPHA = 0.01; BETA = 0.5; n = length(x0); N = size(A,1); m = length(szs)-1; n0 = szs(1); % indsf is an array with indices of the first rows of each block in A % indsl is an array with indices of the last rows of each block in A % E is a matrix s.t. [sum(y0), sum(y1), ..., sum(ym) ]' = E*y E = sparse(m+1,N); indsf = cumsum(szs)-szs+1; indsl= cumsum(szs); for i=1:m+1 E(i,indsf(i):indsl(i)) = ones(1,szs(i)); end; x = x0; y = A*x+b; f = log(E*exp(y)); s = -f(2:m+1); optval = f(1); z = 1./s; gaps = []; resdls = []; for iters = 1:MAXITERS y = A*x+b; f = log(E*exp(y)); s = -f(2:m+1); optval = f(1); gap = s'*z; res = A' * ( exp(y) .* (E'* ([1;z] ./ (E*exp(y)))) ); w = exp(y) .* (E' * (1./(E*exp(y)))); gradsf = A' * (diag(w)*E'); % [gradf0 gradf1 ... gradfm] gaps = [gaps gap]; resdls = [resdls norm(res)]; if ((gap < TOL) & (norm(res) < RESTOL)), return; end; H = A(indsf(1):indsl(1),:)' * ( diag(w(indsf(1):indsl(1))) - ... w(indsf(1):indsl(1)) * w(indsf(1):indsl(1))') * ... A(indsf(1):indsl(1),:); for i=1:m H = H + z(i) * A(indsf(i+1):indsl(i+1),:)' * ... ( diag(w(indsf(i+1):indsl(i+1))) - ... w(indsf(i+1):indsl(i+1)) * w(indsf(i+1):indsl(i+1))') * ... A(indsf(i+1):indsl(i+1),:); end; tinv = gap/(m*MU); sol = -[ H, gradsf(:,2:m+1); gradsf(:,2:m+1)' diag(-s./z) ] \ ... [ res ; -s + tinv*(1./z) ]; dx = sol(1:n); dz = sol(n+[1:m]); dy = A*dx; r = [res; z.*s-tinv]; step = min(1.0, 0.99/max(-dz./z)); newx = x + step*dx; newy = y + step*dy; news = -log(E(2:m+1,:)*exp(newy)); while (min(news) <= 0.0) step = BETA*step; newx = x + step*dx; newy = y + step*dy; news = -log(E(2:m+1,:)*exp(newy)); end; newz = z + step*dz; newres = A' * ( exp(newy) .* (E'* ([1;newz] ./ (E*exp(newy)))) ); newr = [newres; newz.*news-tinv]; while (norm(newr) > (1-ALPHA*step)*norm(r)) step = BETA*step; newx = x + step*dx; newy = y + step*dy; newz = z + step*dz; news = -log(E(2:m+1,:)*exp(newy)); newres = A' * (exp(newy) .* (E'* ([1;newz] ./ (E*exp(newy))))); newr = [newres; newz.*news-tinv]; end; x = x+step*dx; z = z +step*dz; end;