% Generates figures 10.1, 10.2, 10.3, 10.4 of Boyd & Vandenberghe,
% Convex Optimization
%
% Infeasible start Newton method for solving
%
%     minimize   -sum_i log x_i
%     subject to Ax = b


% Generate a feasible problem.

rand('state',3);   
randn('state',4);
n = 100;        
m = 50;        
A = randn(m,n);
x0 = rand(n,1);
b = A*x0;        


% Apply infeasible Newton method to KKT conditions
%
%     -1./x + A'*w = 0 
%      A*x = b

alpha = 0.011;
beta  = 0.5;
tol = 1e-12;
maxiters = 20;

x = ones(n,1);  
w = zeros(m,1); 
resp = zeros(1, maxiters);
resd = zeros(1, maxiters);
steps = zeros(1, maxiters);
for k=1:maxiters
    rp = A*x-b;  rd = -1./x + A'*w;  
    resp(k) = norm(rp);  resd(k) = norm(rd);
    normr = norm([rd; rp]);
    if normr < 1e-12, break; end;
    dw = - (A * diag(x.^2) * A') \ (A*((x.^2).*rd) - rp); 
    dx = - (x.^2) .* (rd + A'*dw);
    s = 1;
    while min(x+s*dx) <= 0,  s = s*beta; end;
    while norm([ -1./(x+s*dx) + A'*(w+s*dw);  A*(x+s*dx)-b ]) > ...
        (1-alpha*s) * normr
        s = s*beta;
    end;
    steps(k) = s;
    x = x + s*dx;
    w = w + s*dw;
end;

% Figure 10.1
figure(1)
semilogy([0:k-1], resp(1:k), 'o', [0:k-1], resp(1:k), '-', ...
    [0:k-1], resd(1:k), 'o', [0:k-1], resd(1:k), '--');
xlabel('iteration number')
ylabel('residuals')

% Figure 10.2 
figure(2);
plot([1:k-1], steps(1:k-1), 'o', [1:k-1], steps(1:k-1), '-');
xlabel('iteration number')
ylabel('step length')
axis([0 12 0.4 1.1]);


% Generate an infeasible problem.

rand('state',3);
randn('state',4);
A = randn(m,n);
x0 = randn(n,1); 
b = A*x0;       


% Apply infeasible Newton method to KKT conditions
%
%     -1./x + A'*w = 0 
%      A*x = b

alpha = 0.011;
beta  = 0.5;
tol = 1e-12;
maxiters = 21;

x = ones(n,1);  
w = zeros(m,1); 
resp = zeros(1, maxiters);
resd = zeros(1, maxiters);
steps = zeros(1, maxiters);
for k=1:maxiters
    rp = A*x-b;  rd = -1./x + A'*w;  
    resp(k) = norm(rp);  resd(k) = norm(rd);
    normr = norm([rd; rp]);
    if normr < 1e-12, break; end;
    dw = - (A * diag(x.^2) * A') \ (A*((x.^2).*rd) - rp); 
    dx = - (x.^2) .* (rd + A'*dw);
    s = 1;
    while min(x+s*dx) <= 0,  s = s*beta; end;
    while norm([ -1./(x+s*dx) + A'*(w+s*dw);  A*(x+s*dx)-b ]) > ...
        (1-alpha*s) * normr
        s = s*beta;
    end;
    steps(k) = s;
    x = x + s*dx;
    w = w + s*dw;
end;

% Figure 10.3
figure(3)
semilogy([0:k-1], resp(1:k), 'o', [0:k-1], resp(1:k), '-', ...
    [0:k-1], resd(1:k), 'o', [0:k-1], resd(1:k), '--');
xlabel('iteration number')
ylabel('residuals')

% Figure 10.4 
figure(4);
plot([1:k-1], steps(1:k-1), 'o', [1:k-1], steps(1:k-1), '-');
xlabel('iteration number')
ylabel('step length')
axis([0 20 0 0.3]);