Equilibrium Networks

In an earlier post we covered some more sophisticated ways of implementing a graded lexicographic index, a/k/a the state index. One of these dealt with producing a superset of the state indices of the neighbors of a given state (here is abusively interpreted as a graph embedded in where edges join any two vertices that are a distance apart). A careful exploitation of this technique leads to a very efficient way of updating the state index for random walks on Basically, the idea is to just determine which element of the neighborhood is the correct one, and restrict the neighborhood algorithm to a single vertex. In practice this means a lot of careful indexing and addition of appropriate binomial coefficients, but the basic idea is really just to omit the extra work in the neighborhood algorithm.

Here’s some MATLAB code that does the trick:

function y = stateindexupdate(a);

% a is an array of states (assumed to be) in a bucketspace
% this illustrates efficient stateindex updating

T = size(a,1);
B = size(a,2);
b = sum(a(1,:));

% Pascal matrix
mb = max(b+1,B+1);
P1 = pascal(mb);
P = P1(1:B-1,:);
PA = [zeros(1,mb);P];

ia(1) = stateindex(a(1,:));

for j = 2:T
 % compute tau:
 tau = 1 + cumsum(fliplr(a(j,:)));
 % compute the ball move that got us here:
 move = a(j,:) - a(j-1,:);
 n = 0; % init the "slot" index
 % compute the corresponding slot index...
 % first, ID the first nonzero element of the move a la Java/C:
 for k = 1:B
  if move(k)
   % now find the second nonzero element m:
   for i = k+1:B
    if move(i)
     m = i;
     break;
    end
   end
   % the slot n is as follows:
   n = move(k) * ( (B-k)*((B-k)+1)/2 - B + m);   % move(i) = sign(move(i))
   break;
  end
 end

 % compute the first (signed) triangular number further from zero:
 for k = 0:B
  n2 = k * (k+1)/2;
  if n2 >= abs(n)
   n2 = sign(n) * n2;
   break;
  end
 end

 d = abs(n2-n);
 l = ( -1 + sqrt(1+8*abs(n2)) )/2;  % l*(l+1)/2 = abs(n2);

 % the sum will go from d+2 to l+1...
 ind = (d+2):(l+1);
 % finally, perform the sum
 lind = length(ind);
 temp = zeros(1,lind);
 if n2 < 0
  for i = 1:lind
   temp(i) = PA(ind(i),tau(ind(i)-1)-1);
  end
  delta = -sum(temp);
 elseif n2 > 0
  for i = 1:lind
   temp(i) = PA(ind(i),tau(ind(i)-1));
  end
  delta = sum(temp);
 else
  delta = 0;
 end

 % finally, update the stateindex:
 ia(j) = ia(j-1) + delta;
end

y=ia;

This runs really fast. A Java version of this code that we produced performs many millions of updates per second on a decent processor core.

One of the nicest things about “grlex” ordering and indexing is that it can be used as a primitive to tackle lots of otherwise difficult problems. For instance, I’ve used it to help (among other things) emulate an older network security monitoring platform, analyze generalized multivariate de Moivre martingales, perform multivariate Lagrange interpolation, study energy functions on equivalence classes of directed graphs, and to construct propagation operators implementing periodic hexagonal, truncated-octahedral, and more generally permutohedral boundary conditions for lattice models (such as spin or lattice Boltzmann models) with surprising ease.

Hopefully some of these more advanced techniques for grlex stuff will be of use elsewhere. If you’d like to use some of this code, just ask: we’ll provide it under an open-source license.