% This code computes a matrix in the neighborhood (not containing zero)
% (Q_ij + eta_ij)>= X_ij >= (Q_ij - eta_ij)
% with minimum condition number, where Q is a positive definite
% n by n matrix and eta is symmetrix. It writes input to SeDuMi to solve
%
% min  lambda_max(X)
% s.t. X - I >= 0, t >= 0, (Q_ij + eta_ij)t>= X_ij >= (Q_ij - eta_ij)t
%% 
% Input
% Q       - a positive definite matrix
% eta     - a symmetric matrix so that the neighborhood does not contain
%           the zero matrix
%%
% Output
% answer  - approximate minimum condition number
% X       - matrix in the convex hull with cond(X) = t
% time    - cpu time
%%
% March 5, 2011
% Written by Zhaosong Lu and Ting Kei Pong

function [answer, time, X] = cond_se_bd_dual2(Q,eta)
% s t X
clear K At b c

n = size(Q,1);

lb = Q - eta;
ub = Q + eta;

At = sparse(1,2+n*(n+1)/2);
c = sparse(1,1);
M = sparse(n^2,n*(n+1)/2);

count = 1;
for i = 1:n
  for j = i:n
    At(2*count-1,[2 2+count]) = [-ub(i,j) 1];
    At(2*count,[2 2+count]) = [lb(i,j) -1];
    c(2*count-1:2*count,1) = [0;0];
    M(n*(i-1)+j,count) = 1;
    M(n*(j-1)+i,count) = 1;
    count = count + 1;
  end
end
At(2*count-1,2) = -1;
count = 2*count;

K.l = n*(n+1)+1;

At(count:count-1+n^2,2+[1:n*(n+1)/2]) = -M;
c(count:count-1+n^2,1) = -reshape(speye(n),n^2,1);

count = count + n^2;

At(count:count-1+n^2,2+[1:n*(n+1)/2]) = M;
At(count:count-1+n^2,1) = -reshape(speye(n),n^2,1);
c(count:count-1+n^2,1) = sparse(n^2,1);

K.s = [n; n];

b = [-1 sparse(1,n*(n+1)/2+1)];

tic
[x,y] = sedumi(At,b,c,K);
time = toc;

t = y(2);
X = reshape(M*y(3:end),n,n);
answer = cond(min(max(X/t,lb),ub));
X = min(max(X/t,lb),ub);




