No BSD License  

Highlights from
Multicanonical Monte Carlo scheme for finding rare growth factors

image thumbnail

Multicanonical Monte Carlo scheme for finding rare growth factors

by

 

05 Oct 2006 (Updated )

A multicanonical Monte Carlo (MMC) scheme developed to approximate the tails of growth factor probab

gf_mmc_driver.m
%  "This is the driver script to be run from the command line. Also 
%  requires functions mcmc1 and mcmc2. This code is referred to by an 
%  upcoming article by T. A. Driscoll and K. Maki in the Education 
%  section of SIAM Review, to appear 2006-2007."
%
%  MMC iteration applied to the growth factors of random matrices.  
%  Runs 'numberchain' instances of MMC simultaneously using the Distributed
%  Computing Toolbox.
%

N = 8;                              % size of nxn matrix
binedge = -.5:1:(2^(N-1)+.5);       % edges of target bins
Nbin = length(binedge) - 1;         % number of bins
M = 5e5;                            % number of realizations per iteration
Niter = 70;                         % number of iterations
burnin = 5e5;                       % burn-in period for MCMC
numberchain = 12;                   % number of independent MMC runs

% ***** Initialization *****

randn('state',sum(100*clock)), rand('state', sum(100*clock)); 

Accept = zeros(Niter,numberchain);  % acceptance rate per iteration per run
sigma = (1/N)*ones(numberchain,1);  % sigma parameter per run
H = zeros(Nbin,Niter,numberchain);  % hits/bin per iteration per run
P = ones(Nbin,numberchain)/Nbin;    % computed histogram per run
g = zeros(Nbin-1,Niter,numberchain);% parameter needed for Berg's iteration
Pnew = zeros(Nbin,1);               % updated inverse weights

GFinitial = zeros(numberchain,1);   % computed initial growth fac for MC
X = zeros(N,N,numberchain);         % initial matrices for MC runs
chaingrowth = 10*(1:numberchain-2); % desired initial growth fac for MC

% Establishing communication for distributed computing
jm = findResource('jobmanager','Name','woprdce');  

% Calculating one initial matrix for Markov Chain
X(:,:,1) = randn(N,N);


% ***** Start MMC iteration *****

for i=1:Niter
    
    % Calculating remaining initial matrices for Markov Chain
    X(:,:,2)=randn(N);
    
    % Implementing Higham and Higham's Theorem
    for z=3:numberchain
        while( GFiter(z)>chaingrowth(z-2)+.5 || GFiter(z) < chaingrowth(z-2)-.5 )
            Mstart = toeplitz([1 -.99999*(-1+chaingrowth(z-2)^(1/(N-1)))*ones(1,N-1)], [1 zeros(1,N-1)]);
            T = triu( (2*rand(N-1)-1)/N );
            T = [T;zeros(1,N-1)];
            Tlc = (chaingrowth(z-2).^((0:N-1)/(N-1)));
            D = diag( sign(randn(N,1)) );
            X(:,:,z) = D*Mstart*[T Tlc'];
            [L,U]=lu(X(:,:,z));
            GFiter(z) = (max(abs(U(:))))/(max(max(abs(X(:,:,z)))));
        end
    end

    % ***** Running MCMC for burn-in *****
    
    % Telling the Distributed Computing Toolbox to complete one job with 
    % 'numberchain' tasks.  Each task is comprised of running a MCMC 
    % for the burnin time (mcmc1.m) with a different initial matrix.

    j = createJob(jm);
    set(j,'Timeout',1200);
    set(j,'FileDependencies',{ 'mcmc1.m'});
    for z = 1:numberchain
        createTask( j, @mcmc1, 2, {X(:,:,z), burnin, sigma(z), P(:,z), binedge});
    end
    submit(j);
    waitForState(j,'finished');
    data = getAllOutputArguments(j);
    
    % Received computed data and recording results 
    for z = 1:numberchain
        Accept(i,z) = data{z,1};
        X(:,:,z) = data{z,2};
    end
    clear j

    % Adjusting the parameter sigma based on burnin results
    for z = 1:numberchain
        if( Accept(i,z)/burnin < .1 )
            sigma(z) = max( sigma(z)/2, eps);
        elseif( Accept(i,z)/burnin > .4 )
            sigma(z) = 2*sigma(z);
        else
            sigma(z) = sigma(z);
        end
    end

    % ***** Running MCMC for M *****
    
    % Again telling the Distributed Computing Toolbox to complete one job 
    % composed of 'numberchain' tasks.  Each task continues running the 
    % independent MCMC iterations (script mcmc2.m) from above.  Now
    % statistics are calculated.
    
    j = createJob(jm);
    set(j,'Timeout',1200);
    set(j,'FileDependencies',{'mcmc2.m'});
    for z=1:numberchain
        createTask(j,@mcmc2, 3, {X(:,:,z), M, sigma(z), P(:,z), binedge});
    end
    submit(j);
    waitForState(j,'finished');
    data = getAllOutputArguments(j);
    
    % Received computed data and recording results
    for z = 1:numberchain
        Accept(i,z) = data{z,1};
        H(:,i,z) = data{z,2};
        X(:,:,z) = data{z,3};
    end
    clear j

    % Tuning the parameter sigma again
    for z = 1:numberchain
        if( Accept(i,z)/(M) < .1 )
            sigma(z) = max( sigma(z)/2, eps);
        elseif( Accept(i,z)/(M) > .4 )
            sigma(z) = 2*sigma(z);
        else
            sigma(z) = sigma(z);
        end
    end
    
    % *****Berg update*****
    for z =1:numberchain
        Pnew = P(:,z);
   
        for k = 1:(Nbin-1)
            if (H(k+1,i,z)*H(k,i,z) > 0)
                g(k,i,z) = H(k+1,i,z)*H(k,i,z)/(H(k+1,i,z)+H(k,i,z));
                ghat(k) = g(k,i,z)/(sum(g(k,1:i,z)));
                Pnew(k+1) = Pnew(k)*(P(k+1,z)/P(k,z))*((H(k+1,i,z)/H(k,i,z))^ghat(k));
            else
                Pnew(k+1) = Pnew(k)*(P(k+1,z)/P(k,z));
            end
        end
        P(:,z) = Pnew/sum(Pnew);
    end    
end




Contact us