% Generates figure 11.19 of Boyd & Vandenberghe, Convex optimization
%
% Barrier method for SDP
%
%     minimize    sum(x)
%     subject to  B + diag(x) >= 0


function figure_11_19()

% Figure 11.19.  Gap versus Newton iterations for three problems.

randn('state',0);
rand('state',0);
mu = 20;

disp('n=50.')
n = 50;  
B = tril(randn(n,n));  
B = B + B' - 2*diag(diag(B));  
B = B/max(abs(eig(B)));
[x, z, iters, gaps] = etp(B, mu, 1e-4);
l = length(gaps);  iters1 = [];  gaps1 = [];
for i=1:l-1
    iters1 = [iters1, iters(i)-1, iters(i+1)-1];  
    gaps1 = [gaps1, gaps(i), gaps(i)]; 
end;
iters1 = [iters1, iters(l)-1] - iters1(1);  
gaps1 = [gaps1, gaps(l)]; 

disp('n=500.')
n = 500;  
B = tril(randn(n,n));  
B = B + B' - 2*diag(diag(B));  
B = B/max(abs(eig(B)));
[x, z, iters, gaps] = etp(B, mu, 1e-4);
l = length(gaps);  iters2 = [];  gaps2 = [];
for i=1:l-1
    iters2 = [iters2, iters(i)-1, iters(i+1)-1];  
    gaps2 = [gaps2, gaps(i), gaps(i)]; 
end;
iters2 = [iters2, iters(l)-1] - iters2(1);  
gaps2 = [gaps2, gaps(l)]; 

disp('n=1000.')
n = 1000;  
B = tril(randn(n,n));  
B = B + B' - 2*diag(diag(B)); 
B = B/max(abs(eig(B)));
[x, z, iters, gaps] = etp(B, mu, 1e-4);
l = length(gaps);  iters3 = [];  gaps3 = [];
for i=1:l-1
    iters3 = [iters3, iters(i)-1, iters(i+1)-1];  
    gaps3 = [gaps3, gaps(i), gaps(i)]; 
end;
iters3 = [iters3, iters(l)-1] - iters3(1);  
gaps3 = [gaps3, gaps(l)]; 

semilogy(iters1, gaps1, iters2, gaps2, '-', iters3, gaps3, '-');
xlabel('Newton iterations');  ylabel('duality gap');
text(iters1(length(iters1)), gaps1(length(iters1)),'n=50');
text(iters2(length(iters2)), gaps2(length(iters2)),'n=500');
text(iters3(length(iters3)), gaps3(length(iters3)),'n=1000');
axis([0 50 1e-5 1e5]);


function [x, Z, iters, gaps] = etp(B,mu,tol)

% [x, Z, iters, gaps] = etp(B, mu, tol)
%
% Solves the SDP
%
%     minimize    sum(x)
%     subject to  B + diag(x) >= 0
%
% B:  nxn symmetric
% iters:  array with number of Newton iterations per outer iteration
% gaps:  array with duality gaps at the end of each outer iteration
%
% The dual problem is
%
%     maximize   -tr(B*Z)
%     subject to diag(Z) = 1 
%                Z = 0
%
%

n = size(B,1);  

MAXITERS = 500;            
ALPHA = 0.01;  
BETA = 0.5;  
NTTOL = 1e-5;               % stop Nt if lambda^2/2 < NTTOL

gaps = [];  iters = [];
x = 1.1*(-min(eig(B)) + 0.1)*ones(n,1);
t = 1;
for k=1:MAXITERS
    X = B+diag(x);  
    R = chol(X);  Rinv = inv(R);  Xinv = Rinv*Rinv';
    val  = t*sum(x) - 2*sum(log(diag(R))); 
    g = t - diag(Xinv);  
    H = Xinv.^2;
    dx = -H\g;  
    fprime = g'*dx;
    if ((-fprime/2) < NTTOL)
        Z = (1/t)*(Xinv - Xinv*diag(dx)*Xinv);
        gap = X(:)'*Z(:);    
        gaps = [gaps, gap];
        iters = [iters, k];  
        if (gap < tol);  return;  end; 
        t = min(t*mu, (n+1)/tol);  
    else
        s = 1;
        dX = Rinv'*diag(dx)*Rinv;
        Xnew = eye(n) + s*dX;  [Rnew,d] = chol(Xnew);
        while (d), 
            s = BETA*s; 
            Xnew = eye(n) + s*dX;
            [Rnew,d] = chol(Xnew);
        end;
        while (t*s*sum(dx) - 2*sum(log(diag(Rnew))) > ALPHA*s*fprime), 
            s = BETA*s; 
            Xnew = eye(n) + s*dX;   [Rnew,d] = chol(Xnew);
        end;
        x = x+s*dx;   
    end;
end;