Code covered by the BSD License  

Highlights from
The JSR toolbox

image thumbnail

The JSR toolbox

by

 

10 Oct 2011 (Updated )

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

jsr_norm_balancedComplexPolytope(M, C, k, myopts)
function [bounds, V, info] = jsr_norm_balancedComplexPolytope(M, C, k, myopts)

% JSR_NORM_BALANCEDCOMPLEXPOLYTOPE Approximates the jsr using b.c.p.'s.
%
%    [BOUNDS, V, INFO] = JSR_NORM_BALANCEDCOMPLEXPOLYTOPE(M, C, K)
%      returns lower and upper bounds on the jsr of M using balanced complex
%      polytope norms.
%      The set M must be a cell array containing the considered matrices.
%      The matrix C is the candidate optimal product.
%      The parameter K corresponds to the length of the candidate product. 
%
%    [BOUNDS, V, INFO] = JSR_NORM_BALANCEDCOMPLEXPOLYTOPE(M, C, K, OPTS)
%      Uses the values of the parameters defined in the structure OPTS.
%      See below and JSRSETTINGS for the available parameters and options.
%
%      BOUNDS are the best lower and upper bounds found.
%
%      V is a matrix from wich the rows correspond to an essential 
%      system of vertices for the polytope.
%
%      INFO is a structure containing various data about the iterations :
%         info.status         - = 0 if normal termination, 1 if maxiter has
%                               been reached, 2 if stopped in iteration
%                               (CTRL-C or opts.maxTime reached)
%         info.elapsedtime 
%         info.niter          - number of iterations performed
%         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.Vhist          - Cell containing the sets of vectors at each
%                               iteration 
%         info.popHist        - Number of vectors in essential set 
%                               at each iteration
%         info.allLb          - lower bounds, (constant = rho(C)^(1/K) )
%                               [1xniter double] 
%         info.allUb          - evolution of the upper bound
%                               [1xniter double] 
%         info.opts           - the structure of options used
%
%  The field opts.bcp (generated by jsrsettings) can be used to
%  tune the method :
%
%      bcp.maxiter        - max number of iterations, (500)
%
%      bcp.plotBounds     - if 1 plots the evolution of the bounds, (0)
%
%      bcp.plotpopHist    - if 1 plots evolution of the size of the
%                           essential set, (0)
%      bcp.plotTime       - if 1 plots evolution of the bounds w.r.t. time 
%                           of computation, (0)
%
% See also JSRSETTINGS
%
% REFERENCES
%    N.Guglielmi and M.Zennaro, 
%      "Finding extremal complex polytope norms for
%       families of real matrices",
%      SIAM J. Matrix Anal. Appl. 31(2):602-620, 2009
%    V.Protasov,
%      "The generalized spectral radius. A geometric approach",
%      Izvestiya Mathematics, 61(5):995-1030, 1997


if (nargin<4)
    myopts = jsrsettings;
    if (nargin<3)
        error('This method requires a candidate optimal product and its length as input.')
    end
end

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

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

msg(logFile,myopts.verbose>1,'\n \n******** Starting jsr_norm_balancedComplexPolytope ******** \n \n')
starttime = cputime;

% Initialization
opts.disp = 0;
options = optimset('Display', 'off', 'LargeScale', 'off','Algorithm','active-set','TolFun', 1e-11);
m = length(M);
N = size(M{1}, 1);
[v, d] = eigs(C, 1, 'LM', opts);
alpha = abs(d)^(1/k);
F = cellDivide(M,alpha);
V = [v.' ; v'];
X = [v.' ; v'];
X = [X ; -X];
w = size(X, 1);
maxiter = myopts.bcp.maxiter;
mumin = Inf;
Vlog = cell(maxiter, 1);
allUb = zeros(1,maxiter);
nv = zeros(1,maxiter);
status = 2;
timeIter = zeros(1,maxiter);

