% The following MATLAB .m file finds the value of the second highest
% eigenvalue in absolute value of the matrix P from the 'frog example'.
% A lower second highest eigenvalue generally means the matrix converges
% more quickly to stationarity.
m = 9; % P(x,y) = 1/(2m+1) for all x, y within m of each other
nodes = 200; % P is a nodes*nodes matrix
beta = 10; % target distribution is proportional to (dist(x,1)+1)^-b
% where dist(x,1) is the number of nodes to get from x to 1
% the following two loops calculate the stationary (target) distribution
% of the matrix P
denom = 0;
for i = 1:nodes
min = i-1;
if nodes-i+1 < min
min = nodes-i+1; % distance from the current node to 1
end
denom = denom + (min + 1)^(-beta);
end
for i = 1:nodes
min = i-1;
if nodes-i+1 < min
min = nodes-i+1;
end
pi(1, i) = (min + 1)^(-beta)/denom; % stationary distribution of P
end
for i = 1:nodes
P(i, i) = 1/(2*m+1);
for j = 1:nodes
if pi(1, j) >= pi(1, i)
alpha = 1; % alpha is the probability that we jump from the
% current node, i, to the proposed node, j
else
alpha = pi(1, j)/pi(1, i);
end
if 0 < abs(i-j) & abs(i-j) < (m+1)
P(i, j) = alpha/(2*m+1);
P(i, i) = P(i, i) + (1-alpha)/(2*m+1);
elseif abs(i-j) >= (nodes-m) % the highest nodes and the lowest
% nodes have edges between them too
P(i, j) = alpha/(2*m+1);
P(i, i) = P(i, i) + (1-alpha)/(2*m+1);
elseif (m+1) <= abs(i-j) & abs(i-j) < (nodes-m)
P(i, j) = 0;
end
end
end
y = eig(P); % finds a vector of the eigenvalues; the following loop
% isolates the second highest one in absolute value
max = 0;
max2 = 0;
for k = 1:nodes
if abs(y(k, 1)) > max
max2 = max;
max = abs(y(k, 1));
elseif abs(y(k, 1))>max2 & max>= abs(y(k, 1))
max2 = abs(y(k, 1));
end
end
max2
% the following can be added on to the above program to find bounds
% on variation distance using all the eigenvalues; a slight modification
% can be made to see bounds on L2(pi) distance
iters = 40;
Q = P'; % we want left eigenvectors of P which will be right eigenvectors
% of Q
sum = 0;
[evecs, evals] = eig(Q);
start = zeros(nodes, 1);
start(1,1) = 1;
am = evecs\start; % finding a sub m for all m
for i = 1:nodes
for j = 2:nodes
sum = sum + abs(am(j,1)*evecs(i,j)) * abs(evals(j,j))^iters;
end
end
sum = sum/2;
sum % bound on variation distance