% Generates figures 9.7 and 9.8 of Boyd & Vandenberghe, Convex
% Optimization
%
% Effect of scaling on convergence of gradient method for the problem
%
%     minimize c'*x - sum_i log(bi-ai'*x)  
%
% x is scaled as T*x with 
%
%     T = diag(1, gamma^(1/n), gamma^(2/n), ..., gamma^((n-1)/n)

randn('state',0);
rand('state',0);   
A = randn(500,100);
b = rand(500,1);         
c = 20*randn(100,1);
m = 500;
n = 100;


% Compute minimum using Newton's method.

alpha = 0.01;
beta = 0.5;
maxiters = 50;
x = zeros(100,1);
for k=1:maxiters
    y = b-A*x;
    val = c'*x - sum(log(y));  
    grad = c + A'*(1./y);
    hess = A'*diag(1./(y.^2))*A;
    if (norm(grad) < 1e-8), break; end;
    v = -hess\grad; 
    fprime = grad'*v;
    s = 1;
    dy = -(A*v)./y;
    dc = c'*v;
    for lsiters = 1:100
        newx = x+s*v;
        if (min(1+s*dy) < -1e-5) 
            s = beta*s;
        else 
            newval = val + s*dc - sum(log(1+s*dy));
            if (newval > val + s*alpha*fprime), s = beta*s;
            else break;
            end;
        end;
    end;
    x =  x+s*v;
end;
fmin = val;
xmin = x;


% Figure 9.8.  Condition numbers at optimum vs gamma.

figure(1)
gammas = logspace(-1.4,1.5,1000);
condnrs = zeros(size(gammas));
for k=1:length(gammas);
   gamma = gammas(k);
   T = diag(gamma.^[0:1/n:(n-1)/n]);
   condnrs(k) = cond(T'*hess*T);
end;
loglog(gammas, condnrs, '-');
axis([0.05 20 10 1e4]); 
xlabel('gamma');
ylabel('kappa');


% Figure 9.7.  Number of iterations versus gamma.

gammas  = logspace(-1.5,1.5,11);
iters = zeros(size(gammas));
alpha = 0.3;
beta = 0.7;
maxiters = 5000;
for i = 1:length(gammas)
    gamma = gammas(i);
    T = diag(gamma.^[0:1/n:(n-1)/n]);
    x = zeros(n,1);
    Asc = A*T;
    csc = T'*c;
    for k=1:maxiters
        y = b-Asc*x;
        val = csc'*x - sum(log(y));  
        grad = csc + Asc'*(1./y);
        if (val - fmin < 1e-5), break; end;
        v = -grad; 
        fprime = grad'*v;
        s = 1;
        dy = -(Asc*v)./y;
        dc = csc'*v;
        for lsiters = 1:1000
            if (min(1+s*dy) < -1e-5) 
                s = beta*s;
            else 
                newval = val + s*dc - sum(log(1+s*dy));
                if (newval > val + s*alpha*fprime), s = beta*s;
                else break;
                end;
            end;
        end;
    x =  x+s*v;
    end;
    iters(i) = k;
end;

figure(2)
loglog(gammas, iters, 'o', gammas, iters,'-');
axis([0.05 20 100 5000]); 
xlabel('gamma');
ylabel('iters');