% Generates figures 11.15 and 11.16 of Boyd & Vandenberghe, Convex 
% Optimization
%
% Barrier method  for a small SOCP.

function figure_11_15()

% Generate a pair of strictly primal and dual feasible SOCPs.
%
% (primal) minimize    c'*x
%          subject to  ||Ai*x + bi||_2 <= ci'*x + di, i=1,...,m
%
% (dual)   maximize    b'*z 
%          subject to  A'*z = c
%                      ||zi||_2 <= wi, i=1,...,m
%
% where 
%
% A = [A1; c1'; A2; c2'; ..., Am; cm']
% b = [b1; d1; b2; d2; ..., bm; dm]
% z = [z1; w1; z2; w2; ... ; zm; wm]

randn('state',0);
rand('state',0);
m = 50; 
n = 50;  
szs = 5*ones(m,1);   % row dimensions of the Ai's
N = sum(szs);

% indsi are the indices of ci' in A
% indsf are the indices of the first row of Ai in A 
% E is a matrix such that [sum(y1); ... ; sum(ym) ] = E*y
indsi = cumsum(szs+1);
indsf = cumsum(szs+1)-szs;
E = sparse(m,N+m);
for i=1:m
    E(i,indsf(i):indsi(i)-1) = ones(1,szs(i));
end;

% Ransom x0;  random zi's;  wi = ||zi|| + random positive slack.
x0 = randn(n,1);  
z = randn(N+m,1);  
z(indsi) = sqrt(E*(z.^2)) + rand(m,1);
A = randn(N+m,n);   
c = A'*z;
b = randn(N+m,1);
b(indsi) = 0;
y = A*x0+b;
b(indsi) = sqrt(E*(y.^2)) - y(indsi) + rand(m,1);

% Compute optimal value and shift b0 so that optimal value is 1.
[x, inniters, gaps] = socp(c, A, b, szs, x0, 1e-6, 10);
s = (1-c'*x)/norm(c)^2;
b = b - s*A*c;
x0 = x0+s*c; 


% xc is the point on the central path with barrier parameter t=1.
% To center, call socp() with tol = 2*m+1.

[xc, inniters, gaps] = socp(c, A, b, szs, x0, 2*m+1, 10);

% Figure 11.15.   Gap versus Newton iterastions for three values of mu.

figure(1)
muvals = [2, 50, 200];
[x, iters, gaps] = socp(c, A, b, szs, xc, 1e-6, muvals(1));
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];  gaps1 = [gaps1, gaps(l)]; 

[x, iters, gaps] = socp(c, A, b, szs, xc, 1e-6, muvals(2));
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];  gaps2 = [gaps2, gaps(l)]; 

[x,iters,gaps] = socp(c, A, b, szs, xc, 1e-6, muvals(3));
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];  gaps3 = [gaps3, gaps(l)]; 

semilogy(iters1,gaps1, iters2,gaps2, '-', iters3,gaps3, '--');
text(iters1(length(iters1)), gaps1(length(iters1)),'mu=2');
text(iters2(length(iters2)), gaps2(length(iters2)),'mu=50');
text(iters3(length(iters3)), gaps3(length(iters3)),'mu=200');
axis([0 90  0.5e-7 1e3]);
xlabel('Newton iterations');  ylabel('duality gap');


% Figure 11.16.  Number of Newton iterations versus mu.

muvals = [1.5, 2:2:8, 10:10:200];
noiters=zeros(1,length(muvals));
for i=1:length(muvals)
    [x, iters, gaps] = socp(c, A, b, szs, xc, 1e-3, muvals(i));
    noiters(i) = iters(length(iters));
end;
figure(2)
plot(muvals,noiters,'o', muvals,noiters,'-');
axis([ 0 200 0 150]);
xlabel('mu');  ylabel('Newton iterations');


function [x, inniters, gaps] = socp(c, A, b, szs, x0, tol, mu)