% Iteration
for iter = 1:maxiter,
    Xsize = 4;
    Xused = 0;
    Xnew = zeros(Xsize, N);
    mu = 1;
    for i = 1:m,
        Y = F{i}*X.';
        for j = 1:w,
            objective = 2*eye(2*w);
            constraintAeq = [real(X.'), -imag(X.') ; imag(X.'), real(X.')];
            constraintBeq = [real(Y(:, j)) ; imag(Y(:, j))];
            [dummy, mutemp, exitflag] = quadprog(objective, zeros(2*w, 1), [], [], constraintAeq, constraintBeq, zeros(2*w, 1), Inf*ones(2*w, 1), [], options);
            if (exitflag == -2) || (exitflag == 0) || isempty(mutemp),
                mutemp = Inf;
            end
            mu = max(mu, mutemp);
            if (mutemp > 1),
                if (Xused == Xsize),
                    Xsize = Xsize*2;
                    Xnew(Xsize, N) = 0;
                end
                Xused = Xused + 1;
                Xnew(Xused, :) = Y(:, j).';
            end
        end
    end
    mumin = min(mumin, mu);
    
    if (iter>1000)
        if (mod(iter,100)==0)
            msg(logFile,myopts.verbose>0,'> Iteration %d - current bounds: [%.15g, %.15g]', iter, alpha, alpha*mumin);
        end
    elseif (iter>150)
        if (mod(iter,50)==0)
            msg(logFile,myopts.verbose>0,'> Iteration %d - current bounds: [%.15g, %.15g]', iter, alpha, alpha*mumin);
        end
    elseif (iter>20)
        if (mod(iter,10)==0)
            msg(logFile,myopts.verbose>0,'> Iteration %d - current bounds: [%.15g, %.15g]', iter, alpha, alpha*mumin);
        end
    else
        msg(logFile,myopts.verbose>0,'> Iteration %d - current bounds: [%.15g, %.15g]', iter, alpha, alpha*mumin);
    end
    
    allUb(iter) = alpha*mumin;
    Xnew = Xnew(1:Xused, :);
    if isempty(Xnew),
        timeIter(iter) = cputime-starttime;
        break;

    end
    X = unique(Xnew, 'rows');
    V = unique([V ; X], 'rows');
    if (size(V, 1) > size(V, 2))
        i = 1;
        w = size(V, 1)-1;
        objective = 2*eye(2*w);
        while (i <= w+1),
            VV = V;
            VV(i, :) = [];
            constraintAeq = [real(VV.'), -imag(VV.') ; imag(VV.'), real(VV.')];
            constraintBeq = [real(V(i, :).') ; imag(V(i, :).')];
            [dummy, mutemp, exitflag] = quadprog(objective, zeros(2*w, 1), [], [], constraintAeq, constraintBeq, zeros(2*w, 1), Inf*ones(2*w, 1), [], options);
            if (exitflag == -2) || (exitflag == 0) || isempty(mutemp),
                mutemp = Inf;
            end
            if (mutemp <= 1),
                V(i, :) = [];
                w = w-1;
                objective = 2*eye(2*w);
            else
                i = i+1;
            end
        end
    end
    Vlog{iter} = V;
    nv(iter) = size(V,1);
    X = intersect(X, V, 'rows');
    w = size(X, 1);
    
    timeIter(iter) = cputime-starttime;
    
    % Save to file option
    if (ischar(myopts.saveinIt))
        if (iter==maxiter);status=1;end
        bounds = [alpha, alpha*mumin];
        elapsedtime = cputime - starttime;

        info.status = status;
        info.elapsedtime = elapsedtime;
        info.niter = iter;
        info.timeIter = timeIter(1:iter);
        info.Vhist = Vlog(1:iter);
        info.popHist = nv(1:iter);
        info.allLb = ones(1,iter)*alpha;
        info.allUb = allUb(1:iter);
        info.opts = myopts;

        save([myopts.saveinIt '.mat'],'bounds','V','info')
    end
    
    if (timeIter(iter)>= myopts.maxTime)
        msg(logFile,myopts.verbose>0,'\nopts.maxTime reached\n');
        break;
    end
end

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

% Post-processing
msg(logFile,myopts.verbose>0,'\n> Bounds on the jsr: [%.15g, %.15g]', alpha, alpha*mumin);
bounds = [alpha, alpha*mumin];
elapsedtime = cputime - starttime;
msg(logFile,myopts.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.Vhist = Vlog(1:indLastV);
info.popHist = nv(1:indLastV);
info.allLb = ones(1,iter)*alpha;
info.allUb = allUb(1:iter);
info.opts = myopts;

% Save output option
if (ischar(myopts.saveEnd))
   save([myopts.saveEnd '.mat'],'bounds','V','info')
end

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

if (myopts.bcp.plotpopHist)
   figure
   plot(1:indLastV,nv(1:indLastV),'-*');
   title('bcp : Evolution of size of essential system');
   xlabel('Iterations')
   xlim([1 indLastV])
end

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

end

Contact us