Equilibrium Networks

While browsing through my notes recently I came across this cute and not-quite trivial, but entirely elementary result. Since it doesn’t seem to be available anywhere else, I present it as blog fodder. (As usual, I will offer $.50 and a Sprite to the first person who provides a usable reference in the event that this result is already known.)

Consider the second-order matrix difference equation

with given. Some algebra shows that (e.g.)

etc.

where the bosonic symmetrized product for two “species” of order is defined as the sum of all distinguishable products with occurences of the first species and of the second (if or then the BSP is a pure power of the nontrivial species).

For example, equals

It is not hard to verify that

and more generally for that

The same algebra that facilitated writing down the first few terms for the difference equation quickly leads to the general solution, viz.

for odd, and

for even.

I used this to understand what happens when one of the internal states of a classical Bose gas becomes “more fermionic” in the sense that fewer particles can occupy that state, but the microscopic transition probabilities are otherwise unchanged. (The underlying motivation was to apply this to get a handle on the effects of some complex configurations for an older mathematically-oriented network monitoring framework in a tractable traffic regime.) It turns out that this is completely uninteresting unless the remaining “bosonic” states are close to equally likely to occur: in this event you get a slight but also quite counter-intuitive effect on the state distributions. Pseudocolor figures of the energy functions corresponding to these distributions will illustrate:

3 totally bosonic interal states. The distribution is multigeometric, and the energy is affine.

Capping the occupancy of one internal state

A more restrictive cap

A "fermionic" internal state

Note that the color scales differ for each figure. I think the pictures indicate how surprising the behavior of this (completely classical) “poor man’s exclusion principle” is: the effect is as unexpected as it is slight.

Hopefully posting the main result will save some other folks a head-scratch or two. I’d be interested to know if it’s applied elsewhere.

The previous post on this topic covered basic code for producing a graded lexicographic index (a/k/a the state index) and its inverse. This post will cover some more advanced techniques. The first one of these is simple: as de Boor pointed out (and it took me quite a while to catch on) the state index of can be given by a simple formula. Let denote the matrix given by the following MATLAB code:

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

For example, with this is

(Identifying the elements of the Pascal matrix with binomial coefficients would introduce some minor technical complications w/r/t the top row, and so we avoid this.) Now let

Then

where

and an empty sum is evaluated as zero. This is embodied in the following MATLAB code:

J = find(alpha);
K = length(J);

for k = 1:K
   temp1 = sum(alpha(J(1:(k-1))));
   for m = 1:alpha(J(k))
      stateindex = stateindex + PA(B+1-J(k),b+2-temp1-m);
   end
end

stateindex = stateindex + 1;

This is hard to beat in general. However, for some applications it is sufficient to compute the state index at a neighboring site. And this can actually be done much more effectively in practice.

A deceptively simple algorithm always produces a superset of the neighborhood of a state in bucketspace. Write

and for a finite subset of let and respectively denote the least and greatest elements of Now set and given for form

where

and

Then the state indices of the nearest neighbors of are contained in the set It is not hard to see that this algorithm requires simple operations (e.g., lookup, addition/subtraction, min/max) if we store the (fairly small) Pascal matrix. On “interior” states (i.e., states with no zero entries), this algorithm produces the exact neighborhood. However, on “boundary” states (i.e., states with at least one zero entry), some spurious indices are returned. It is not obvious (especially given its structure) how the algorithm might be refined to deal with this problem directly.

However, the spurious indices can themselves be checked efficiently. Given and the following MATLAB code illustrates how to invert the state index:

PA2 = flipud(fliplr(PA));
alpha = zeros(1,B);
ind0 = 0; % init
for j = 1:B
   temp1 = cumsum(PA2(j,(ind0+1):(mb-1)));
   ind1 = max(find(temp1<s));
   if ind1
      temp2 = temp1(ind1);
      s = s - temp2;
      alpha(j) = ind1;
      ind0 = ind0 + ind1;
   end
end

The inversion essentially takes operations. With this tool in hand, we can trivially establish which state indices are spurious and which are not.

In the next post on this topic, we’ll produce a technique which is actually much faster than any of these for updating state indices in both MATLAB and Java, even after accounting for considerable care in the optimization (such as accounting for the order in which bitshifts are executed).