%
% [x, inniters, gaps] = socp(c, A, b, szs, x0, tol, mu)
%
% Solves the SOCP
%
%     minimize    c'*x
%     subject to  ||Ai*x+bi||^2/(ci'*x+di) <= ci'*x+di, i=1,...,m
%                 ci'*x + di >= 0                           
%
% given a strictly primal feasible x0.
%
% A = [A1; c1'; A2'; c2'; ... ; Am; cm']
% b = [b1; d1; b2; d2; ...; bm; dm]
% szs = [n1; n2; ... ; nm] where Ai = nixm
%
% The dual problem is
%
%     maximize    b'*z 
%     subject to  A'*z = c
%                 ||zi|| <= wi
%              
%  with variable z = [z1; w1; z2; w2; ... ; zm; wm].


MAXITERS = 200;  
NTTOL = 1e-5;    % stop centering when Nt decrement <= NTTOL 
ALPHA = 0.01; 
BETA = 0.5;      % parameters for backtracking line search

n = length(x0);  
m = length(szs);   
N = sum(szs);

% indsi are the indices of ci' in A
% indsf are the indices of the first row of Ai in A 
% E is a matrix s.t. [sum(y1);  ... ; sum(ym) ]' = E*y
% E2 is a matrix s.t. [ sum(y1)-w1;  ... ' sum(ym)-wm ]' = E2*y

indsi = cumsum(szs+1);  
indsf = cumsum(szs+1)-szs;
E = sparse(m,N+m);  
E2 = sparse(m,N+m);
for i=1:m
    E(i,indsf(i):indsi(i)-1) = ones(1,szs(i));
    E2(i,indsf(i):indsi(i)) = [ones(1,szs(i)) -1];
end;

x = x0;  
y = A*x+b;  
primalobj = c'*x;  
dualobj = -Inf;  
gap = Inf;
t = 1;  
inniters = [];  gaps = []; 
for iters = 1:MAXITERS
    y = A*x+b;  
    v = y(indsi);  
    s = v.^2 - E*(y.^2);

    % phi = t*c'*x + sum_i f(ci'*x+di, Ai*x+bi) where  
    % f(w,v) = - sum_i log(v^2 - w'*w) 
    % grad f(w,v) = (2/(v^2 - w'*w)) * [w; -v];
    % hess f(w,v) = (2*(v^2-w'*w)*I + 4*[w;-v]*[w;-v]')/(v^2 - w'*w)^2
    val = t*c'*x - sum(log(s));
    d = 2*(y./(E2'*s));  
    g = t*c + A'*d;
    H = A' * (spdiags(2./(E2'*s),0,N+m,N+m) * A);
    for i=1:m
        u = A(indsf(i):indsi(i),:)'*d(indsf(i):indsi(i));
        H = H + u*u';
    end;
    dx = -H\g;
    fprime = g'*dx;
    dy = A*dx;
    s = 1;
    newx = x+s*dx;  newy = A*newx+b;  newv = newy(indsi);
    while (min(newv < 0) | (max(sqrt(E*(newy.^2)) - newv) > 0)) 
        s = BETA*s; 
        newx = x+s*dx;  newy = A*newx+b;  newv = newy(indsi);
    end;
    newval = t*c'*newx - sum(log(newv)) - ...
        sum(log(newv- (E*newy.^2) ./ newv));
    while (newval > val + s*ALPHA*fprime)
        s = BETA*s; 
        newx = x+s*dx;  newy = A*newx+b;  newv = newy(indsi);
        newval = t*c'*newx - sum(log(newv)) - ...
            sum(log(newv- (E*newy.^2) ./ newv));
    end;
    x = x+s*dx;
    if (-fprime/2 < NTTOL), 
        gap = 2*m/t; 
        gaps = [gaps, gap];
        inniters = [inniters, iters]; 
        if (gap < tol), return; end;
        t = min(t*mu, (2*m+1)/tol);  
    end; 
end;