% Generates figure 10.5 of Boyd & Vandenberghe, Convex Optimization % % Newton's method applied to a convex-concave game % % f(x,y) = x'*A*y + c'*x + d'*y - log(1-x'*x) + log(1-y'*y); randn('state',0); n = 100; A = randn(n,n); c = randn(n,1); d = randn(n,1); BETA = .5; ALPHA = .01; x = .01*ones(n,1); y = .01*ones(n,1); resx = []; resy = []; res = []; % Infeasible start Newton method for solving r(x,y) = 0 where % % r(x,y) = ( A*y + c + 2*x/(1-x'*x), A'*x + d - 2*y/(1-y'*y) ) for iters =1:100 r = [ A*y + (2/(1-x'*x))*x + c; A'*x - (2/(1-y'*y))*y + d ]; resx = [resx, norm(r(1:n))]; resy = [resy, norm(r(n+[1:n]))]; res = [res, norm(r)]; if (norm(r) < 1e-8), break; end; Dr = [ ((2/(1-x'*x))*eye(n) + (4/(1-x'*x)^2)*x*x'), A ; A', (-(2/(1-y'*y))*eye(n) - (4/(1-y'*y)^2)*y*y') ]; step = -Dr\r; dx = step(1:n); dy = step(n+[1:n]); s = 1; newx = x+s*dx; newy = y+s*dy; while ((norm(newx) >= 1) | (norm(newy) >= 1)), s = BETA*s; newx = x+s*dx; newy = y+s*dy; end; newr = [ A*newy + (2/(1-newx'*newx))*newx + c; A'*newx - (2/(1-newy'*newy))*newy + d ]; while (norm(newr) > (1-ALPHA*s)*norm(r)) s = BETA*s; newx = x+s*dx; newy = y+s*dy; newr = [ A*newy + (2/(1-newx'*newx))*newx + c; A'*newx - (2/(1-newy'*newy))*newy + d]; end; x = x+s*dx; y = y+s*dy; end; semilogy(0:(length(res)-1), res, 'o', 0:(length(res)-1), res, '-'); xlabel('iteration'); ylabel('residual'); axis([0 8 1e-15 1e5]);