function [bounds, P, info] = jsr_conic_ellipsoid(M,varargin)
% JSR_CONIC_ELLIPSOID Approximates the jsr using ellipsoidal norms.
%
% [BOUNDS, P, INFO] = JSR_CONIC_ELLIPSOID(M)
% returns lower and upper bounds on the jsr of M using ellipsoid
% norms. The set of matrices M must be a cell array of matrices.
%
% [BOUNDS, P, INFO] = JSR_CONIC_ELLIPSOID(M, LOWER, UPPER)
% Uses the parameters UPPER and LOWER as starting bounds for the
% bisection.
%
% [BOUNDS, P, INFO] = JSR_CONIC_ELLIPSOID(M, OPTS)
% Does the same but with value of the parameters defined in the
% structure OPTS. See below and help JSRSETTINGS for available
% parameters.
%
% BOUNDS contains the lower and upper bounds on the JSR
%
% P is the matrix defining the optimal norm
%
% INFO is a structure containing various data about the iterations :
% info.status - = 0 if normal termination, 1 if maxiter
% reached, 2 if stopped in iteration (ctr-c
% or maxTime reached)
% info.elapsedtime
% info.bissec - [lower upper] last bisection interval
% info.niter - number of iterations (bisection)
% info.opts - the structure of options used
%
%
% The field opts.ellips (generated by jsrsettings) can be used to
% tune the method :
%
% ellips.initBounds - = [LOWER UPPER], gives starting bounds for
% bisection method, if not specified uses a
% few iterations of jsr_prod_bruteForce
% ellips.maxiter - maximum number of iterations (for bisection),
% (1000)
% ellips.tol - asked precision for bisection. Stopping
% condition is
% upper-lower < options.ellips.tol*upper
% Influences the asked SeDuMi precision;
% params.eps = options.ellips.tol/10, (1e-8)
% See also JSRSETTINGS
%
% Requires SeDuMi ( http://sedumi.ie.lehigh.edu/ )
%
% REFERENCES
% V.D.Blondel, Y.Nesterov, J.Theys,
% "On the accuracy of the ellipsoid norm approximation
% of the joint spectral radius",
% Linear Algebra and its Applications, 394(1):91-107, 2005
% T.Ando, M.-h. Shih,
% "Simultaneous contractibility",
% SIAM J. Matrix Anal. Appl. 19(2):487-498, 1998
if (nargin > 1)
if (length(varargin) == 2)
options = jsrsettings('ellips.initBounds',[varargin{1} varargin{2}]);
elseif isstruct(varargin{1})
options = varargin{1};
else
warning('Invalid optional arguments, using default parameters !');
options = jsrsettings;
end
else
options = jsrsettings;
end
% logfile opening
close =1;
if (ischar(options.logfile) )
logFile = fopen(options.logfile,'w');
if (logFile == -1)
warning(sprintf('Could not open file %s',options.logfile));
end
elseif isnumeric(options.logfile)
if (options.logfile==0)
logFile = -1;
elseif options.logfile==1
logFile = fopen('log_ellipsoid','w');
if (logFile == -1)
warning('Could not open logfile')
end
else
logFile = options.logfile;
close =0;
end
else
logFile = fopen('log_ellipsoid','w');
if (logFile == -1)
warning('Could not open logfile')
end
end
if (logFile~=-1)
fprintf(logFile,[datestr(now) '\n\n']);
end
msg(logFile,options.verbose>1,'\n \n******** Starting jsr_conic_ellipsoid ******** \n \n')
starttime = cputime;
% Parameters
maxiter = options.ellips.maxiter;
% Initialization
m = length(M); % number of matrices
n = length(M{1}); % dimension of the matrices
opts.disp = 0;
iter = 0;
status = 2;
% Starting bounds
if(isempty(options.ellips.initBounds))
maxdepth = ceil(log(42)/log(m));
optsbrute = jsrsettings('verbose',0,'logfile',0,'bruteforce.maxdepth',maxdepth);
bounds = jsr_prod_bruteForce(M, optsbrute);
msg(logFile,options.verbose>1,'Bounds from bruteforce with depth %2.0f : [%.15g, %.15g]',maxdepth,bounds(1),bounds(2));
upper = bounds(2);
lower = bounds(1);
else
lower = options.ellips.initBounds(1);
upper = options.ellips.initBounds(2);
end
if (upper-lower <= options.ellips.tol*upper),
P = [];
bounds = [lower upper];
msg(logFile,options.verbose>0,'No bisection needed, initial bounds are tight enough');
elapsedtime = cputime - starttime;
info.status = 0;
info.elapsedtime = elapsedtime;
info.niter = iter;
info.opts = options;
return;
end
candidate = (upper+lower) / 2;
% Semidefinite variables
K.s = ones(1, m+1) * n;
NB_VAR = (m+1) * n^2;
% Objective function
c = sparse(NB_VAR, 1);
% Constraints
A = sparse(m*n*(n+1)/2, NB_VAR);
b = sparse(m*n*(n+1)/2, 1);
idx = 0;
for i = 1:m,
for j = 1:n,
for k = j:n,
idx = idx + 1;
A(idx, i*n^2 + (j-1)*n + k) = -1;
A(idx, 1:n^2) = - vec(M{i}(:, k) * M{i}(:, j)')';
A(idx, j*n-n+k) = A(idx, j*n-n+k) + candidate^2;
end
end
end
params.fid = 0;
params.maxiter = 100;
params.eps = options.ellips.tol/10;
params.bigeps = 1e-3;
stop = 0;
% bisection
msg(logFile,options.verbose>1,'\nStarting bisection method');
while (stop == 0 && iter < maxiter)
iter = iter+1;
valid = 0;
if (iter>1000)
if (mod(iter,100)==0)
msg(logFile,options.verbose>0,'\n> Upper bound in (bisection iteration %d): [%.15g, %.15g], testing %.15g',iter, lower, upper, candidate);
end
elseif (iter>150)
if (mod(iter,50)==0)
msg(logFile,options.verbose>0,'\n> Upper bound in (bisection iteration %d): [%.15g, %.15g], testing %.15g',iter, lower, upper, candidate);
end
elseif (iter>20)
if (mod(iter,10)==0)
msg(logFile,options.verbose>0,'\n> Upper bound in (bisection iteration %d): [%.15g, %.15g], testing %.15g',iter, lower, upper, candidate);
end
else
msg(logFile,options.verbose>0,'\n> Upper bound in (bisection iteration %d): [%.15g, %.15g], testing %.15g',iter, lower, upper, candidate);
end
[result dual infosed] = sedumi(A, b, c, K, params);
msg(logFile,options.verbose>1,'> SeDuMi : feasratio = %.9g (%g sec.)', infosed.feasratio, infosed.cpusec);
% Check feasibility
P = mat(result(1:n^2), n);
for i = 1:m,
if (min(eig(candidate^2 * P - M{i}' * P * M{i})) < 0),
valid = -1;
break;
end
end
if (valid == -1),
msg(logFile,options.verbose>1,' Failure - Infeasible solution');
valid = 0;
elseif (infosed.pinf == 1),
msg(logFile,options.verbose>1,' Failure - Infeasible problem');
elseif (infosed.numerr ~= 0),
msg(logFile,options.verbose>1,' Failure - Numerical problems');
elseif (infosed.feasratio < 0.5) || (infosed.feasratio > 2),
msg(' Failure - Invalid feasibility ratio');
else
msg(logFile,options.verbose>1,' Success - Found feasible solution');
valid = 1;
end
% Update the bounds
old = candidate;
if valid,
upper = candidate;
else
lower = candidate;
end
% Update the constraints
if (upper-lower < upper*options.ellips.tol),
candidate = upper;
% End of bisection loop
stop = 1;
else
candidate = (upper+lower) / 2;
end
idx = 0;
for i = 1:m,
for j = 1:n,
for k = j:n,
idx = idx + 1;
A(idx, j*n-n+k) = A(idx, j*n-n+k) - old^2 + candidate^2;
end
end
end
% Save to file option
if (ischar(options.saveinIt))
if(iter==maxiter); status = 1; end;
elapsedtime = cputime-starttime;
info.status = status;
info.elapsedtime = elapsedtime;
info.bissec = [lower upper];
info.niter = iter;
info.opts = options;
tempLowerEllipsEst = max(1/sqrt(n)*upper,1/sqrt(m)*upper);
bounds = [tempLowerEllipsEst upper];
save([options.saveinIt '.mat'],'bounds','info');
end
if ( (cputime-starttime)>=options.maxTime )
stop = 1;
msg(logFile,options.verbose>0,'\nopts.maxTime reached')
end
end
if (iter==maxiter)
status = 1;
elseif (cputime-starttime)>=options.maxTime
status = 2;
else
status = 0;
end
% Post-processing
result = sedumi(A, b, c, K, params);
P = mat(result(1:n^2), n);
lowerEllipsEst = max(1/sqrt(n)*upper,1/sqrt(m)*upper);
msg(logFile,options.verbose>0,'\n> Bounds on the jsr: [%.15g, %.15g]', lowerEllipsEst, upper);
bounds = [lowerEllipsEst, upper];
elapsedtime = cputime - starttime;
msg(logFile,options.verbose>1,'\n End of algorithm after %5.2f s',elapsedtime)
if (logFile~=-1 && close)
fclose(logFile);
end
info.status = status;
info.elapsedtime = elapsedtime;
info.bissec = [lower, upper];
info.niter = iter;
info.opts = options;
