% Generates figure 9.23 of Boyd & Vandenberghe, Convex Optimization
%
% Newton method for a problem with 10000 variables.
% minimize c'*x - sum_i log(bi-ai'*x);
randn('state',0);
rand('state',0);
m = 100000;
n = 10000;
N = 200;
A = [ speye(n,n); -speye(n,n); sprandn(m-2*n,n,0.000001) ];
m1 = (m-2*n)/N; n1 = n/N;
A1 = sprandn(m1,n1,0.05);
A = A + [sparse(2*n,n); kron(speye(N,N),sparse(A1))];
b = [ones(n,1); ones(n,1); rand(m-2*n,1)];
figure(1);
spy(A'*A);
title('sparsity pattern of hessian');
maxiters = 50;
alpha = 0.01;
beta = 0.5;
vals = zeros(1,maxiters);
steps= zeros(1,maxiters);
x = zeros(n,1);
for iter=1:maxiters
y = b-A*x; yinv = 1./y;
val = -sum(log(y)); vals(iter) = val;
grad = A'*yinv;
D = spdiags(yinv.^2, 0, m, m);
hess = A'*D*A;
hess1 = (hess+hess')/2;
v = -hess1\grad;
lambdasq = -grad'*v;
disp([int2str(iter), ': Nt. decr. = ', num2str(sqrt(lambdasq))]);
if (sqrt(lambdasq) < 1e-5), break; end;
fprime = -lambdasq;
s = 1;
dy = -(A*v)./y;
while (min(1+s*dy) < 0), s = beta*s; end;
newx = x+s*v; newval = val - sum(log(1+s*dy));
while (newval > val + s*alpha*fprime)
s = beta*s;
newx = x+s*v;
newval = val - sum(log(1+s*dy));
end;
x = x+s*v;
end;
fmin = vals(iter);
figure(2)
semilogy([0:iter-1], vals(1:iter) - fmin, '-', ...
[0:iter-1], vals(1:iter) - fmin, 'o');
axis([0 20 1e-8 1e5]);
xlabel('k')
ylabel('error')
评论0