% 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')