% Generates figures 11.17 and 11.18 of Boyd & Vandenberghe, Convex % Optimization % % Barrier method for a small SDP. function figure_11_17() % Generate a pair of strictly primal and dual feasible SDPs % % (primal) minimize c'*x % subject to x1*A1 + ... + xn*An <= B % % (dual) maximize -tr(B*Z) % subject to tr(Ai*Z) + ci = 0, i=1,...,m % Z >= 0 randn('state',0); rand('state',0); m = 100; n = 50; % Random symmetric Ai's A = randn(m^2,n); for i=1:n mtrx = reshape(A(:,i),m,m); mtrx = (mtrx + mtrx')/2; A(:,i) = mtrx(:); end; % Random positive definite B B = randn(m,m); B = B*B'; [V,D] = eig(B); B = V*diag(rand(m,1))*V'; % ci = -tr(Ai*Z) for a random positive definite Z Z = randn(m,m); Z = Z*Z'; [V,D] = eig(Z); Z = V*diag(rand(m,1))*V'; c = -A'*Z(:); % Compute optimal vaue and shift B so that optimal value is 1. x0 = zeros(n,1); [x, Z, inniters, gaps] = sdp(A, B, c, x0, 10, 1e-6); s = (1-c'*x)/norm(c)^2; B = B + s*reshape(A*c,m,m);; x0 = s*c; % Figure 11.17. Gap versus Newton iterations for three values of mu. figure(1) muvals = [2, 50, 150]; [x, Z, iters, gaps] = sdp(A, B, c, x0, 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] - iters1(1); gaps1 = [gaps1, gaps(l)]; [x, Z, iters, gaps] = sdp(A, B, c, x0, 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] - iters2(1); gaps2 = [gaps2, gaps(l)]; [x, Z, iters, gaps] = sdp(A, B, c, x0, 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] - iters3(1); gaps3 = [gaps3, gaps(l)]; semilogy(iters1,gaps1, iters2,gaps2, '-', iters3,gaps3, '--'); axis([0 100 0.5e-7 1e3]); xlabel('Newton iterations'); ylabel('duality gap'); 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=150'); % Figure 11.19. Number of Newton iterations versus mu. figure(2) muvals = [1.5 2 4 6 8 10:10:120]; noiters=zeros(1,length(muvals)); for i=1:length(muvals) [x, Z, iters, gaps] = sdp(A, B, c, x0, muvals(i), 1e-3); noiters(i) = iters(length(iters))-iters(1); % skip initial centering disp(['mu = ', num2str(muvals(i)), ': ', num2str(noiters(i)), ... ' Newton iterations.']) end; plot(muvals,noiters,'o', muvals,noiters,'-'); axis([ 0 120 0 140]); xlabel('mu'); ylabel('Newton iterations'); function [x, Z, inniters, gaps] = sdp(A, B, c, x0, mu, tol) % [x, Z, inniters, gaps] = sdp(A, B, c, x0, mu, tol) % % Solves the SDP % % minimize c'*x % subjec to x1*A1 + ... + xn*An <= B % % x0: strictly feasible starting point, not necesarily on central path % B: symmetric nxn matrix % A: [A1(:), A2(:) ... An(:)] % inniters: no of Newton iters per outer iteration % gaps: duality gap at the end of each outer iterations n = length(c); m = size(B,1); MAXITERS = 500; ALPHA = 0.01; BETA = 0.5; NTTOL = 1e-5; % stop inner iteration if lambda^2/2 < NTTOL t = 1; x = x0; gaps = []; inniters = []; Asc = zeros(m^2,1); for k=1:MAXITERS D = B - reshape(A*x,m,m); R = chol(D); Rinv = inv(R); %D = R'*R val = t*c'*x - 2*sum(log(diag(R))); for i=1:n Asc(:,i) = reshape(Rinv'*reshape(A(:,i),m,m)*Rinv, m^2,1); end; g = t*c + Asc'*reshape(eye(m),m^2,1); H = Asc'*Asc; dx = -H\g; fprime = g'*dx; if ((-fprime/2) < NTTOL) Z = Rinv*Rinv'/t; gap = Z(:)'*D(:); gaps = [gaps,gap]; inniters = [inniters, k]; if (gap < tol), return; end; t = min(t*mu, (m+1)/tol); else s = 1; [Rnew, nonpd] = chol(B- reshape(A*(x+s*dx),m,m)); while (nonpd), s = BETA*s; [Rnew,nonpd] = chol(B- reshape(A*(x+s*dx),m,m)); end; while (t*c'*(x+s*dx) - 2*sum(log(diag(Rnew))) > ... val + ALPHA*s*fprime) s = BETA*s; Rnew = chol(B- reshape(A*(x+s*dx),m,m)); end; x = x+s*dx; end; end;