Code covered by the BSD License

# The JSR toolbox

### Raphael Jungers (view profile)

10 Oct 2011 (Updated )

Gathers and compares the best methods for the joint spectral radius computation

jsr_norm_maxRelaxation2D(M, varargin)
function [bounds, R, Phi, info] = jsr_norm_maxRelaxation2D(M, varargin)

% JSR_NORM_MAXRELAXATION2D Approximates the jsr using Max-Relaxation in 2D.
%
%    [BOUNDS, R, PHI, INFO] = JSR_NORM_MAXRELAXATION2D(M)
%      returns lower and upper bounds on the jsr of M using Max-Relaxation
%      scheme (heuristically). M must be a cell array containing the 2x2
%      matrices. Uses default values (see below) for all the parameters.
%
%    [BOUNDS, R, PHI, INFO] = JSR_NORM_MAXRELAXATION2D(M, STEP, AVGFUN)
%      The parameter STEP corresponds to the discretization step of the
%      angular interval [-pi, +pi]. Default is 2*pi/10000.
%      The parameter AVGFUN is either a handler to an averaging function with
%      2 inputs and 1 output, or 'a', 'g', 'h', 'q' for arithmetic,
%      geometric, harmonic or quadratic mean. Default is 'a'. If not
%      specified, default is 'a'.
%
%    [BOUNDS, R, PHI, INFO] = JSR_NORM_MAXRELAXATION2D(M, OPTS)
%      Does the same but with the values of the parameters in the structure
%      OPTS. See below and help JSRSETTINGS for available parameters.
%
%      BOUNDS contains the lower and upper bounds on the JSR
%
%      The vectors R and PHI are such that we have Phi = -pi:step:+pi and
%      R(i) = ||(cos Phi(i), sin Phi(i))||.
%
%      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
%                               (CTRL-C or opts.maxTime reached)
%         info.elapsedtime
%         info.niter          - number of iterations
%         info.timeIter       - cputtime - start time in s. at the end of
%                               each iteration. Note : diff(info.timeIter)
%                               gives the time taken by each iteration.
%         info.allLb          - evolution of the lower bound  [1x(niter) double]
%         info.allUb          - evolution of the upper bound  [1x(niter) double]
%         info.opts           - the structure of options used
%
%
%  The field opts.maxrel2D (generated by jsrsettings) can be used to
%  tune the method :
%
%      maxrel2D.step           - discretization step of the angular interval
%                                [-pi, +pi]. (2*pi/10000)
%      maxrel2D.avgfun         - either a handler to an averaging function with
%                                2 inputs and 1 output, either 'a', 'g', 'h', 'q'
%                                for arithmetic, geometric, harmonic or
%      maxrel2D.maxiter        - maximum number of iterations, (1000)
%
%      maxrel2D.tol            - tolerance for stopping condition, stops when
%                                upperBound-lowerBound < tol, (1e-10)
%
%      maxrel2D.plotBounds     - if 1 plots the evolution of the bounds, (0)
%
%      maxrel2D.plotEllips     - if 1 plots the unit ball, (0)
%
%      maxrel2D.plotTime       - if 1 plots evolution of the bounds w.r.t.
%                                time of computation, (0)
%
%
%
% REMARK : This program returns heuristic estimates based on [1]
%
% REFERENCES
% [1] V.Kozyakin,
%       "Iterative building of Barabanov norms and computation of the
%        joint spectral radius for matrix sets",
%       arXiv:0810.2154v2 [math.RA], 22 Oct. 2008

if (nargin > 1)
if (length(varargin) > 1)
opts = jsrsettings('maxrel2D.step',varargin{1},'maxrel2D.avgfun',varargin{2});
elseif isnumeric(varargin{1})
opts = jsrsettings('maxrel2D.step',varargin{1});
else
opts = varargin{1};
end
else
opts = jsrsettings;
end

avgfun = opts.maxrel2D.avgfun;
if ischar(avgfun),
switch avgfun(1),
case 'g',  avgfun = @geomMean;
case 'h',  avgfun = @harmMean;
otherwise, avgfun = @arithMean;
end
end

% logfile opening
close =1;
if (ischar(opts.logfile) )
logFile = fopen(opts.logfile,'wt');
if (logFile == -1)
warning(sprintf('Could not open file %s',opts.logfile));
end
elseif isnumeric(opts.logfile)
if (opts.logfile==0)
logFile = -1;
elseif (opts.logfile==1)
logFile = fopen('log_maxRelaxation2D','wt');
if (logFile == -1)
warning('Could not open logfile')
end
else
logFile = opts.logfile;
close =0;
end
else
logFile = fopen('log_maxRelaxation2D','wt');
if (logFile == -1)
warning('Could not open logfile')
end
end

if (logFile~=-1)
fprintf(logFile,[datestr(now) '\n\n']);
end

msg(logFile,opts.verbose>1,'\n \n******** Starting jsr_norm_maxRelaxation2D ******** \n \n')

starttime = cputime;

% Parameters
step = opts.maxrel2D.step;
tol = opts.maxrel2D.tol;
maxiter = opts.maxrel2D.maxiter;
m = length(M);

% Initialization
status = 2;
lowerBounds = zeros(1,maxiter);
upperBounds = zeros(1,maxiter);
timeIter = zeros(1,maxiter);

