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

Kara Maki (view profile)

 

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