% Generates figure 9.25, Boyd & Vandenberghe, Convex Optimization % % Number of Newton iterations required to minimize self-concordant % functions. warning off % analytic center of Ax <= b, where b>0 randn('seed',0); rand('seed',0); mvalues = [100, 1000, 1000]; nvalues = [ 50, 500, 50]; noproblems = 5; % number of problems for each dimension nostartingpts = 10; % number of starting points for each problem niters_1 = []; initvals_1 = []; niters_2 = []; initvals_2 = []; niters_3 = []; initvals_3 = []; alpha = 0.1; beta = 0.8; for probsize = 1: length(mvalues) m = mvalues(probsize); n = nvalues(probsize); prob = 0; while prob < noproblems A = randn(m,n); b = rand(m,1); % Find exact minimum. x = zeros(n,1); maxiters = 20; unbounded = 1; for iters=1:maxiters y = b-A*x; g = A'*(1./y); H = A'*diag(1./y.^2)*A; v = -H\g; lambda = sqrt(-g'*v); if (lambda < 1e-8), unbounded = 0; break; end; dy = -A*v; s = 1; ysc = dy./y; while (min(1+s*ysc) < 0), s = beta*s; end; f = -sum(log(y)); fnew = f - sum(log(1+s*ysc)); while (fnew > f - alpha*s*lambda^2), s = beta*s; fnew = f - sum(log(1+s*ysc)); end; x = x+s*v; end; if (unbounded == 1) disp('Unbounded.'), continue; % generate another problem else prob = prob+1; end; xmin = x; fmin = -sum(log(b-A*xmin)); y = b-A*xmin; for k = 1:nostartingpts % Generate a random direction v and compute % % smax = sup{s | xmin + sv > 0} % fmax = f(xmin + (1-1e-6)*smax*v). % % Take a starting point x = xmin + s*v such that % % f(x) = fmin + (k/nostartingpoints)*(fmax - fmin). v = randn(n,1); dy = -A*v; inds = find(dy<0); smax = min( -y(inds)./dy(inds) ); maxinitdiff = -sum( log(1+(1-1e-6)*smax*dy./y) ); initdiff = (k/nostartingpts) * maxinitdiff; l = 0; u = smax; s = exp((log(l+1)+log(u+1))/2)-1; f0 = -sum(log(1+s*dy./y)); while ((abs(f0 - initdiff) > 1e-3) & (u-l > 1e-10)) if f0 > initdiff, u=s; else l=s; end; s = exp((log(l+1)+log(u+1))/2)-1; f0 = -sum(log(1+s*dy./y)); end; x = xmin + s*v; if (probsize==1) initvals_1 = [initvals_1, -sum(log(b-A*x)) - fmin]; elseif (probsize==2) initvals_2 = [initvals_2, -sum(log(b-A*x)) - fmin]; else initvals_3 = [initvals_3, -sum(log(b-A*x)) - fmin]; end; % Run Newton's method. maxiters = 100; for iters=1:maxiters y = b-A*x; g = A'*(1./y); H = A'*diag(1./y.^2)*A; v = -H\g; if (g'*v > 0), iters = Inf; break; end; lambda = sqrt(-g'*v); if (lambda < 1e-5), break; end; dy = -A*v; s = 1; ysc = dy./y; while (min(1+s*ysc) < 0), s = beta*s; end; f = -sum(log(y)); fnew = f - sum(log(1+s*ysc)); while (fnew > f - alpha*s*lambda^2), s = beta*s; fnew = f - sum(log(1+s*ysc)); end; x = x+s*v; end; if (iters == maxiters), iters = Inf; end; disp(['Problem size ', int2str(probsize), ', problem ',... int2str(prob), ', starting point ', int2str(k), ... '. Number of Iterations: ', int2str(iters), '.']); if (probsize==1) niters_1 = [niters_1, iters]; elseif (probsize==2) niters_2 = [niters_2, iters]; else niters_3 = [niters_3, iters]; end; end; end; end; d = plot(initvals_1, niters_1, 'o', initvals_2, niters_2, 's', ... initvals_3, niters_3, 'd'); xlabel('f(x0)-pstar'); ylabel('iterations'); axis([0 35 0 26]);