% Discretization
n = 2*ceil(round(2*pi/step)/2);
Phi = linspace(-pi, pi, n+1);
R = ones(1, n+1);
keyZERO = n/2 + 1;

H = zeros(m, n+1);
PHI = zeros(m, n+1);
keyPhi = zeros(m, n+1);
msg(logFile,opts.verbose>1,'Starting to compute the images of the points');
for i = 1:m,
H(i, :) = sqrt( (M{i}(1,1)*cos(Phi) + M{i}(1,2)*sin(Phi)).^2 + (M{i}(2,1)*cos(Phi) + M{i}(2,2)*sin(Phi)).^2 );
PHI(i, :) = atan2(  M{i}(2,1)*cos(Phi) + M{i}(2,2)*sin(Phi), M{i}(1,1)*cos(Phi) + M{i}(1,2)*sin(Phi)  );
keyPhi(i, :) = round((PHI(i, :) + pi)*n/(2*pi) + 1);
end

% Iteration
msg(logFile,opts.verbose>1,'Starting iteration \n');
for iter = 1:maxiter,
Rall = H .* R(keyPhi);
Rstar = max(Rall, [], 1);
RHOup = max(Rstar./R);
RHOdown = min(Rstar./R);
invGamma = 1/avgfun(RHOup, RHOdown);
R = max(R, Rstar * invGamma);
R = R / R(keyZERO);

lowerBounds(iter) = RHOdown;
upperBounds(iter) = RHOup;
timeIter(iter) = cputime-starttime;

if (iter>1000)
if (mod(iter,100)==0)
msg(logFile,opts.verbose>0,'> Iteration %3.0f - current bounds: [%.15g, %.15g] \n', iter, RHOdown, RHOup);
end
elseif (iter>150)
if (mod(iter,50)==0)
msg(logFile,opts.verbose>0,'> Iteration %3.0f - current bounds: [%.15g, %.15g] \n', iter, RHOdown, RHOup);
end
elseif (iter>20)
if (mod(iter,10)==0)
msg(logFile,opts.verbose>0,'> Iteration %3.0f - current bounds: [%.15g, %.15g] \n', iter, RHOdown, RHOup);
end
else
msg(logFile,opts.verbose>0,'> Iteration %3.0f - current bounds: [%.15g, %.15g] \n', iter, RHOdown, RHOup);
end

% Save to file option
if (ischar(opts.saveinIt))
if (iter==maxiter);status=1;end
bounds = [RHOdown, RHOup];
elapsedtime = cputime - starttime;

info.status = status;
info.elapsedtime = elapsedtime;
info.niter = iter;
info.timeIter = timeIter(1:iter);
info.allLb = lowerBounds(1:iter);
info.allUb = upperBounds(1:iter);
info.opts = opts;

save([opts.saveinIt '.mat'],'bounds','R','Phi','info')
end

if (RHOup - RHOdown < tol),
break;
end

if (timeIter(iter)>=opts.maxTime)
msg(logFile,opts.verbose>0,'\nopts.maxTime reached\n');
break;
end
end

% Post-processing
msg(logFile,opts.verbose>0,'\n> Bounds on the jsr: [%.15g, %.15g]', RHOdown, RHOup);

if(iter==maxiter)
status = 1;
elseif (timeIter(iter)>=opts.maxTime)
status = 2;
else
status = 0;
end

bounds = [RHOdown, RHOup];
elapsedtime = cputime - starttime;

msg(logFile,opts.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.niter = iter;
info.timeIter = timeIter(1:iter);
info.allLb = lowerBounds(1:iter);
info.allUb = upperBounds(1:iter);
info.opts = opts;

% Save output option
if (ischar(opts.saveEnd))
save([opts.saveEnd '.mat'],'bounds','R','Phi','info')
end

% Figures
if(opts.maxrel2D.plotBounds)
figure
it = 1:iter;
plot(it,info.allUb,'-*r',it,info.allLb,'-g+')
title('maxRel2D : Evolution of the bounds on the JSR')
legend('Upper bound','Lower bound')
xlabel('Iterations')
end

if (opts.maxrel2D.plotEllips)
figure
axis equal;
hbar = plot(cos(Phi)./R, sin(Phi)./R, 'r-', 'Linewidth', 3);
legend(hbar, '||x||_{MR} = 1');
end

if (opts.maxrel2D.plotTime)
figure
if (info.timeIter(end) > 600)
time = info.timeIter/60;
plot(time,info.allLb,'-+g',time,info.allUb,'-*r','MarkerSize',5);
xlabel('time from start in min')
else
plot(info.timeIter,info.allLb,'-+g',info.timeIter,info.allUb,'-*r','MarkerSize',5)
xlabel('time from start in s')
end
title('maxRel2D : Evolution of the bounds w.r.t. to time')
legend('Lower bound','Upper bound')
end

%%%%%%%%%%%%%%%%%%%%%%%%%
% AVERAGING FUNCTIONS   %
%%%%%%%%%%%%%%%%%%%%%%%%%

function z = arithMean(x, y)
z = (x+y)/2;

function z = geomMean(x, y)
z = sqrt(x*y);

function z = harmMean(x, y)
z = 2*x*y/(x+y);