# Multicanonical Monte Carlo scheme for finding rare growth factors

### 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
% 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
% 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

```